* fichier : xfem03.dgibi
********************************************************************
*
* xfem03.dgibi
*
* calcul avec elements xfem (XQ4R)
* d'une plaque en traction modele de Rousselier avec fissure droite
*
* création : as, le 24.11.2009 : Résultats moyens à cause des fonctions
* de pointes de fissure adaptées à la rupture fragile sans prise en compte
* de la plasticité étendue !!!
*
********************************************************************
*opti echo 1 ;
graph = faux;
logmess = faux;
* si on veut tester la décharge:
logdech = faux;
*******************************************************
*** Options de calcul : 2D
opti dime 2 elem qua4 mode plan defo;
*******************************************************
*** Maillage : Plaque rectangulaire saine
* dimensions
* bord droit, gauche, hauteur
b= 25. ; bg = -0. * b ; hb1= 0.5 * b ;
hgoup= 0.35 * b;
* finesse, pour la comparaison FEM/XFEM, un nombre d'éléments
* pair est nécessaire
nx1 = 6; dx1 = (b-bg) / (flot nx1) ;
ha1 = (flot nx1) * dx1;
* demi-zone centrale
DENS (dx1) ;
p1 = bg 0. ; p2 = b 0. ;
p3 = b hb1 ; p4 = bg hb1 ;
p11 = (bg + ((b-bg)/2.5)) 0. ; p33 = (b - ((b-bg)/2.5)) (hb1/5.);
p44 = p11 plus (0. (coor p33 2)); p22 = p33 moin (0. (coor p33 2));
pa1 = p1 plus (((b-bg)/2.) 0.);
d1g = droi (nx1/2) p1 p11; d11g = droi (nx1/2) p11 pa1;
d11d = droi (nx1/2) pa1 p22; d1d = droi (nx1/2) p22 p2;
d11 = d11d et d11g; d22 = droi (nx1/2) p22 p33;
d33 = droi nx1 p33 p44; d44 = droi (nx1/2) p44 p11;
d2 = droi (nx1/2) p2 p3; d4 = droi (nx1/2) p4 p1;
d3 = droi nx1 p3 p4;
d333 = droi (nx1/2) p33 p3;
saco = dall d11 d22 d33 d44;
sah = regl d33 (nx1/2) d3; sag = regl d44 (nx1/2) d4;
sad = dall d22 d1d d2 d333;
satot = saco et sah et sad et sag;
stot0 = satot et (syme droi p1 p2 satot);
ELIM stot0 (1.E-5*dx1);
si(logmess);
mess 'maillage central ext total' ;
mess 'nbno= ' (nbno satot) (nbno stot0);
mess 'nbel= ' (nbel satot) (nbel stot0);
finsi;
*******************************************************
*** Maillage de la fissure
* longueur de la fissure
a0 = (b - bg) /2. ;
* finesse du maillage de la fissure
na0 = 10;
da0 = a0 / (flot na0) ;
* maillage
DENS (da0) ;
ptip0 = (bg) (0.);
ptip1 = (a0 + bg) (0.);
lcrack0 = (ptip0 droi na0 ptip1) coul 'ROUG';
* table (+ pratique pour gérer la propa)
TXFEM = tabl;
TXFEM . 0 = tabl;
TXFEM . 0 . 'POINTE1' = ptip1;
TXFEM . 0 . 'FISSURE' = lcrack0;
* creation des level set
psi0 phi0 = PSIPHI stot0 lcrack0 'DEUX' ptip1;
isolv7 = ( prog (-1.*a0) (-5.*dx1) PAS dx1 (-1.*dx1)
PAS (0.2*dx1) dx1 PAS dx1 (5.*dx1) a0 );
si(graph);
*opti isov lign;
TRAC phi0 (stot0 et lcrack0) isolv7 ;
TRAC psi0 (stot0 et lcrack0) isolv7;
*opti isov surf;
finsi;
*******************************************************
*** Modeles et materiaux
* definition des models elementaires
mod1a = MODE stot0 mecanique elastique isotrope
plastique_endom rousselier xq4R;
* enrichissement
Che1X = tabl;
Che1X . 0 = TRIE mod1a psi0 phi0;
* apres cela on assemble
mod1tot = mod1a ;
* on crée le ou les materiau et on les assemble
Ey1 = 5.520E+2 '/' 2.627E-3 ;
nu1 = 0.3 ;
rho1= 7.8e-6 ;
Epsprog0 = (PROG
0.000E+0 2.627E-3 2.967E-3 3.276E-3
3.585E-3 5.129E-3 8.221E-3 1.132E-2 1.442E-2
1.754E-2 3.334E-2 5.006E-2 6.914E-2 9.391E-2
1.314E-1 1.956E-1 3.134E-1 5.336E-1) ;
Sigprog0 = (PROG
0.000E+0 5.520E+2 5.530E+2 5.540E+2
5.550E+2 5.600E+2 5.700E+2 5.800E+2 5.900E+2
6.000E+2 6.500E+2 7.000E+2 7.500E+2 8.000E+2
8.500E+2 9.000E+2 9.500E+2 1.000E+3);
Tcurve = EVOL 'MANU' 'Epsilon' Epsprog0 'Sigma' Sigprog0 ;
mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
*******************************************************
*** Matrices + CL
*Rig1tot = RIGI mod1tot mat1tot;
* on bloque les depl de corps rigides
pgoup1 = stot0 POIN 'PROC' (0. hgoup);
pgoup2 = stot0 POIN 'PROC' (0. (-1.*hgoup));
pgouptot = pgoup1 et pgoup2;
cl1x = BLOQ 'UX' (pgouptot);
cl1y = BLOQ 'UY' (pgouptot);
cl1tot = cl1x et cl1y;
*******************************************************
*** deplacement impose ***
* evolution
du0 = 0.04 ;
ufin0 = 0.44 ;
si logdech;
uval0 = PROG 0. 'PAS' du0 ufin0 'PAS' (0. - du0) (0.8*ufin0);
tval0 = (PROG 0. 'PAS' du0 ufin0 'PAS' du0 (1.2*ufin0))* 1.E-2 ;
sinon;
uval0 = PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0;
tval0 = (PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0)* 1.E-2;
finsi;
Nmax0 = DIME tval0 ;
uev0 = EVOL 'ROUG' 'MANU' 't' tval0 'u' uval0 ;
* chpoint
uchp0 = (MANU 'CHPO' pgoup1 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS')
ET (MANU 'CHPO' pgoup2 2 'UX' 0. 'UY' -1. 'NATU' 'DIFFUS');
fchp0 = DEPI cl1tot uchp0;
* chargement
uimp1 = CHAR 'DIMP' fchp0 uev0 ;
tdess1 = tabl;
tdess1 . 1 = mot 'MARQ PLUS';
si(graph);
dess uev0 tdess1;
trac (vect uchp0 'DEPL' 'VERT') (stot0 et lcrack0);
finsi;
*******************************************************
*** RESOlution avec PASaPAS
* Initialisation calcul
TAB1 = TABL;
TAB1.'MODELE' = mod1tot;
TAB1.'CARACTERISTIQUES' = mat1tot;
TAB1.'BLOCAGES_MECANIQUES' = cl1tot;
TAB1.'CHARGEMENT' = uimp1;
TAB1.'PRECISION' = 1E-05;
TAB1.'PRECISINTER' = 1E-05;
*TAB1. 'TEMPS0' = 0. ;
*TAB1. 'PROCESSEURS' = 'MONO';
TAB1.'TEMPS_CALCULES' = tval0 ;
TAB1.'TEMPS_SAUVES' = tval0 ;
* Resolution non-lineaire
PASAPAS TAB1;
****************
* Cas fissure maillée :
sbtot = satot;
* definition des models elementaires
mod1b = MODE sbtot mecanique elastique isotrope
plastique_endom rousselier ;
mat1b = MATE mod1b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
* on bloque les depl de corps rigides
dcl0 = d11d et d1d;
cl0 = BLOQ 'UY' dcl0;
pgoup1b = sbtot POIN 'PROC' (0. hgoup);
cl1xb = BLOQ 'UX' (pgoup1b);
cl1yb = BLOQ 'UY' (pgoup1b);
cl1totb = cl1xb et cl1yb ;
cltot = cl1totb et cl0;
*** deplacement impose ***
* chpoint
uchp0b = (MANU 'CHPO' pgoup1b 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS');
fchp0b = DEPI cl1totb uchp0b;
* chargement
uimp1b = CHAR 'DIMP' fchp0b uev0 ;
tdess1 = tabl;
tdess1 . 1 = mot 'MARQ PLUS';
si(graph);
dess uev0 tdess1;
trac (vect uchp0b 'DEPL' 'VERT') (sbtot et lcrack0);
finsi;
*******************************************************
* Initialisation calcul
TAB2 = TABL;
TAB2.'MODELE' = mod1b;
TAB2.'CARACTERISTIQUES' = mat1b;
TAB2.'BLOCAGES_MECANIQUES' = cltot;
TAB2.'CHARGEMENT' = uimp1b;
TAB2.'PRECISION' = 1E-05;
TAB2.'PRECISINTER' = 1E-05;
*TAB2. 'PROCESSEURS' = 'MONO';
TAB2.'TEMPS_CALCULES' = tval0 ;
TAB2.'TEMPS_SAUVES' = tval0 ;
* Resolution non-lineaire
PASAPAS TAB2;
******************************************************
* post traitement de pasapas
* courbes force ouverture ****************************
cod1 = prog ; for1 = prog;
cod2 = prog ; for2 = prog;
i = -1;
repe BPOST1 (Nmax0);
i = i + 1;
u1i = TAB1 . 'DEPLACEMENTS' . i;
cod1i = extr u1i 'UY' p4;
cod1 = cod1 et (prog cod1i);
frea1i = TAB1 . 'REACTIONS' . i;
for1i = extr frea1i 'FY' pgoup1;
for1 = for1 et (prog for1i);
u2i = TAB2 . 'DEPLACEMENTS' . i;
cod2i = extr u2i 'UY' p4;
cod2 = cod2 et (prog cod2i);
frea2i = TAB2 . 'REACTIONS' . i;
for2i = extr frea2i 'FY' pgoup1b;
for2 = for2 et (prog for2i);
fin BPOST1;
evfcod1 = EVOL 'BLEU' 'MANU' 'ouverture' cod1 'force' for1 ;
evfcod2 = EVOL 'ROUG' 'MANU' 'ouverture' cod2 'force' for2 ;
si(graph);
dess (evfcod1 et evfcod2);
finsi;
* tracés a la fin de la mise en charge **************
sig1i = TAB1 . 'CONTRAINTES' . i;
var1i = TAB1 . 'VARIABLES_INTERNES' . i;
svm1 = VMIS sig1i mod1tot;
sig2i = TAB2 . 'CONTRAINTES' . i;
var2i = TAB2 . 'VARIABLES_INTERNES' . i;
svm2 = VMIS sig2i mod1b;
si(graph);
TRAC svm1 mod1tot (DEFO u1i (stot0 et lcrack0) 3.);
TRAC svm2 mod1b (DEFO u2i (sbtot et lcrack0) 3.);
TRAC var1i mod1tot (DEFO u1i stot0 3.);
TRAC var2i mod1b (DEFO u2i sbtot 3.);
finsi;
* Reconstruction du déplacement réel pour avoir la déformée :
chp2UX = XFEM 'RECO' tab1.deplacements. i mod1tot ;
si graph;
stot0d = (defo (coul stot0 roug) chp2UX) ;
sbtotd = defo (coul sbtot jaun) tab2.deplacements. i;
trac stot0d;
trac (stot0d et (defo stot0 tab1.deplacements. i));
finsi;
*******************************************************
* propagation élémentaire de fissure
* ... a venir
*==============================================================*
*
err6 = (ABS (((maxi for2) - (maxi for1))/(maxi for2)));
si logmess;
mess 'erreur sur l effort : ' err6;
finsi;
SI (err6 >EG 1.);
ERRE 5;
FINSI;
*opti donn 5 ;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer