* fichier : xfem02.dgibi
********************************************************************
*
* xfem02.dgibi
*
* calcul avec elements xfem (XQ4R)
* d'une plaque elasto plastique en traction
* avec fissure droite
*
* création : bp, le 24.09.2009
* ajout g_theta : bp, 08.06.2010
*
********************************************************************
*opti echo 1 ;
*graph = vrai; opti trac PSC ; opti EPTR 5;
graph = faux;
fllong = faux;
*******************************************************
*** Options de calcul : 2D
opti dime 2 elem qua4 mode plan defo;
*******************************************************
*** Maillage : Plaque rectangulaire saine
* dimensions
* bord droit, gauche, hauteur
b= 0.050 ; bg = -0.25 * b ; hb1= 0.6 * b ;
hgoup= 0.35 * b;
* finesse
* ici, nx1 doit etre un multiple de 2 et de 5
nx1 = 20;
si(fllong); nx1 = 40; fins;
dx1 = (b-bg) / (flot nx1) ;
na1 = enti (b / dx1);
na6 = nx1 - na1;
na2 = 2;
ha1 = (flot na2) * dx1;
* demi-zone centrale
DENS (dx1) ;
pa1 = (0.) (0.*ha1) ;
pa2 = (b ) (0.*ha1) ;
pa3 = (b) (ha1) ;
pa4 = (0.) (ha1) ;
pa5 = (bg) (ha1) ;
pa6 = (bg) (0.*ha1) ;
da1 = pa1 droi na1 pa2;
da2 = pa2 droi na2 pa3;
da3 = pa3 droi na1 pa4;
da4 = pa4 droi na6 pa5;
da5 = pa5 droi na2 pa6;
da6 = pa6 droi na6 pa1;
da61 = (da6 et da1);
da34 = (da3 et da4);
sa1 = dall da61 da2 da34 da5;
* zone de transition
dx2 = (2.*dx1);
DENS (dx2);
pb3 = pa3 plus (0. dx2);
pb4 = pa4 plus (0. dx2);
pb5 = pa5 plus (0. dx2);
nb3 = na1 / 2;
nb4 = na6 / 2;
db3 = pb3 droi (nb3) pb4;
db4 = pb4 droi (nb4) pb5;
db34 = db3 et db4;
nb34 = nb3 + nb4;
* procedure maillage 2 pour 1 entre 2 lignes : lgros et lfine
lfine = da34;
lgros = db34;
nbmotif1= nbel lgros;
if1 = 1; ig1 = 1;
logpm2 = vrai;
repe bmotif1 nbmotif1;
pf1 = lfine POIN if1;
if1 = if1 + 1;
pf3 = lfine POIN (if1);
if1 = if1 + 1;
pf2 = lfine POIN (if1);
pg1 = lgros POIN ig1;
ig1 = ig1 + 1;
pg2 = lgros POIN (ig1);
pg3 = 0.5 * (pg1 plus pg2);
pm3 = 0.5 * (pf3 plus pg3);
si(logpm2);
logpm2 = faux;
pm2 = 0.5 * (pf2 plus pg2);
geo1 = (manu 'QUA4' pf1 pf3 pm3 pg1)
et (manu 'QUA4' pf3 pf2 pm2 pm3)
et (manu 'QUA4' pm3 pm2 pg2 pg1);
sino;
logpm2 = vrai;
geo1 = (manu 'QUA4' pf1 pf3 pm3 pm2)
et (manu 'QUA4' pf3 pf2 pg2 pm3)
et (manu 'QUA4' pm2 pm3 pg2 pg1);
finsi;
si(&bmotif1 ega 1); motif1 = geo1;
sino; motif1 = motif1 et geo1;
fins;
fin bmotif1;
* on recupere motif1
sb1 = motif1;
* zone "grossiere"
DENS (dx2);
pc3 = (b) (hgoup);
pc5 = (bg) (hgoup);
dc34 = pc3 droi nb34 pc5;
sc1 = regl db34 dc34;
pd3 = ( b) (hb1);
pd5 = (bg) (hb1);
dd34 = pd3 droi nb34 pd5;
sd1 = regl dc34 dd34;
* symetrie et assemblages
sa2 sb2 sc2 sd2 = sa1 sb1 sc1 sd1 SYME 'DROIT' pa6 pa2 ;
* satot = sa1 et sa2;
* sbtot = (sb1 et sc1 et sd1) et (sb2 et sc2 et sd2);
satot = (sa1 et sb1) et (sa2 et sb2);
sbtot = (sc1 et sd1) et (sc2 et sd2);
stot0 = satot et sbtot;
ELIM stot0 (1.E-5*dx1);
mess 'maillage central ext total' ;
mess 'nbno= ' (nbno satot) (nbno stot0);
mess 'nbel= ' (nbel satot) (nbel stot0);
*******************************************************
*** Maillage de la fissure
* longueur de la fissure
a0 = 0.025 ;
* finesse du maillage de la fissure
na0 = 10;
da0 = a0 / (flot na0) ;
* maillage
DENS (da0) ;
ptip0 = (bg) (0.);
ptip1 = (a0) (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 satot 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 satot mecanique elastique isotrope plastique xq4R;
mod1b = MODE sbtot mecanique elastique isotrope;
* enrichissement
Che1X = tabl;
Che1X . 0 = TRIE mod1a psi0 phi0;
* apres cela on assemble
mod1tot = mod1a et mod1b;
* on crée le ou les materiau et on les assemble
Ey1 = 5.520E+8 '/' 2.627E-3 ;
nu1 = 0.3 ;
rho1= 7800. ;
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) * 1.E6 ;
Tcurve = EVOL 'MANU' 'Epsilon' Epsprog0 'Sigma' Sigprog0 ;
mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve;
*******************************************************
*** CL et chargement
* on bloque les depl de corps rigides
pgoup1 = sbtot POIN 'PROC' (0. hgoup);
pgoup2 = sbtot 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.10E-3 ;
*ufin0 = 0.5E-3 ;
ufin0 = 1.E-3 ;
uval0 = PROG 0. 'PAS' du0 ufin0 ;
tval0 = uval0 * 1.E-2 ;
dtval0 = du0 * 1.E-2 ;
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;
*******************************************************
*** MISE EN CHARGE AVEC PASAPAS ***
* temps de mise en charge
tcalc0 = PROG 0. 'PAS' dtval0 (3.*dtval0);
Ncalc0 = DIME tcalc0;
* Initialisation de la table
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' = tcalc0 ;
TAB1.'TEMPS_SAUVES' = tcalc0 ;
* Resolution non-lineaire ******
PASAPAS TAB1;
* Mise a Jour de la table TXFEM (a terme devra etre inclus dans PASAPAS)
i = 0 ;
repe BOUX0 (Ncalc0 - 1);
i = i + 1;
TXFEM . i = TXFEM . (i - 1);
Che1X . i = Che1X . (i - 1);
fin BOUX0;
* calcul de l integrale J *****
* * J local
* si (fllong);
* SUPTAB = TABL ;
* SUPTAB.'OBJECTIF' = MOT 'J';
* SUPTAB.'PSI' = psi0;
* SUPTAB.'PHI' = phi0;
* SUPTAB.'FRONT_FISSURE' = ptip1 ;
* SUPTAB.'SOLUTION_PASAPAS' = TAB1;
* SUPTAB.'COUCHE' = 2;
* G_THETA SUPTAB;
* si(graph);
* q7 = (SUPTAB . 'CHAMP_THET1') ;
* vq7 = VECT q7 AZUR ;
* trac vq7 (stot0 et lcrack0);
* fins;
* fins;
* calcul de l integrale J global
JTAB = TABL;
JTAB.'OBJECTIF' = MOT 'J';
JTAB.'PSI' = psi0;
JTAB.'PHI' = phi0;
JTAB.'FRONT_FISSURE' = ptip1 ;
JTAB.'SOLUTION_PASAPAS' = TAB1;
*definition dun champ_theta perso pour calcul de j global
sared = ELEM satot 'APPUYE' 'LARGEMEMENT' (CONT satot);
sared = (satot DIFF sared) coul vert;
monch1 = MANU 'CHPO' sared 'UX' 1. 'UY' 0.;
monvec1 = vect monch1 vert;
si(graph);
trac monvec1 (stot0 et lcrack0);
* trac (exco 'UX' monch1) (stot0 et lcrack0);
fins;
JTAB.'CHAMP_THETA' = monch1;
G_THETA JTAB;
*******************************************************
*** PROPAGATION élémentaire de fissure ***
*rem: + tard, prevoir procedure surtout pour le 3d
NBPROPA = 2;
da1 = 1.5 * dx1;
repe BPROPA NBPROPA;
i = i + 1;
*** Maillage de la fissure ***
ptip1 = ptip1 plus (da1 0.);
lcrack1 = (lcrack0 droi ptip1) coul 'ROUG';
TXFEM . i = tabl;
TXFEM . i . 'POINTE1' = ptip1;
TXFEM . i . 'FISSURE' = lcrack1;
*** Actualisations ***
* des level set
psi1 phi1 = PSIPHI satot lcrack1 'DEUX' ptip1 ;
si(graph);
*opti isov lign;
trac (stot0 et lcrack1);
TRAC phi1 (stot0 et lcrack1) isolv7 ;
TRAC psi1 (stot0 et lcrack1) isolv7;
*opti isov surf;
finsi;
* du modele, de la rigidite,... (inutile pour le materiau)
Che1X . i = TRIE mod1a psi1 phi1;
mod1tot = mod1a et mod1b;
TAB1.'MODELE' = mod1tot;
*** resolution non-lineaire PASaPAS ***
tlast1 = extr TAB1.'TEMPS_CALCULES' (dime TAB1.'TEMPS_CALCULES');
tcalc1 = TAB1.'TEMPS_CALCULES' et (PROG (tlast1 + dtval0));
Ncalc1 = dime tcalc1;
TAB1.'TEMPS_CALCULES' = tcalc1 ;
TAB1.'TEMPS_SAUVES' = tcalc1 ;
PASAPAS TAB1;
*** calcul de l integrale J global (le contour ne change pas) ***
JTAB.'PSI' = psi1;
JTAB.'PHI' = phi1;
JTAB.'FRONT_FISSURE' = ptip1 ;
*JTAB.'CHAMP_THETA' = JTAB.'CHAMP_THET1';
JTAB.'CHAMP_THETA' = monch1;
* on utilise la reprise car le contour ne change pas
G_THETA JTAB;
* si (fllong);
* *** calcul de l integrale J local, de K1, de K2 ***
* * a chaque propa il faut une nouvelle table, car le contour change...
* SUPTAB = TABL ;
* SUPTAB.'OBJECTIF' = MOT 'J';
* SUPTAB.'PSI' = psi1;
* SUPTAB.'PHI' = phi1;
* SUPTAB.'FRONT_FISSURE' = ptip1 ;
* SUPTAB.'SOLUTION_PASAPAS' = TAB1;
* SUPTAB.'COUCHE' = 2;
*
* G_THETA SUPTAB;
* q7 = (SUPTAB . 'CHAMP_THET1') ;
*
* si(graph);
* vq7 = VECT q7 AZUR ;
* trac vq7 (stot0 et lcrack1);
* fins;
* fins;
* i etant l indice de la table de PASAPAS TAB1,
* il faut l'incrementer avant et apres chaque appel a pasapas
i = i + 1;
TXFEM . i = TXFEM . (i - 1);
Che1X . i = Che1X . (i - 1);
fin BPROPA;
******************************************************
*** POST TRAITEMENT ***
* courbe J global - t (le contour ne change pas) ***
si(graph);
dess JTAB . EVOLUTION_RESULTATS 'TITR' 'Jglobal - t';
finsi;
* courbe force ouverture ****************************
cod1 = prog;
for1 = prog;
ipost = -1;
Npost = i;
repe BPOST1 (Npost + 1);
ipost = ipost + 1;
u1i = TAB1 . 'DEPLACEMENTS' . ipost;
cod1i = extr u1i 'AY' pa1;
cod1 = cod1 et (prog cod1i);
frea1i = TAB1 . 'REACTIONS' . ipost;
for1i = extr frea1i 'FY' pgoup1;
for1 = for1 et (prog for1i);
fin BPOST1;
evfcod1 = EVOL 'VERT' 'MANU' 'ouverture' cod1 'force' for1 ;
si(graph);
dess evfcod1 tdess1 'TITR' 'Force - Ouverture';
finsi;
* quelques tracés ****************************
u1i = TAB1 . 'DEPLACEMENTS' . i;
u1phy = XFEM 'RECO' u1i mod1tot ;
var1i = TAB1 . 'VARIABLES_INTERNES' . i;
sig1i = TAB1 . 'CONTRAINTES' . i;
def1i = ELAS mod1tot sig1i mat1tot;
si(graph);
TRAC var1i mod1tot (DEFO u1phy stot0 20.) 'TITR' 'EPSE';
TRAC (exco EPYY def1i) mod1tot (DEFO u1phy stot0 20.) 'TITR' 'EPYY';
finsi;
*
*opti donn 5 ;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer