* fichier : xfem_gd.dgibi
opti DIME 2 ELEM QUA4 mode plan defo;
********************************************************************
* xfemcapipica.dgibi
*
* Test de l'opérateur de passage des contraintes(déformations)
* PK2 aux contraintes de Cauchy pour la XFEM
* d'une plaque elastique en traction avec fissure droite
*
* création : as, le 14.01.2010
********************************************************************
opti echo 0;
* ajout de option epsilon lineaire pour la precision des test!
OPTION epsilon lineaire;
* Options :
*########
* pour afficher les messages et les graphes + choix de la sortie en .ps:
lmess = faux ; lgraph = faux ;
si lgraph ; lps=vrai; sinon; lps= faux; finsi;
si lps; opti trac psc; opti ftra 'xfem_gd.ps'; finsi;
* Maillage :
* ########
dx1 = 0.06 ;
blim = 1. + 5 * dx1 ; nbr1 = enti (blim / dx1);
* densité sur la hauteur de la bande :
nbh1 = 3 ; hh1 = ((flot nbh1)*dx1);
p1 = 'POIN' 0. hh1 ; p2 = 'POIN' blim hh1 ;
pm3 = 'POIN' blim (0.) ; pm4 = 'POIN' 0. (0.) ;
p3 = 'POIN' blim (-1.*hh1) ; p4 = 'POIN' 0. (-1.*hh1) ;
lh = droi nbr1 p1 p2 ;
ldh = droi p2 pm3 DINI dx1 DFIN dx1;
lgh = droi p1 pm4 DINI dx1 DFIN dx1;
lm = droi nbr1 pm4 pm3;
ldb = droi pm3 p3 DINI dx1 DFIN dx1;
lgb = droi pm4 p4 DINI dx1 DFIN dx1;
lb = droi nbr1 p4 p3;
s1h = 'DALLER' lm ldh lh lgh ;
s1b = 'DALLER' lb ldb lm lgb ;
* fissure initiale :
* #################
a0 = 3.*dx1;
pfis0 = -0.01 0.; pfisi = a0 0.;
fis0 = coul (droi 1 pfis0 pfisi) roug;
* Modèle et materiau :
* ####################
Ey1 = 2.e5 ; nu1=0.3; rho1 = 7.e-6 ;
* CHARGEMENT : déplacement imposé :
* #################################
tf = 1.; NT = ENTI (3.*tf);
dt0 = tf / (flot NT) ; l_t = PROG 0. pas dt0 tf;
ufin = tf ; l_uimp = PROG 0. pas (ufin / (flot NT)) ufin;
evou = evol 'MANU' 'temps' l_t 'U_imp' l_uimp ;
coeff0 = 0.2;
**###################
** Calcul MASSIF **
**###################
s1q = s1h ;
elim s1q 0.0001;
* Modèle et materiau :
* ####################
Mod1q = s1q MODE mecanique elastique isotrope ;
mat1q = MATE Mod1q 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
* Conditions aux limites, Déplacements imposés :
* ################################################
ptrach = poin s1q proc (p1 plus ((2.*dx1) (-2.*dx1) ));
cl1hx = 'BLOQUE' ptrach 'UX' ; cl1hy = 'BLOQUE' ptrach 'UY' ;
psyms = elem lm comp (poin s1q proc pfisi) pm3;
cl1sym = 'BLOQUE' psyms 'UY' ;
si lgraph ; trac (s1q et (coul psyms roug)); finsi;
cl0 = (cl1hx 'ET' cl1hy 'ET' cl1sym) ;
dep1 = DEPI cl1hy coeff0;
uimp0 = CHAR 'DIMP' dep1 evou ;
* Resolution iterative du probleme
* ################################################
TAB0 = TABL; TAB0.TEMPS0 = 0. ;
TAB0.'MODELE' = Mod1q ;
TAB0.'CARACTERISTIQUES' = mat1q ;
TAB0.'CHARGEMENT' = uimp0;
TAB0.'BLOCAGES_MECANIQUES' = cl0 ;
TAB0.'TEMPS_CALCULES' = L_t ; TAB0.'TEMPS_SAUVES' = L_t ;
TAB0.'CONVERGENCE_FORCEE' = faux;
TAB0.'DELTAITER' = 190;
TAB0.'MAXITERATION' = 199;
PASAPAS tab0;
*opti donn 5;
**##################
** Calcul XFEM **
**##################
s1 = s1h et s1b;
elim s1 0.0001;
si lgraph; trac (s1 et fis0); finsi;
* Modèle et materiau :
* ####################
Mod1a = s1 MODE mecanique elastique isotrope xq4R;
Mod1 = Mod1a ;
mat1 = MATE Mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
** Table XFEM : fisssure initiale, critere d'initiation, modele zone cohesive
** #################################################################
tab_xfem = tabl;
TAB_xfem.fissure = fis0 ;
* creation des level set
psi1 phi1 = PSIPHI s1 fis0 'DEUX' pfisi;
* Appel a TRIELE :
cheX1 = trie mod1 psi1 phi1 'SAUT';
* Conditions aux limites, Déplacements imposés :
* ################################################
ptrach = poin s1 proc (p1 plus ((2.*dx1) (-2.*dx1) ));
ptracb = poin s1 proc (p4 plus ((2.*dx1) (2.*dx1) ));
cl1hx = 'BLOQUE' ptrach 'UX' ; cl1bx = 'BLOQUE' ptracb 'UX' ;
cl1hy = 'BLOQUE' ptrach 'UY' ; cl1by = 'BLOQUE' ptracb 'UY' ;
cl1 = (cl1hx 'ET' cl1bx 'ET'cl1hy 'ET' cl1by) ;
* CHARGEMENT : déplacement imposé :
dep1 = (DEPI cl1hy coeff0) et (DEPI cl1by (-1.*coeff0)) ;
uimp0 = CHAR 'DIMP' dep1 evou ;
* Resolution iterative du probleme
* ################################################
TAB1 = TABL; TAB1.TEMPS0 = 0. ;
TAB1.'MODELE' = Mod1 ;
TAB1.'CARACTERISTIQUES' = mat1 ;
TAB1.'CHARGEMENT' = uimp0;
TAB1.'BLOCAGES_MECANIQUES' = cl1 ;
TAB1.'TEMPS_CALCULES' = L_t ; TAB1.'TEMPS_SAUVES' = L_t ;
TAB1.'CONVERGENCE_FORCEE' = faux;
TAB1.'DELTAITER' = 190;
TAB1.'MAXITERATION' = 199;
*opti donn 5 ;
pasapas tab1;
**######################
** Post-traitement **
*########################
naff = mini (lect ((dime tab1.deplacements) - 1)
((dime tab0.deplacements) - 1));
si lgraph;
trac ((defo s1q tab1.deplacements.naff 1.)
et (defo (coul s1h roug) tab0.deplacements.naff 1.));
finsi;
conf0=form;
**#####################
* TEST de PICA et CAPI
**######################
CAU1 = pica tab1.contraintes.naff tab1.deplacements.naff
tab1.deplacements. 0 tab1.modele;
PIC1 = capi CAU1 tab1.deplacements.naff tab1.deplacements. 0
tab1.modele;
CAU0 = pica tab0.contraintes.naff tab0.deplacements.naff
tab0.modele;
PIC0 = capi CAU0 tab0.deplacements.naff tab0.modele;
si lmess ;
mess ' ';
mess 'Test de CAPI et PICA en grandes defs ';
mess ' ';
mess 'Configuration de départ = configuration initiale';
mess '----------------------';
mess ' ';
mess 'Elements enrichis :';
finsi;
el1 = 3; pg1 = 61;
sig01 = extr tab1.contraintes.naff 'SMYY' 1 el1 pg1;
cau11 = extr CAU1 'SMYY' 1 el1 pg1;
pic11 = extr PIC1 'SMYY' 1 el1 pg1;
si lmess; mess 'XFEM:' 'ini' sig01 'pica1' cau11 'inicapi' pic11; fins;
pg0 = 3;
sig00 = extr tab0.contraintes.naff 'SMYY' 1 el1 pg0;
cau01 = extr CAU0 'SMYY' 1 el1 pg0;
pic01 = extr PIC0 'SMYY' 1 el1 pg0;
si lmess; mess 'QUA4:' 'ini' sig00 'pica1' cau01 'inicapi' pic01; fins;
err6001 = abs ((sig01-sig00)/sig00);
err6002 = abs ((cau11-cau01)/cau01);
err6003 = abs ((pic11-pic01)/pic01);
si lmess; mess ' '; fins;
el3 = 2;
sig013 = extr tab1.contraintes.naff 'SMYY' 1 el3 pg1;
cau13 = extr CAU1 'SMYY' 1 el3 pg1;
pic13 = extr PIC1 'SMYY' 1 el3 pg1;
si lmess; mess 'XFEM:' 'ini' sig013 'pica1' cau13 'inicapi' pic13;fins;
sig003 = extr tab0.contraintes.naff 'SMYY' 1 el3 pg0;
cau03 = extr CAU0 'SMYY' 1 el3 pg0;
pic03 = extr PIC0 'SMYY' 1 el3 pg0;
si lmess; mess 'QUA4:' 'ini' sig003 'pica1' cau03 'inicapi' pic03;fins;
err5001 = abs ((sig013-sig003)/sig003);
err5002 = abs ((cau13-cau03)/cau03);
err5003 = abs ((pic13-pic03)/pic03);
si lmess; mess ' '; fins;
si lmess; mess 'Element non enrichi:'; fins;
el2 = 11;
sig012 = extr tab1.contraintes.naff 'SMYY' 1 el2 pg1;
cau12 = extr CAU1 'SMYY' 1 el2 pg1;
pic12 = extr PIC1 'SMYY' 1 el2 pg1;
si lmess; mess 'XFEM:' 'ini' sig012 'pica1' cau12 'inicapi' pic12;fins;
sig002 = extr tab0.contraintes.naff 'SMYY' 1 el2 pg0;
cau02 = extr CAU0 'SMYY' 1 el2 pg0;
pic02 = extr PIC0 'SMYY' 1 el2 pg0;
si lmess; mess 'QUA4:' 'ini' sig002 'pica1' cau02 'inicapi' pic02;fins;
err26001 = abs ((sig012-sig002)/sig002);
err26002 = abs ((cau12-cau02)/cau02);
err26003 = abs ((pic12-pic02)/pic02);
errmax1 = maxi (prog err6001 err6002 err6003 err5001 err5002
err5003 err26001 err26002 err26003 );
* Test de CAPI et PICA en grandes defs
*
* Configuration de départ = configuration initiale';
* -------------------------------------------------
* XFEM :
uvr1 = XFEM 'RECO' tab1.deplacements. naff mod1 ;
FORM uvr1;
PIC111 = capi CAU1
(tab1.deplacements. 0 - tab1.deplacements.naff)
tab1.deplacements. naff tab1.modele;
CAU111 = pica pic111 (tab1.deplacements. 0 - tab1.deplacements.naff)
tab1.deplacements. naff tab1.modele;
* Standard :
form conf0;
form tab0.deplacements.naff;
PIC011 = capi CAU0 (tab0.deplacements. 0 - tab0.deplacements.naff)
tab0.modele;
CAU011 = pica pic011 (tab0.deplacements. 0 - tab0.deplacements.naff)
tab0.modele;
si lmess;
mess ' ';
mess 'Configuration de départ = configuration déformée';
mess '----------------------';
mess ' ';
mess 'Elements enrichis :';
fins;
el1 = 3; pg1 = 61;
sig01 = extr CAU1 'SMYY' 1 el1 pg1;
cau11 = extr CAU111 'SMYY' 1 el1 pg1;
pic11 = extr PIC111 'SMYY' 1 el1 pg1;
si lmess; mess 'XFEM:' 'ini' sig01 'pica1' cau11 'inicapi' pic11;fins;
pg0 = 3;
sig00 = extr CAU0 'SMYY' 1 el1 pg0;
cau01 = extr CAU011 'SMYY' 1 el1 pg0;
pic01 = extr PIC011 'SMYY' 1 el1 pg0;
si lmess; mess 'QUA4:' 'ini' sig00 'pica1' cau01 'inicapi' pic01;fins;
err6001 = abs ((sig01-sig00)/sig00);
err6002 = abs ((cau11-cau01)/cau01);
err6003 = abs ((pic11-pic01)/pic01);
si lmess; mess ' ';fins;
el3 = 2;
sig013 = extr CAU1 'SMYY' 1 el3 pg1;
cau13 = extr CAU111 'SMYY' 1 el3 pg1;
pic13 = extr PIC111 'SMYY' 1 el3 pg1;
si lmess; mess 'XFEM:' 'ini' sig013 'pica1' cau13 'inicapi' pic13;fins;
sig003 = extr CAU0 'SMYY' 1 el3 pg0;
cau03 = extr CAU011 'SMYY' 1 el3 pg0;
pic03 = extr PIC011 'SMYY' 1 el3 pg0;
si lmess; mess 'QUA4:' 'ini' sig003 'pica1' cau03 'inicapi' pic03;fins;
err5001 = abs ((sig013-sig003)/sig003);
err5002 = abs ((cau13-cau03)/cau03);
err5003 = abs ((pic13-pic03)/pic03);
si lmess; mess ' '; fins;
si lmess; mess 'Element non enrichi:'; fins;
el2 = 11;
sig012 = extr CAU1 'SMYY' 1 el2 pg1;
cau12 = extr CAU111 'SMYY' 1 el2 pg1;
pic12 = extr PIC111 'SMYY' 1 el2 pg1;
si lmess; mess 'XFEM:' 'ini' sig012 'pica1' cau12 'inicapi' pic12;fins;
sig002 = extr CAU0 'SMYY' 1 el2 pg0;
cau02 = extr CAU011 'SMYY' 1 el2 pg0;
pic02 = extr PIC011 'SMYY' 1 el2 pg0;
si lmess; mess 'QUA4:' 'ini' sig002 'pica1' cau02 'inicapi' pic02;fins;
err26001 = abs ((sig012-sig002)/sig002);
err26002 = abs ((cau12-cau02)/cau02);
err26003 = abs ((pic12-pic02)/pic02);
errmax2 = maxi (prog err6001 err6002 err6003 err5001 err5002
err5003 err26001 err26002 err26003 );
errmax = maxi (prog errmax1 errmax2);
si lmess;
mess 'Erreur entre calcul standard et XFEM:' errmax;
fins;
* Erreur initiale entre calcul standard et XFEM: 0.10829
* Attention, on fait des erreurs assez importantes car le nombre de
* pts de Gauss n'est pas le même dans les deux cas :
* -> différences sur le calcul des contraintes initiales
* -> le pt de gauss xq4r est le plus proche possible du pt QUA4
* mais reste different, 'où une erreur importante pour les éléments
* en pointe de fissure.
* pr etre sur, comparer avec éléments xq4r à 4pts de Gauss ->erreur~e-15
SI (errmax >EG 20.);
ERRE 5;
FINSI;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer