* fichier : xfem3d_02.dgibi
*******************************************************
**-------------------------------------------------------------------
MESS '' ;
MESS ' CALCUL ELASTIQUE DYNAMIQUE ';
MESS ' en 3D XFEM ' ;
MESS ' Unites: N MPa mm g ms' ;
MESS '' ;
**--------------------------------------------------------------------
*******************************************************
*** Options de calcul
opti dime 3 elem CUB8 mode trid;
* ajout de option epsilon lineaire pour la precision des test!
OPTION epsilon lineaire;
opti 'TRAC' PSC ;
* opti 'TRAC' X ;
*******************************************************
*** options de trace
isov1 = prog -0.7 pas 0.1 0.7;
*******************************************************
*** maillage
pa0 = 0. 0. 0. ; pb0 = 3. 0. 0. ;
liab = pa0 droi 12 pb0;
su1a = liab trans 16 (0. 0. -4.);
liam = COTE 3 su1a;
su1b = liam trans 8 (0. 0. -2.);
liah = COTE 3 su1b;
vol1a = su1a volu tran 12 (0. 6. 0.);
vol1b = su1b volu tran 12 (0. 6. 0.);
vol1 = vol1a et vol1b;
elim 0.001 vol1;
*
*trac cach (vol1 et (liab coul vert) et (liah coul rouge));
*mess (nbno vol1) (nbel vol1) ;
*
pc1 = (enve vol1a) poin proch (3. 6. 0.);
pc2 = (enve vol1b) poin proch (3. 6. -6.);
**********************************************************
*** fissure 1 = debouchante (=1 seule pointe) rectiligne
opti elem tri3;
a0 = 3.;
pf1 = 0. 0. -3.;
pf2 = 3. 0. -3.;
pf3 = 3. a0 -3;
pf4 = 0. a0 -3;
ll12 = droi 8 pf1 pf2;
fro23 = droi 10 pf2 pf3;
pcrack2 = droi 8 pf3 pf4;
ll41 = droi 10 pf4 pf1;
crack1 = (dall ll12 fro23 pcrack2 ll41 'PLAN') coul roug;
psi1 phi1 = PSIPHI vol1a crack1 'DEUX' pcrack2;
*isov1 = prog -0.001 pas 0.0001 0.001;
* opti isov surf ;
* trac phi1 (vol1a et vol1b) isov1 ((vol1a et vol1b ) et crack1);
* trac psi1 (vol1a et vol1b) isov1 ((vol1a et vol1b ) et crack1);
**********************************************************
*** MODELE & MATERIAU ***
**********************************************************
Ey0 = 200000.;
nu0 = 0.3 ;
rho0= 7.8e-3 ;
**********************************************************
*** Modele & materiau ***
*
* elements standards
mod1b = 'MODELISER' vol1b 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'CUB8';
mat1b = 'MATERIAU' mod1b'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
*
**********************************************************
* zone de propagation (X-FEM)
mod1a = 'MODELISER' vol1a 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'XC8R';
* che1X = TRIE Mod1a psi1 phi1 'SAUT' ;
che1X = TRIE Mod1a psi1 phi1 ;
*
modtot = mod1a et mod1b;
mat1a = 'MATERIAU' mod1a 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
**********************************************************
**********************************************************
*** CL et DEPLACEMENTS IMPOSES ***
**********************************************************
cl1 = bloq (pc1 et pc2) 'UZ';
cl2 = bloq liab 'UX' 'UY';
cl3 = bloq liah 'UX' 'UY';
cltot = (cl1 'ET' cl2 'ET' cl3) ;
*
*... chargement
for1 = FORC 'FZ' liah (-3.E3);
for2 = FORC 'FZ' liab (3.E3);
*
xx1 = prog 0. 7.5e-4 10.;
yy1 = prog 0. 1. 1.;
evochar = evol manu 'temps' xx1 'charge' yy1;
char1 = char meca (for1 et for2) evochar;
*
*trac (vect (for1 et for2) 'FX' 'FY' 'FZ' 1.e-2) ((enve vol1)
*et (liah coul vert) et (liab coul rose));
*******************
**********************************************************
*** MATRICES ***
**********************************************************
*** de rigidite ***
* Xfem
*Rig1a = RIGI Mod1a mat1a;
* standard
*Rig1b = RIGI Mod1b mat1b;
*** de masse ***
* Xfem
mas1a = mass mod1a mat1a;
* standard
mas1b = mass mod1b mat1b ;
*
*=======================================================*
* *
* RESOLUTION *
* *
*=======================================================*
dt1 = 2.5e-4;
xx0 = prog 0. pas dt1 'NPAS' 4;
*
tab1 = table;
tab1.dynamique = vrai;
tab1.modele = modtot;
tab1.caracteristiques = mat1a et mat1b;
tab1.blocages_mecaniques = cltot;
tab1.chargement = char1;
tab1.temps_calcules = xx0;
pasapas tab1;
*
*post traitement
ind0 = index (tab1.deplacements);
nnn1 = dime ind0;
lwela1 = prog;
lwcin1 = prog;
ltemp1 = prog;
ev1 = table;
def1 = table;
*
ii = 0 ;
repeter bouc1 nnn1;
ii = ii + 1;
pass = ind0.ii;
sig1 = tab1.contraintes.pass;
ut1 = tab1.deplacements.pass;
vt1 = tab1.vitesses.pass;
tt1 = tab1.temps.pass;
etot1 = EPSI ut1 modtot;
*
ev1.ii = evol chpo ut1 'UZ' liah;
def1.ii = defo ut1 vol1 30.;
ltemp1 = ltemp1 et (prog tt1);
*...elastique
Wela0 = 'INTG' (0.5 *(ENER sig1 etot1 modtot)) modtot ;
lwela1 = lwela1 et (prog wela0);
*
*...cinetique
Wcin0 = 0.5 * (XTMX vt1 (mas1a 'ET' mas1b)) ;
lwcin1 = lwcin1 et (prog wcin0);
fin bouc1;
*
*
**********************************************************
errmax = 0.;
lerr1= prog 0. -1.20735e-07 2.44871e-06 2.08117e-06
-9.9959e-07 3.38064e-04 9.78166e-04;
lresu = prog 1.e-8 3.49142e-03 1.09155e-02 1.8968e-02 2.44906e-02
73.8862 27.5316;
*comparaison XFEM et STANDARD 3D
ii = 0;
repeter bouc3 nnn1;
ii = ii + 1 ;
ord2 = extr lresu ii;
val0 = extr lerr1 ii;
ord1 = extr ev1.ii 'ORDO';
val1 =( (abs(extr ord1 3)) - ord2)/ ord2;
message 'ord1 ord2 val1' (extr ord1 3) ord2 val1;
si (ii > 1);
message ' valeurs a comparer' val1 val0;
errmax = maxi(prog errmax val1);
finsi;
fin bouc3;
*
evn1 = evol manu 'temps' ltemp1 'ener' lwela1;
evn2 = evol manu 'temps' ltemp1 'ener' lwcin1;
ii = ii + 1 ;
ord2 = extr lresu ii;
val0 = extr lerr1 ii;
ord1 = extr evn1 'ORDO';
val1 = abs(( (abs(extr ord1 5)) - ord2)/ ord2);
message 'ord1 ord2 val1' (extr ord1 5) ord2 val1;
message ' valeurs a comparer' val1 val0;
errmax = maxi(prog errmax val1);
*
ii = ii + 1 ;
ord2 = extr lresu ii;
val0 = extr lerr1 ii;
ord1 = extr evn2 'ORDO';
val1 = abs(( (abs(extr ord1 5)) - ord2)/ ord2);
message 'ord1 ord2 val1' (extr ord1 5) ord2 val1;
message ' valeurs a comparer' val1 val0;
errmax = maxi(prog errmax val1);
*
si(errmax >EG 0.001);
ERRE 5;
finsi;
*
*opti donn 5;
fin;
**-------------------------------------------------------------------
**-------------------------------------------------------------------
MESS '' ;
MESS ' CALCUL ELASTIQUE DYNAMIQUE ';
MESS ' en 3D standard ' ;
MESS ' Unites: N MPa mm g ms' ;
MESS '' ;
**--------------------------------------------------------------------
*
opti echo 1 ;
*******************************************************
*** Options de calcul
opti dime 3 elem CUB8 mode trid;
* opti 'TRAC' PSC ;
opti 'TRAC' X ;
*******************************************************
*** options de tracé
*isov1 = prog -0.7 pas 0.1 0.7;
*******************************************************
*** maillage 1
pa1 = 0. 0. 0. ; pb1 = 3. 0. 0. ;
liab2 = pa1 droi 12 pb1;
*
sy1a = liab2 trans 12 (0. 0. -3.);
bor1 = COTE 3 sy1a;
vol2a = sy1a volu tran 12 (0. 6. 0.);
*
pa2 = 0. 0. -3. ; pb2 = 3. 0. -3. ;
bor2 = pa2 droi 12 pb2;
sy1b = bor2 trans 12 (0. 0. -3.);
vol2b = sy1b volu tran 12 (0. 6. 0.);
*
liah2 = inve ( COTE 3 sy1b);
*
vol2 = vol2a et vol2b;
*
*longueur fissure 3 mm
xfro1 = 2.8;
*
px = poin vol2a plan (1. 0. -3.) (0. 1. -3.) (0. 0. -3.) 0.001;
pla1 = elem (enve vol2a) appuye strictement px;
elim 0.001 (pla1 et bor1);
cha1 = 'COOR' 2 pla1;
geo1 = cha1 'POIN' 'EGSUPE' xfro1;
pla2 = elem pla1 appuye strictement geo1;
*
px = poin vol2b plan (1. 0. -3.) (0. 1. -3.) (0. 0. -3.) 0.001;
plb1 = elem (enve vol2b) appuye strictement px;
elim 0.01 (plb1 et bor2);
chb1 = 'COOR' 2 plb1;
geo2 = chb1 'POIN' 'EGSUPE' xfro1;
plb2 = elem plb1 appuye strictement geo2;
*
**********************************************************
*** fissure 1 = debouchante (=1 seule pointe) rectiligne
pf1 = (cont pla2) poin proch (0. 3. -3.);
pf2 = (cont pla2) poin proch (3. 3. -3.);
*
pcrack2 = (cont pla2) elem compris pf1 pf2;
crack1 = (diff pla1 pla2);
*
**********************************************************
*** Modele
*
Mod2a = vol2a MODE mecanique elastique isotrope ;
Mod2b = vol2b MODE mecanique elastique isotrope;
*
**********************************************************
Ey1 = 200000. ;
nu1 = 0.3 ;
rho1= 7.8e-3 ;
mat2a = MATE Mod2a 'YOUN' Ey1 'NU' nu1 'RHO' rho1 ;
mat2b = MATE Mod2b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 ;
**********************************************************
*CL
pc3 = (enve vol2) poin proch (3. 6. 0.);
pc4 = (enve vol2) poin proch (3. 6. -6.);
*
cl11 = bloq (pc3 et pc4) 'UZ' ;
cl12 = bloq liab2 'UX' 'UY';
cl13 = bloq liah2 'UX' 'UY';
clto2 = cl11 et cl12 et cl13;
*
nb1 = nbno pla2;
pla3 = chan 'POI1' pla2;
ii = 1;
pp0 = pla3 poin ii;
plb3 = plb2 poin proch pp0;
repeter bouc1 (nb1 - 1);
ii = ii + 1 ;
pp0 = pla3 poin ii;
plb3 = plb3 et (plb2 poin proch pp0);
fin bouc1;
*
cond1 = 'RELA' 'UX' pla3 - 'UX' plb3;
cond2 = 'RELA' 'UY' pla3 - 'UY' plb3;
cond3 = 'RELA' 'UZ' pla3 - 'UZ' plb3;
**********************************************************
*CHARGEMENT
fo11 = FORC 'FZ' liah2 (-3.E3);
fo12 = FORC 'FZ' liab2 (3.E3);
*
*trac (vect (fo11 et fo12) 'FX' 'FY' 'FZ' 1.e-2) ((enve vol2)
*et (liah2 coul vert) et (liab2 coul rose));
*opti donn 5;
*
xx2 = prog 0. 7.5e-4 10.;
yy2 = prog 0. 1. 1.;
evocha2 = evol manu 'temps' xx2 'charge' yy2;
char2 = char meca (fo11 et fo12) evocha2;
*
**********************************************************
*** MATRICES ***
**********************************************************
*** de masse ***
mas2a = mass mod2a mat2a;
mas2b = mass mod2b mat2b ;
**********************************************************
*=======================================================*
* *
* RESOLUTION *
* *
*=======================================================*
dt1 = 2.5e-4;
xx0 = prog 0. pas dt1 'NPAS' 4;
*
tab2 = table;
tab2.dynamique = vrai;
tab2.modele = mod2a et mod2b;
tab2.caracteristiques = mat2a et mat2b;
tab2.blocages_mecaniques = clto2 et cond1 et cond2 et cond3;
tab2.chargement = char2;
tab2.temps_calcules = xx0;
pasapas tab2;
*post traitement
ind0 = index (tab2.deplacements);
nnn2 = dime ind0;
*
lwela2 = prog;
lwcin2 = prog;
ltemp2 = prog;
ev2 = table;
def2 = table;
*
ii = 0;
repeter bouc2 nnn2;
ii = ii + 1;
pass = ind0.ii;
sig2 = tab2.contraintes.pass;
ut2 = tab2.deplacements.pass;
vt2 = tab2.vitesses.pass;
tt2 = tab2.temps.pass;
etot2 = EPSI ut2 (mod2a et mod2b);
*
ev2.ii = evol rose chpo ut2 'UZ' liah2;
def2.ii = defo ut2 vol2 30. rose;
ltemp2 = ltemp2 et (prog tt2);
*...elastique
Wela2 = 'INTG' (0.5 *(ENER sig2 etot2 (mod2a et mod2b)))
(mod2a et mod2b) ;
lwela2 = lwela2 et (prog wela2);
*
*...cinetique
Wcin2 = 0.5 * (XTMX vt2 (mas2a 'ET' mas2b)) ;
lwcin2 = lwcin2 et (prog wcin2);
*
fin bouc2;
**********************************************************
errmax = 0.;
lerr1= prog 0. 4.28804e-08 3.46393e-7 2.36123e-6 1.3508e-5 6.27995e-5;
*comparaison
ii = 0;
repeter bouc3 nnn2;
ii = ii + 1 ;
*comparaison XFEM et STANDARD 3D
val0 = extr lerr1 ii;
*
*trac cach (def1.ii et def2.ii);
*
*dess (ev1.ii et ev2.ii);
*
ord1 = extr ev1.ii 'ORDO';
ord2 = extr ev2.ii 'ORDO';
si (abs(extr ord2 3) < 1.e-8);
val1 = val0;
message 'ord1 ord2 val1' (extr ord1 3) (extr ord2 3) val1;
sinon;
val1 = ((extr ord1 3) -(extr ord2 3))/(extr ord2 3);
message 'ord1 ord2 val1' (extr ord1 3) (extr ord2 3) val1;
message ' valeurs a comparer' val1 val0;
finsi;
errmax = maxi(prog errmax val1);
*
fin bouc3;
*
evn1 = evol manu 'temps' ltemp1 'ener' lwela1;
evn2 = evol manu 'temps' ltemp1 'ener' lwcin1;
evn3 = evol rose manu 'temps' ltemp2 'ener' lwela2;
evn4 = evol rose manu 'temps' ltemp2 'ener' lwcin2;
*dess (evn1 et evn3);
*dess (evn2 et evn4);
val0 = 3.37828e-4;
ord1 = extr evn1 'ORDO';
ord2 = extr evn3 'ORDO';
val1 = abs (((extr ord1 5) -(extr ord2 5))/(extr ord2 5));
message 'ord1 ord2 val1' (extr ord1 5) (extr ord2 5) val1;
message ' valeurs a comparer' val1 val0;
errmax = maxi(prog errmax val1);
*
val0 = 9.78566e-4;
ord1 = extr evn2 'ORDO';
ord2 = extr evn4 'ORDO';
val1 = ((extr ord1 5) -(extr ord2 5))/(extr ord2 5);
message 'ord1 ord2 val1' (extr ord1 5) (extr ord2 5) val1;
message ' valeurs a comparer' val1 val0;
errmax = maxi(prog errmax val1);
si(errmax >EG 0.001);
ERRE 5;
finsi;
*
opti donn 5;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer