******************************************************************* ** raff03.dgibi ******************************************************************* * * Calcul en mecanique de la rupture avec un maillage raffine par * l'operateur RAFF et des elements X-FEM. * * d'une plaque elastique en traction en 2D déformations planes * avec fissure horizontale * test avec des éléments QUA4 QUA6 TRI3 et TRI6 * * Dimentions du quart de plaque: * l1= 0.5 m * l2= 0.5 m * Longueur de la fissure : * a=0.005 m * Inclinaison : * alpha=20 degre * *création : gg, le 15.03.2017 * * Comparaison de la simulations avec une solution de reference * pour g,k1 et k2 donne par : * [TADA, STRESS ANALYSIS HANDBOOK,1973] * [ISIDA,IJF(7),301-316,1971] * ************************************************************************ *** Options de calcul : 2D opti dime 2 elem qua4 mode plan defo; opti trac x; *opti trac psc; *graph = vrai; graph = faux; ******************************************************* *** Maillage : Plaque rectangulaire saine * dimension et finesse Lb1 = 0.50 ; nb1 = 20; dx2 = (2.*Lb1) / (flot nb1) ; dx1 = dx2/16; P1 = (-1*Lb1) (-1*Lb1) ; P2 = (-1*Lb1) Lb1 ; P3 = Lb1 Lb1 ; P4 = Lb1 (-1*Lb1) ; L1 = DROITE nb1 P1 P2 ; L2 = DROITE nb1 P2 P3 ; L3 = DROITE nb1 P3 P4 ; L4 = DROITE nb1 P4 P1 ; stot1 = DALL L1 L2 L3 L4 ; ******************************************************* *** 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; ******************************************************* * raffinement du maillage proche des fronts * rayon min de la zone raffiné dmin= 1.5*dx2; dmax=dmin+(dx2); x1 y1 = coor stot1 ; DF1= (((x1 - (a0 * (cos alpha0)))**2)+ ((y1 - (a0 * (sin alpha0)))**2))**0.5; DF2= (((x1 + (a0 * (cos alpha0)))**2)+ ((y1 + (a0 * (sin alpha0)))**2))**0.5; *DF3= (0.5*DF1)+(0.5*DF2); DF3= mini DF1 DF2; aff = (((dx2 - dx1)*(DF3 - dmax))/(dmax - dmin)) + (dx2) ; dens1 = born aff 'SCAL' 'COMPRIS' (dx1*1.000001) (dx2*1.000001) ; si (graph); TRAC dens1 stot1 TITR 'Fonction de densite pour RAFF'; finsi; stot2 = raff stot1 dens1 ; PTX = STOT2 POIN 'DROI' ptip1 ptip2 (dx2*0.5); PTX = PTX POIN 'DROI' (0.0 0.0) ((-1.0*(sin alpha0)) (cos alpha0)) (a0+(dx2*0.5)); SURFX = STOT2 ELEM 'APPUYE' 'LARGEMENT' PTX ; SURFC = STOT2 DIFF SURFX ; si (graph); TRAC 'FACE' ((SURFC coul jaune) ET (SURFX COUL BLEU) ET lcrack0 ) TITR 'Surface FEM Jaune et X-FEM BLEU'; finsi ; * creation des level set psi0 phi0 = PSIPHI SURFX lcrack0 'DEUX' ptip1 ptip2; isolv7 = prog (-2.*a0) PAS (a0 / 4.) (2.*a0); si (graph); TRAC (SURFX et lcrack0) TITR 'ELEMENTS X-FEM'; TRAC phi0 (SURFX et SURFC ) isolv7 ; TRAC psi0 (SURFX et SURFC ) isolv7; finsi ; *opti donn 5 ; ******************************************************* *** Modeles et materiaux * definition des models elementaires mod1x = MODE surfx mecanique elastique xq4R; mod1f = MODE surfc mecanique elastique; *mod1tot = MODE stot2 mecanique elastique isotrope xq4R; * enrichissement Che1X = tabl; Che1X . 0 = TRIE mod1x psi0 phi0; * constructionsion des blocages des ddl X-fem non actifs dans * les éléments de transition et des relation de conformitées dues au * raffinent. Rel1X = tabl; mod1tot = mod1x et mod1f; Rel1X.0 = RELA mod1tot; list resu Rel1X.0; *list resu Che1X . 0; * apres cela on assemble mod1tot = mod1x et mod1f; * 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; list resu RIG1TOT; * on bloque les depl de corps rigides pb1 = surfc POIN 'PROC' (0.0 Lb1); pb2 = surfc POIN 'PROC' (0.0 (-1.0*Lb1)); pb3 = surfc POIN 'PROC' (Lb1 0.0); pb4 = surfc POIN 'PROC' ((-1.0*Lb1) 0.0); cl1x = BLOQ 'UX' (pb1 et pb2); cl1y = BLOQ 'UY' (pb3 et pb4); cl1tot = cl1x et cl1y; K1tot = Rig1tot et cl1tot et Rel1X.0 ; syy0 = 1.E9; pre1 = PRES 'MASSIF' mod1tot (L4 et L2) (-1.*syy0); si (graph); TRAC (VECT pre1 'FORC' 'BLEU') (stot2 et lcrack0) TITR 'Chargement'; finsi; list resu cl1tot; list resu Rel1X.0; ******************************************************* *** RESOlution avec RESO u0 = RESO K1tot pre1; ******************************************************* *** Post traitement * verif du fonctionnement des principaux operateurs defel0 = EPSI mod1tot u0 'LINE'; sig0 = SIGM mod1tot u0 mat1tot 'LINE'; 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 ; si (graph); def0 = DEFO u0phy stot2 20.; TRAC sig0 mod1tot def0 'TITR' 'Contraintes'; ucrk0u ucrk0d = XFEM 'FISS' lcrack0 u0 mod1x ; def0 = def0 et (DEFO ucrk0u lcrack0 20. coul rouge) et (DEFO ucrk0d lcrack0 20. coul vert); TRAC def0 TITR 'DEFORMEE'; * trac defel0 mod1tot (lcrack0); * TRAC svm0 mod1tot (stot1 et lcrack0); * TRAC wdefel0 mod1tot (stot1 et lcrack0); * TRAC gru0 mod1tot (stot1 et lcrack0); finsi; ******************************************************* * 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 calcule J pour differents contours prj1 = prog ; prcou1= prog ; lecou1 = lect 1 pas 1 5; 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 TITR 'J en fonction du contour'; fins; * extraction d un J moyen jmoy1 = jmoy1 / ncou1; *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; 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)) ; ********************************************************* 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 ); Si (errG1 < 0.01); 'MESS' ' ----------------------' ; 'MESS' ' SUCCES DU CAS-TEST !' ; 'MESS' ' ----------------------' ; Erre 0; Sinon; 'MESS' ' ---------------------' ; 'MESS' ' ECHEC DU CAS-TEST !' ; 'MESS' ' ---------------------' ; Erre 5; finsi; fin;