* fichier : xfem3d_01.dgibi
*******************************************************
*opti echo 0;
MESS '' ;
MESS ' CALCUL ELASTIQUE DE PROPAGATION STATIQUE DE FISSURE ';
MESS ' en 3D XFEM ' ;
MESS '' ;
*******************************************************
*** Options
*** options de calcul
opti dime 3 elem CUB8 mode trid;
opti epsilon lineaire;
*** options de trace
graph = faux ;
*graph = vrai ;
opti 'TRAC' PSC ;
*opti 'TRAC' X ;
*******************************************************
*** Maillage : cube sain 1 x 1 x 1
* finesse de la discretisation
* n1 = 36;
* n1 = 24;
* n1 = 12;
n1 = 8;
* geoemtrie
L1 = 1. ;
dx1 = (L1) / (flot n1) ;
isov1 = dx1 * (prog -6. PAS 0.5 6.);
* face avant (a)
DENS (dx1) ;
pa1 = (-0.5*L1) 0. (-0.5*L1) ;
pa2 = ( 0.5*L1) 0. (-0.5*L1) ;
pa3 = ( 0.5*L1) 0. ( 0.5*L1) ;
pa4 = (-0.5*L1) 0. ( 0.5*L1) ;
dra1 = pa1 droi n1 pa2;
* repartition 0.25 // 0.5 // 0.25
vz1 = 0. 0. (0.25*L1);
sa1 = dra1 TRAN (n1/4) vz1;
dra10 = sa1 cote 3;
sa2 = dra10 TRAN (n1/2) (2.*vz1);
dra20 = sa2 cote 3;
sa3 = dra20 TRAN (n1/4) vz1;
dra3 = sa3 cote 3;
* volume = cube
ex = 1. 0. 0.;
ey = 0. 1. 0.;
ez = 0. 0. 1.;
ny1 = n1;
vol1 = sa2 VOLU 'TRAN' ny1 (L1*ey);
vol2 = sa1 VOLU 'TRAN' ny1 (L1*ey);
vol3 = sa3 VOLU 'TRAN' ny1 (L1*ey);
dx1elim = 1.E-5*dx1;
ELIM vol1 vol2 dx1elim;
ELIM vol1 vol3 dx1elim;
vol23 = vol2 et vol3;
vol0 = vol1 et vol23;
* qq recup
pa0 = poin vol1 'PROCH' (0. 0. 0.);
pa0x = poin vol1 'PROCH' (0. L1 0.);
pa0z = (poin vol1 'PROCH' ((-0.5*L1) L1 0.))
et (poin vol1 'PROCH' ((0.5*L1) L1 0.));
sb1 = (dra3 TRAN ny1 (L1*ey)) coul 'BLEU';
ELIM sb1 vol3 (1.E-4 * dx1);
sc1 = (dra1 TRAN ny1 (L1*ey)) coul 'OCEA';
ELIM sc1 vol2 (1.E-4 * dx1);
con0 = aret vol0;
env0 = ENVE vol0;
eye1 = 100. -50. 30. ;
eye2 = 0. -0. 100. ;
si(graph);
trac eye1 (vol0 et sb1 et sc1) 'CACH';
fins;
mess (nbno vol23) (nbno vol1) (nbno vol0);
mess (nbel vol23) (nbel vol1) (nbel vol0);
**********************************************************
* tables et indice de propagation
it = 0;
pcrack1 = tabl ;
crack1 = tabl;
che1X = tabl;
***********************************************************
*** Maillage de la fissure
*** fissure semi elliptique aplatie
a0 = 0.50;
rapab0 = 0.2;
b0 = a0 * rapab0;
bbord = 0.20;
* finesse du maillage de la fissure
da0 = dx1 / 4.;
na0 = enti (a0 / da0);
nb0 = enti (b0 / da0);
&ELEM = VALE 'ELEM';
opti elem TRI3;
* da000 = rapab0 * da0;
da000 = da0;
* points de base
DENS (da000);
pf1 = a0 0. 0. ;
pf2 = a0 bbord 0. ;
pf4 = (-1.*a0) bbord 0. ;
pf5 = (-1.*a0) 0. 0. ;
DENS (da0);
pf0 = 0. 0. 0.;
pf3 = 0. bbord 0.;
* lignes de base
lf1 = (pf5 droi pf0) droi pf1;
lf2 = (pf4 droi pf3) droi pf2;
* profil elliptique
un2 = manu lf2 'CHPO' 'SCAL' 1.;
x2 = coor lf2 1;
vy2 = b0 * ((un2 - ((x2/a0)**2))**0.5);
vy2 = vy2 nomc 'UY';
lfron1 = (lf2 plus vy2) coul 'ORAN';
* surface reglee
sfron1 = (regl lfron1 lf1) coul 'ROUG';
* tres important
opti elem &ELEM ;
si(graph);
trac eye2 sfron1;
trac eye1 (con0 et sfron1 ) 'CACH';
trac eye1 (con0 et lfron1) 'CACH';
fins;
pcrack1.it = lfron1;
crack1.it = sfron1 ;
psi1 phi1 = PSIP vol1 crack1.it 'DEUX' pcrack1.it;
si(graph);
opti isov surf ;
trac phi1 vol0 isov1 'COUP' pf0 (0. 0. 1.) pf3;
trac psi1 vol0 isov1 'COUP' pf0 pf1 pf3;
fins;
*opti donn 5 ;
**********************************************************
*** MODELE & MATERIAU ***
**********************************************************
* coef elastique
Ey0 = 2.E11;
nu0 = 0.3 ;
rho0= 7800. ;
mu0 = (0.5 * Ey0) / (1. + nu0) ;
Estar0 = Ey0 / (1. - (nu0**2));
* Loi de Paris (on convertit C : mm -> m)
mParis = 3.;
CParis = 1.E-9 * 1.E-3;
**********************************************************
*** Modele & materiau ***
*
* zone de propagation (X-FEM)
mod1 = MODE vol1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'XC8R';
mat1 = MATE mod1 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
che1X.it = TRIE mod1 psi1 phi1 ;
*
* elements standards
mod23 = MODE vol23 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'CUB8';
mat23 = MATE mod23 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
*
mod1tot = mod1 et mod23;
mat1tot = mat1 et mat23 ;
**********************************************************
*** CL et DEPLACEMENTS IMPOSES ***
**********************************************************
syy0 = -100.E6;
*** MODE I ***
* CL mode I
cl1 = bloq pa0 'UX' 'UY' 'UZ';
cl2 = bloq pa0x 'UZ';
cl3 = bloq pa0z 'UX' 'UZ';
cl1tot = (cl1 ET cl2 ET cl3) ;
* chargement mode I
for1 = PRES 'MASS' mod23 sb1 syy0;
for2 = PRES 'MASS' mod23 sc1 syy0;
for1tot = for1 et for2;
si(graph);
trac (vect for1tot 'FORC' 'BLEU')
(con0 et pcrack1.it) 'CACH';
fins;
*** matrice elastique et assemblage
Rig1tot = RIGI mod1tot mat1tot;
K1tot = Rig1tot et cl1tot;
*******************************************************
* INITIALISATIONS POUR LA BOUCLE TEMPORELLE
*******************************************************
logda = FAUX;
u1 = tabl;
u1 . 0 = RESO K1tot (0.*for1tot);
TKchpo1 = tabl; TK1ev = tabl; TK2ev = tabl; TK3ev = tabl;
*******************************************************
* BOUCLE TEMPORELLE POUR FAIRE PROPAGER LA FISSURE
*******************************************************
*nt1=n1;
nt1=2;
REPE BT1 nt1;
it = it + 1;
saut lign;
mess '***********************' it '*************************';
TITR (chai 'iteration ' it);
*******************************************************
* PROPAGATION ÉLÉMENTAIRE DE FISSURE
*******************************************************
si(logda);
*** Angle de propagation ***
un1 = manu 'CHPO' pcrack1 . (it - 1) 'SCAL' 1.;
* angle
sgnK2 = (masq K2G 'EGSUPE' 0.) - (masq K2G 'INFERIEUR' 0.);
* modif pour etre tranquille
masK2 = masq (abs K2G) 'EGINFE' (1.E-4*(abs K1G));
K2G = (masK2) + ((un1 - masK2) * K2G);
tet0c = K1G / K2G;
tet8c = (tet0c**2) + (8. * un1);
tet1c = tet0c - ( sgnK2 * (tet8c**0.5) );
tet1c = 2. * (atg (0.25 * tet1c));
tet1c = (0. * masK2 ) + ((un1 - masK2) * tet1c);
* on recupere la direction du repere local de la fissure (CHPOINT)
ETIP1 = KTAB . 'UTILTET1' . 'DIRECTION1';
ETIP2 = KTAB . 'UTILTET1' . 'DIRECTION2';
* on fait tourner de l'angle critique
DA1VEC = ((cos tet1c) * ETIP1) + ((sin tet1c) * ETIP2);
*trac (vect DA1VEC) (pcrack1 . it et con0);
*** vitesse de propagation ***
K1eff = ( (K1G**2) + (K2G**2) )**0.5;
K1eff = ( (K1G**2) + (K2G**2) + ((Estar0/mu0)*(K3G**2)) )**0.5;
* Passage en MPa.m^0.5
K1eff = K1eff * 1.E-6;
* quel est le dn correspondant a da1min ?
K1min= mini K1eff 'ABS';
dadn = CParis * (K1min**mParis);
dn = da1min / dadn;
dadn = (CParis * dn) * (K1eff**mParis) ;
DA1VEC = dadn * DA1VEC;
*** on cree un maillage par translation ****
trac (vect DA1VEC 1) (pcrack1 . (it - 1) et con0);
pcrk1 = pcrack1 . (it - 1) PLUS DA1VEC;
*** quelques chpo utiles ***
* on calcul la distance signee a la paroi
dis0 = PSIP pcrack1 . (it - 1) env0;
dis1 = PSIP pcrk1 env0;
* on cree un chpoint permettant de savoir si un point est a l interieur
unvol1 = manu 'CHML' vol0 'SCAL' 1.;
uncrk1 = PROI pcrk1 unvol1;
*** Etape 1 : traitement des segments "a cheval" sur le bord ***
* boucle sur les segments de pcrk1
dismall1 = 1.E-6 * da0;
ncrk1 = nbel pcrk1;
icrk1 = 0;
repe Bcrk1 ncrk1;
icrk1 = icrk1 + 1;
* recup des segments et des points
pcrk1i = ELEM pcrk1 icrk1;
pt1a = poin pcrk1i 'INITIAL'; pt1b = poin pcrk1i 'FINAL';
* test pour savoir si a cheval
un1a = extr uncrk1 pt1a 'SCAL'; un1b = extr uncrk1 pt1b 'SCAL';
si((abs(un1a - un1b)) >eg 0.5);
phi1a = EXTR dis1 pt1a 'PHI'; phi1b = EXTR dis1 pt1b 'PHI';
si(un1a < 0.5);
* il faut ramener le point pt1a
v1a = phi1a / (phi1a - phi1b);
v1a = v1a * (pt1b moin pt1a);
DEPL pt1a 'PLUS' v1a;
sino;
* il faut ramener le point pt1b
v1b = phi1b / (phi1b - phi1a);
v1b = v1b * (pt1a moin pt1b);
DEPL pt1b 'PLUS' v1b;
fins;
fins;
fin Bcrk1;
*** Etape 2 : traitement des points encore a l'exterieur ***
* mise a jour de uncrk1
uncrk1 = PROI pcrk1 unvol1;
icrk1 = 0;
repe Bcrk2 ncrk1;
icrk1 = icrk1 + 1;
* recup des segments et des points
pcrk1i = ELEM pcrk1 icrk1;
pt1a = poin pcrk1i 'INITIAL';
* pt1a a l exterieur ?
un1a = EXTR uncrk1 pt1a 'SCAL';
si(un1a < 0.5);
pcrk0i = ELEM pcrack1 . (it - 1) icrk1;
pt0a = poin pcrk0i 'INITIAL';
* distance signée associée
phi0a = EXTR dis0 pt0a 'PHI';
phi1a = EXTR dis1 pt1a 'PHI';
* il faut ramener le point pt1a
v1a = phi1a / (phi0a - phi1a);
v1a = v1a * (pt1a moin pt0a);
DEPL pt1a 'PLUS' v1a;
fins;
fin Bcrk2;
* sans oublier le dernier point du dernier segment
pt1a = poin pcrk1i 'FINAL';
* pt1a a l exterieur ?
un1a = EXTR uncrk1 pt1a 'SCAL';
si(un1a < 0.5);
* recup des segments et des points du front precedent
pcrk0i = ELEM pcrack1 . (it - 1) ncrk1;
pt0a = poin pcrk0i 'FINAL';
* distance signée associée
phi0a = EXTR dis0 pt0a 'PHI';
phi1a = EXTR dis1 pt1a 'PHI';
* il faut ramener le point pt1a
v1a = phi1a / (phi0a - phi1a);
v1a = v1a * (pt1a moin pt0a);
DEPL pt1a 'PLUS' v1a;
fins;
*** si tous les points du front sont sortis
* c est qu on a tout cassé !!! => on quitte
meshun1 = extr uncrk1 'MAILLAGE';
si(ega (TYPE meshun1) 'ENTIER'); quit BT1; fins;
*** Etape 3 : points sur le bord qui sont entrés a linterieur ***
* mise a jour de dis1
dis1 = PSIP pcrk1 env0;
* 1er point du 1er segment
pcrk0i = ELEM pcrack1 . (it - 1) 1;
pt0a = poin pcrk0i 'INITIAL';
phi0a = EXTR dis0 pt0a 'PHI';
si((abs phi0a) dismall1);
pt0b = poin pcrk0i 'FINAL';
phi0b= EXTR dis0 pt0a 'PHI';
pt1b = poin pcrk1i 'FINAL';
phi1b= EXTR dis1 pt1b 'PHI';
v1a = phi1a / (phi1b - phi1a) ;
v1a = v1a * (pt1a moin pt1b);
DEPL pt1a 'PLUS' v1a;
fins;
fins;
* dernier point du dernier segment
pcrk0i = ELEM pcrack1 . (it - 1) ncrk1;
pt0b = poin pcrk0i 'FINAL';
phi0b= EXTR dis0 pt0a 'PHI';
si((abs phi0b) dismall1);
pt0a = poin pcrk0i 'INITIAL';
phi0a = EXTR dis0 pt0a 'PHI';
pt1a = poin pcrk1i 'INITIAL';
phi1a = EXTR dis1 pt1a 'PHI';
v1b = phi1b / (phi1a - phi1b) ;
v1b = v1b * (pt1b moin pt1a);
DEPL pt1b 'PLUS' v1b;
fins;
fins;
*** on recupere le vrai nouveau front ***
pcrack1 . it = ELEM pcrk1 'APPUYE' 'STRICT' meshun1;
si((flot (it / 2)) ega ((flot it) / 2.));
pcrack1 . it = pcrack1 . it coul 'ROUG';
sinon;
pcrack1 . it = pcrack1 . it coul 'ORAN';
fins;
*** on maille la surface fissuree ***
&ELEM = VALE 'ELEM';
opti 'ELEM' TRI3;
nregl1 = 1 + (enti (da1min / da0));
sfiss1 = regl pcrk1 nregl1 pcrack1 . (it - 1) ;
*** on retraite la surface fissuree cree au debut ***
ELIM sfiss1 dismall1;
sfiss1 = (REGE sfiss1) elem 'TRI3';
crack1 . it = sfiss1 et crack1 . (it - 1) ;
opti 'ELEM' &ELEM;
*** fin de l avancee du maillage du front **********************
*** Actualisations ***
* ...des level set
psi1 phi1 = PSIP vol1 crack1.it 'DEUX' pcrack1.it;
si(graph);
trac phi1 vol0 isov1 'COUP' pf0 (0. 0. 1.) pf3;
trac psi1 vol0 isov1 'COUP' pf0 pf1 pf3;
fins;
* ...du modele, de la rigidite,... (inutile pour le materiau)
Che1X . it = TRIE mod1tot psi1 phi1;
Rig1tot = RIGI mod1tot mat1tot;
K1tot = Rig1tot et cl1tot;
*******************************************************
* PAS DE PROPA
*******************************************************
sino;
crack1 . it = crack1 . (it - 1);
pcrack1 . it = pcrack1 . (it - 1);
fins;
trac (con0 et crack1 . it) eye2;
****************************
*** resolution ***
u1 . it = RESO K1tot for1tot;
u1phy = XFEM 'RECO' u1 . it mod1tot ;
sig1 = SIGM mod1tot u1 . it mat1tot;
ucrk1u ucrk1d = XFEM 'FISS' crack1 . it u1 . it mod1tot ;
si(graph);
* amp1 = 25.;
amp1 = 10.;
def1 = (DEFO u1phy vol0 amp1) ;
TRAC sig1 mod1tot def1 eye1;
def2 = (DEFO u1phy con0 amp1)
et (DEFO ucrk1u crack1 . it amp1)
et (DEFO ucrk1d crack1 . it amp1);
TRAC def2 eye1;
def3 = def1
et (DEFO ucrk1u crack1 . it amp1)
et (DEFO ucrk1d crack1 . it amp1);
TRAC def3 eye1 'COUP' pf0 (0. 0. 1.) pf3 'CACH';
finsi;
*** Calcul des FIC par methode integrale G_THETA ***
KTAB = tabl ;
KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
KTAB .'PSI' = psi1;
KTAB .'PHI' = phi1;
KTAB .'FRONT_FISSURE' = pcrack1 . it ;
KTAB .'MODELE' = mod1tot;
KTAB .'CARACTERISTIQUES' = mat1tot;
KTAB .'SOLUTION_RESO' = u1 . it;
KTAB .'CHARGEMENTS_MECANIQUES' = for1tot;
KTAB . 'COUCHE' = 1;
* KTAB . 'COUCHE' = 2;
G_THETA KTAB;
Kchpo1 = KTAB . 'CHPO_RESULTATS' ;
TKchpo1 . it = Kchpo1;
K1ev = EVOL BLEU 'CHPO' Kchpo1 'K1' pcrack1 . it;
K2ev = EVOL DEFA 'CHPO' Kchpo1 'K2' pcrack1 . it;
K3ev = EVOL ROUG 'CHPO' Kchpo1 'K3' pcrack1 . it;
TK1ev . it = K1ev;
TK2ev . it = K2ev;
TK3ev . it = K3ev;
S1pr = EXTR K1ev 'ABSC';
K1pr = EXTR K1ev 'ORDO';
K2pr = EXTR K2ev 'ORDO';
K3pr = EXTR K3ev 'ORDO';
* K1eqpr = ((K1pr**2) + (K2pr**2))**0.5;
K1eqpr = ( (K1pr**2) + (K2pr**2) + ((Estar0/mu0)*(K3pr**2)) )**0.5;
K1eqev = EVOL 'VERT' 'MANU' 's' S1pr 'K1eq' K1eqpr;
si(it ega 1);
tdess1 = tabl;
tdess1 . 1 = 'MOT' 'MARQ PLUS TIRR ';
tdess1 . 2 = 'MOT' 'MARQ CROI TIRR ';
tdess1 . 3 = 'MOT' 'TIRM ';
fins;
dess (K1ev et K2ev et K3ev et K1eqev) tdess1 'TITR' 'SIF';
dadnev = CParis * (K1eqev**mParis);
dess dadnev 'TITR' 'DADN';
K1G = EXCO Kchpo1 'K1';
K2G = EXCO Kchpo1 'K2';
K3G = EXCO Kchpo1 'K3';
GTAB = tabl ;
GTAB .'OBJECTIF' = MOT 'J ';
GTAB .'PSI' = psi1;
GTAB .'PHI' = phi1;
GTAB .'FRONT_FISSURE' = pcrack1 . it ;
GTAB .'MODELE' = mod1tot;
GTAB .'CARACTERISTIQUES' = mat1tot;
GTAB .'SOLUTION_RESO' = u1 . it;
GTAB .'CHARGEMENTS_MECANIQUES' = for1tot;
GTAB .'COUCHE' = 1;
G_THETA GTAB;
Jchpo1 = GTAB . CHPO_RESULTATS;
Jev1 = EVOL ROSE 'CHPO' Jchpo1 'J' pcrack1 . it;
dess Jev1;
*** Increment de propagation ***
logda = VRAI;
* increment manuel
da1min = 1.4 * dx1;
* OPTI DONN 5 echo 1 trac X;
FIN BT1;
*******************************************************
*******************************************************
*** erreur? (avant propa elementaire) ***
*******************************************************
K1ev1 = TK1ev . 1;
K2ev1 = TK2ev . 1;
K3ev1 = TK3ev . 1;
K1ev1m = maxi (extr K1ev1 'ORDO');
K2ev1m = maxi (extr K2ev1 'ORDO');
K3ev1m = maxi (extr K3ev1 'ORDO');
* mess K1ev1m K2ev1m K3ev1m;
K1ref1m = 1.59717E+08;
K2ref1m = 0.;
K3ref1m = 0.;
errK1 = abs ((K1ev1m - K1ref1m) / K1ref1m);
errK2 = abs ((K2ev1m - K2ref1m) / K1ref1m) ;
errK3 = abs ((K2ev1m - K2ref1m) / K1ref1m);
errK1 = errK1 + errK2 + errK3;
*******************************************************
*** erreur? (apres propa elementaire) ***
*******************************************************
Jmaxref1 = 8.14855E+05;
Jminref1 = 4.75952E+05;
Jmin1 = mini Jchpo1;
Jmax1 = maxi Jchpo1;
errJ1 = (abs ((Jmin1 - Jminref1) /Jminref1))
+ (abs ((Jmax1 - Jmaxref1) /Jmaxref1));
mess 'erreur relative sur K1+K2+K3 et sur J';
mess ' ' errK1 errJ1;
*******************************************************
*2% d ecart acceptee
*******************************************************
Si ((errJ1 < 2.E-2) ou (errK1 < 2.E-2));
Erre 0;
Sinon;
Erre 5;
Finsi;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer