* fichier : xfem01.dgibi
opti echo 1 ;
*******************************************************************
*
* 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 0 ;
*** 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;
mod1b = MODE sb1 mecanique elastique isotrope;
* apres cela on assemble
mod1tot = mod1a et mod1b;
* enrichissement
Che1X = tabl;
Che1X . 0 = TRIE mod1tot psi0 phi0;
* 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;
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);
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)
© Cast3M 2003 - All rights reserved.
Disclaimer