* fichier : xfem01.dgibi ************************************************************************ * Section : Mecanique Endommagement ************************************************************************ ******************************************************************* * * xfem01.dgibi * * calcul avec elements xfem (XQ4R) * d'une plaque elastique en traction * avec fissure inclinée * * création : bp, le 24.09.2009 * ajout g_theta : bp, 08.06.2010 * solution de reference pour g,k1 et k2 donne par : * [TADA, STRESS ANALYSIS HANDBOOK,1973] * [ISIDA,IJF(7),301-316,1971] -> facteur de forme * propagation elementaire ******************************************************************** opti echo 1 ; *** Traces *graph = faux; *graph = vrai; graph = vrai; opti trac PSC ; opti EPTR 5; *** finesse du maillage (=flfin) et cas test complet (=fllong) *flfin = faux; flfin = vrai; fllong = faux; *fllong = vrai; ******************************************************* *** Options de calcul : 2D opti dime 2 elem qua4 mode plan defo; ******************************************************* *** Maillage : Plaque rectangulaire saine * dimension et finesse *La1 = 0.1 ; La1 = 0.15 ; Lb1 = 0.5 ; na1 = 20; si(flfin); La1 = 0.2 ; na1 = 40; *na1 = 100; * na1 = 200; * na1 = 400; fins; dx1 = (2.*La1) / (flot na1) ; nb1 = na1 / 2; dx2 = (2.*Lb1) / (flot nb1) ; * zone centrale DENS (dx1) ; pa1 = (-1.*La1) (-1.*La1) ; pa2 = ( La1) (-1.*La1) ; pa3 = ( La1) ( La1) ; pa4 = (-1.*La1) ( La1) ; da1 = pa1 droi na1 pa2; da2 = pa2 droi na1 pa3; da3 = pa3 droi na1 pa4; da4 = pa4 droi na1 pa1; sa1 = dall da1 da2 da3 da4; * zone exterieure (pour avoir un milieu quasi-infini) nb1 = na1 / 2; dx2 = (2.*Lb1) / (flot nb1) ; DENS (dx2) ; pb1 = (-1.*Lb1) (-1.*Lb1) ; pb2 = ( Lb1) (-1.*Lb1) ; pb3 = ( Lb1) ( Lb1) ; pb4 = (-1.*Lb1) ( Lb1) ; pb11 = (-1.*La1) (-1.*Lb1) ; pb12 = ( 1.*La1) (-1.*Lb1) ; pb21 = ( Lb1) (-1.*La1) ; pb22 = ( Lb1) ( 1.*La1) ; pb31 = ( La1) ( Lb1) ; pb32 = (-1.*La1) ( Lb1) ; pb41 = (-1.*Lb1) ( La1) ; pb42 = (-1.*Lb1) (-1.*La1) ; db1 = pb11 droi na1 pb12; db2 = pb21 droi na1 pb22; db3 = pb31 droi na1 pb32; db4 = pb41 droi na1 pb42; sa1b1 = (regl da1 db1 'DINI' dx1 'DFIN' dx2); sa2b2 = (regl da2 db2 'DINI' dx1 'DFIN' dx2); sa3b3 = (regl da3 db3 'DINI' dx1 'DFIN' dx2); sa4b4 = (regl da4 db4 'DINI' dx1 'DFIN' dx2); sabtot= sa1b1 et sa2b2 et sa3b3 et sa4b4; dc1 = (cote sa1b1 4) coul bleu; dd1 = (cote sa1b1 2) coul vert; dc2 = (cote sa2b2 4) coul bleu; dd2 = (cote sa2b2 2) coul vert; dc3 = (cote sa3b3 4) coul bleu; dd3 = (cote sa3b3 2) coul vert; dc4 = (cote sa4b4 4) coul bleu; dd4 = (cote sa4b4 2) coul vert; *1er coin dc2ext = dc2 plus (pb2 moin pb21); dd1ext = dd1 plus (pb2 moin pb12); dsb2 = dc2ext et dd1ext et dc2 et dd1; elim dsb2 (1.E-4*dx1); sb2 = dall dc2ext dd1ext dc2 dd1; *2eme coin dc3ext = dc3 plus (pb3 moin pb31); dd2ext = dd2 plus (pb3 moin pb22); dsb3 = dc3ext et dd2ext et dc3 et dd2; elim dsb3 (1.E-4*dx1); sb3 = dall dc3ext dd2ext dc3 dd2; *3eme coin dc4ext = dc4 plus (pb4 moin pb41); dd3ext = dd3 plus (pb4 moin pb32); dsb4 = dc4ext et dd3ext et dc4 et dd3; elim dsb4 (1.E-4*dx1); sb4 = dall dc4ext dd3ext dc4 dd3; *4eme coin dc1ext = dc1 plus (pb1 moin pb11); dd4ext = dd4 plus (pb1 moin pb42); dsb1 = dc1ext et dd4ext et dc1 et dd4; elim dsb1 (1.E-4*dx1); sb1 = dall dc1ext dd4ext dc1 dd4; sbtot = sb1 et sb2 et sb3 et sb4; *petite modif pour propager + loin si(flfin); sa1 = sa1 et sa2b2 et sa4b4; sb1 = sa1b1 et sa3b3 et sbtot; sino; * attention on reutilise sb1 !!! sb1 = sabtot et sbtot ; fins; stot1 = sa1 et sb1; mess 'maillage central ext total' ; mess 'nbno= ' (nbno sa1) (nbno sb1) (nbno stot1); mess 'nbel= ' (nbel sa1) (nbel sb1) (nbel stot1); *opti donn 5; ******************************************************* *** Maillage de la fissure * longueur et angle de la fissure *a0 = 0.1; a0 = 0.05; *alpha0 = 0.; alpha0 = 20.; *alpha0 = 40.; *alpha0 = 45.; *alpha0 = 50.; * finesse du maillage de la fissure *na0 = 10; *da0 = a0 / (flot na0) ; da0 = dx1 / 10.; na0 = enti (a0 / da0); * maillage DENS (da0) ; ptip1 = (a0 * (cos alpha0)) (a0 * (sin alpha0)); ptip2 = (-1.*a0 * (cos alpha0)) (-1.*a0 * (sin alpha0)); lcrack0 = (ptip2 droi na0 ptip1) coul 'ROUG'; * table (+ pratique pour gérer la propa) TXFEM = tabl; TXFEM . 0 = tabl; TXFEM . 0 . 'POINTE1' = ptip1; TXFEM . 0 . 'POINTE2' = ptip2; TXFEM . 0 . 'FISSURE' = lcrack0; * creation des level set psi0 phi0 = PSIPHI sa1 lcrack0 'DEUX' ptip1 ptip2; isolv7 = prog (-2.*a0) PAS (a0 / 4.) (2.*a0); si(graph); *opti isov lign; TRAC (stot1 et lcrack0) ; TRAC phi0 (stot1 et lcrack0) isolv7 ; TRAC psi0 (stot1 et lcrack0) isolv7; *opti isov surf; finsi; *opti donn 5 ; ******************************************************* *** Modeles et materiaux * definition des models elementaires mod1a = MODE sa1 mecanique elastique isotrope xq4R 'CONSTITUANT' 'A'; mod1b = MODE sb1 mecanique elastique isotrope 'CONSTITUANT' 'B'; * apres cela on assemble mod1tot = mod1a et mod1b; * enrichissement Che1X = tabl; Che1X . 0 = TRIE mod1tot psi0 phi0; Rel1X = tabl; * constructionsion des blocages des ddl X-fem non actifs dans * les éléments de transition. Rel1X . 0 = RELA mod1tot; * on crée le ou les materiau et on les assemble Ey1 = 2.E11 ; nu1 = 0.3; rho1 = 7800.; mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1; ******************************************************* *** Matrices + CL Rig1tot = RIGI mod1tot mat1tot; * on bloque les depl de corps rigides pb12 = sb1 POIN 'PROC' (0.5 * (pb1 PLUS pb2)); pb23 = sb1 POIN 'PROC' (0.5 * (pb2 PLUS pb3)); pb34 = sb1 POIN 'PROC' (0.5 * (pb3 PLUS pb4)); pb41 = sb1 POIN 'PROC' (0.5 * (pb4 PLUS pb1)); cl1x = BLOQ 'UX' (pb12 et pb34); cl1y = BLOQ 'UY' (pb41 et pb23); cl1tot = cl1x et cl1y; K1tot = Rig1tot et cl1tot et (Rel1X . 0); syy0 = 1.E9; dbup = dc4ext et db3 et dd2ext; dbdo = dd4ext et db1 et dc2ext; pre1 = PRES 'MASSIF' mod1tot (dbup et dbdo) (-1.*syy0); si(graph); TRAC (VECT pre1 'FORC' 'BLEU') (stot1 et lcrack0); finsi; ******************************************************* *** RESOlution avec RESO u0 = RESO K1tot pre1; ******************************************************* *** Post traitement * verif du fonctionnement des principaux operateurs defel0 = EPSI mod1tot u0; sig0 = SIGM mod1tot u0 mat1tot; svm0 = VMIS sig0 mod1tot; gru0 = GRAD mod1tot u0; def0 = ELAS mod1tot sig0 mat1tot; ddef0 = def0 - defel0; * calcul du residu Fint0 = BSIG sig0 mod1tot; Freac0 = REAC cl1tot u0; Res0 = Fint0 - (Freac0 + pre1); * * calcul energie elastique = energie exterieure * wdefel0= 0.5 * (ENER sig0 defel0 mod1tot); * Iwdef0 = INTG wdefel0 mod1tot; * wext = 0.5 * (PSCA u0 pre1 (mots UY) (mots FY)); * Iwext = maxi (RESU wext); * * travail des efforts internes de deformations * wdef0 = WORK mod1tot sig0 gru0; * messages mess 'erreur relative / Residu =' ((maxi Res0) / (maxi pre1)); * mess 'erreur relative / Energie =' ((Iwdef0-Iwext) / (abs Iwext)); mess 'ecart relatif EPSI / SIGM+ELAS= ' (maxi ddef0 abs) ; * traces u0phy = XFEM 'RECO' u0 mod1tot ; ucrk0u ucrk0d = XFEM 'FISS' lcrack0 u0 mod1tot ; si(graph); def0 = DEFO u0phy stot1 20.; TRAC sig0 mod1tot def0; def0 = def0 et (DEFO ucrk0u lcrack0 20.) et (DEFO ucrk0d lcrack0 20.); TRAC def0; * trac defel0 mod1tot (lcrack0); * TRAC svm0 mod1tot (stot1 et lcrack0); * TRAC wdefel0 mod1tot (stot1 et lcrack0); * TRAC gru0 mod1tot (stot1 et lcrack0); * test de l'option CHAN CHPO SUPP : sig0_p = CHAN 'CHPO' sig0 mod1a 'SUPP'; geo1 = EXTR sig0_p 'MAILLAGE' ; TRAC (EXCO sig0_p 'SMYY') geo1 (sa1 COUL gris); TRAC (EXCO sig0 'SMYY') mod1a; finsi; * OPTI DONN 5 echo 1 trac X; ******************************************************* * calcul de l integrale J, de K1, de K2 *** valeur theorique de J,K1,K2 *** k0th1 = syy0 * ((a0 * pi) ** 0.5) * (cos alpha0); Estar1 = (1. - (nu1**2)) / Ey1; Gth1 = Estar1 * (k0th1 ** 2); * facteur de forme pour a/b=0.1 et h/b=Inf *facfor1 = 1.006; * facteur de forme pour a/b=0.1 et h/b=1 facfor1 = 1.014; Gth1 = (facfor1 ** 2) * Gth1; K1ref = k0th1 * (cos alpha0) * facfor1; K2ref = k0th1 * (sin alpha0) * facfor1; *** G_THETA -> J *** SUPTAB = TABL ; SUPTAB.'OBJECTIF' = MOT 'J'; SUPTAB.'PSI' = psi0; SUPTAB.'PHI' = phi0; SUPTAB.'FRONT_FISSURE' = ptip1 ; SUPTAB.'MODELE' = mod1tot; SUPTAB.'CARACTERISTIQUES' = mat1tot; SUPTAB.'SOLUTION_RESO' = u0; SUPTAB.'CHARGEMENTS_MECANIQUES' = pre1; *on teste toutes les options si(fllong); SUPTAB.'COUCHE' = 0; G_THETA SUPTAB; * on teste l'option ou l on choisit soi meme le champ theta q7 = (SUPTAB . 'CHAMP_THET1') ; si(graph); vq7 = VECT q7 BLEU ; trac vq7 (stot1 et lcrack0); finsi; OTER SUPTAB 'COUCHE'; SUPTAB . 'CHAMP_THETA' = q7; G_THETA SUPTAB; * on calcule J pour differents contours prj1 = prog (SUPTAB . 'RESULTATS'); prcou1= prog 0.; si(flfin); si (na1 > 10); lecou1 = lect 1 pas 1 15; sino; lecou1 = lect 1 pas 1 (na1/4); fins; sino; lecou1 = lect 1 pas 1 3; fins; ncou1 = dime lecou1; jmoy1 = 0.; repe BCOU1 ncou1; icou1 = extr lecou1 (&BCOU1); SUPTAB.'COUCHE' = icou1; G_THETA SUPTAB; prj1 = prj1 et (prog (SUPTAB . 'RESULTATS')); prcou1=prcou1 et (prog (flot icou1)); jmoy1=jmoy1 + (SUPTAB . 'RESULTATS'); fin BCOU1; evjcou1 = evol BLEU manu 'CONTOUR' prcou1 'J' prj1 ; si(graph); dess evjcou1; fins; * extraction d un J moyen jmoy1 = jmoy1 / ncou1; *cas rapide (on ne teste poas toutes les options) sino; si(flfin); ncou1 = 3; sino; ncou1 = 1; fins; SUPTAB.'COUCHE' =ncou1; G_THETA SUPTAB; jmoy1 = SUPTAB . 'RESULTATS'; fins; *erreur relative / J errG1 = abs (1. - (jmoy1 / Gth1)); *** Calcul des FIC par methode integrale G_THETA *** KTAB = TABL; KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE'; KTAB .'PSI' = psi0; KTAB .'PHI' = phi0; KTAB .'FRONT_FISSURE' = ptip1 ; KTAB .'MODELE' = mod1tot; KTAB .'CARACTERISTIQUES' = mat1tot; KTAB .'SOLUTION_RESO' = u0; KTAB .'CHARGEMENTS_MECANIQUES' = pre1; ncou1 = 1; KTAB . 'COUCHE' = ncou1; G_THETA KTAB; K1G = KTAB . 'RESULTATS' . 'I'; K2G = KTAB . 'RESULTATS' . 'II'; errK1G = ABS (1 - (ABS ( K1G/K1ref))); errK2G = ABS (1 - (ABS ( K2G/K2ref))); GK1K2 = Estar1 * ((K1G**2) + (K2G**2)) ; ******************************************************* *** erreur? *** mess '------------------------------------------------'; mess (chai ' G K1 K2 '); mess (chai 'reference ' Gth1 ' ' K1ref ' ' K2ref ); mess (chai ' G_THETA ' jmoy1 ' ' K1G ' ' K2G ); mess (chai ' ERR_REL ' errG1 ' ' errK1G ' ' errK2G ); *OPTI DONN 5 echo 1 trac X; * Si (errG1 < 0.05); * Erre 0; * Sinon; * Erre 5; * Finsi; *OPTI DONN 5 echo 1 trac X; lcrack1 = lcrack0; ******************************************************* * BOUCLE TEMPORELLE POUR FAIRE PROPAGER LA FISSURE ******************************************************* it1 = 0 ; si(flfin); si(fllong); nt1 = 18; sino; nt1 = 3; fins; sino; si(fllong); nt1 = 7; sino; nt1 = 3; fins; fins; REPE BT1 nt1; it1 = it1 + 1; saut lign; mess '***********************' it1 '*************************'; TITR (chai 'iteration ' it1); ******************************************************* * propagation élémentaire de fissure *** Increment et angle de propagation *** * increment manuel * da1 = 1.5 * dx1; * increment auto pptip1 = poin sa1 'PROCH' ptip1; dx1tip = coor 3 pptip1; da1 = 1.5 * dx1tip; * angle si((abs K2G)