Télécharger la plaquette Cast3M


fr en

Télécharger xfem02.dgibi

Retour à la liste

Numérotation des lignes : 

  1. ************************************************************************
  2. ************************************************************************
  3. ********************************************************************
  4. * 
  5. *              xfem02.dgibi
  6. *
  7. * calcul avec elements xfem (XQ4R) 
  8. * d'une plaque elasto plastique en traction
  9. * avec fissure droite 
  10. * 
  11. * création : bp, le 24.09.2009 
  12. * ajout g_theta : bp, 08.06.2010
  13. *
  14. ********************************************************************
  15. *opti echo 1 ;
  16. *graph = vrai; opti trac PSC ; opti EPTR 5;
  17. graph = faux;
  18. fllong = faux;
  19.  
  20. *******************************************************
  21. *** Options de calcul : 2D
  22. opti dime 2 elem qua4 mode plan defo;
  23.  
  24.  
  25. *******************************************************
  26. *** Maillage : Plaque rectangulaire saine
  27.  
  28. * dimensions
  29. * bord droit, gauche, hauteur 
  30.   b= 0.050 ;     bg = -0.25 * b ; hb1= 0.6 * b ;
  31.   hgoup= 0.35 * b;
  32.  
  33. * finesse
  34. * ici, nx1 doit etre un multiple de 2 et de 5
  35. nx1 = 20;
  36. si(fllong); nx1 = 40; fins;
  37. dx1 = (b-bg) / (flot nx1) ;
  38. na1 = enti (b / dx1); 
  39. na6 = nx1 - na1;
  40. na2 = 2;
  41. ha1 = (flot na2) * dx1;
  42.  
  43. * demi-zone centrale
  44. DENS (dx1) ;
  45. pa1 = (0.) (0.*ha1) ;
  46. pa2 = (b ) (0.*ha1) ;
  47. pa3 = (b)  (ha1) ;
  48. pa4 = (0.) (ha1) ;
  49. pa5 = (bg) (ha1) ;
  50. pa6 = (bg) (0.*ha1) ;
  51. da1 = pa1 droi na1 pa2;
  52. da2 = pa2 droi na2 pa3;
  53. da3 = pa3 droi na1 pa4;
  54. da4 = pa4 droi na6 pa5;
  55. da5 = pa5 droi na2 pa6;
  56. da6 = pa6 droi na6 pa1;
  57. da61 = (da6 et da1);
  58. da34 = (da3 et da4);
  59. sa1 = dall da61 da2 da34 da5;
  60.  
  61. * zone de transition
  62. dx2 = (2.*dx1);
  63. DENS (dx2);
  64. pb3 = pa3 plus (0. dx2);
  65. pb4 = pa4 plus (0. dx2);
  66. pb5 = pa5 plus (0. dx2);
  67. nb3 = na1 / 2;
  68. nb4 = na6 / 2;
  69. db3 = pb3 droi (nb3) pb4;
  70. db4 = pb4 droi (nb4) pb5;
  71. db34 = db3 et db4;
  72. nb34 = nb3 + nb4;
  73.  
  74. * procedure maillage 2 pour 1 entre 2 lignes : lgros et lfine
  75. lfine = da34;
  76. lgros = db34;
  77. nbmotif1= nbel lgros;
  78. if1 = 1;  ig1 = 1;
  79. logpm2 = vrai;
  80. repe bmotif1 nbmotif1;
  81.   pf1 = lfine POIN if1; 
  82.   if1 = if1 + 1;
  83.   pf3 = lfine POIN (if1);
  84.   if1 = if1 + 1;
  85.   pf2 = lfine POIN (if1);
  86.   pg1 = lgros POIN ig1; 
  87.   ig1 = ig1 + 1;
  88.   pg2 = lgros POIN (ig1);
  89.   pg3 = 0.5 * (pg1 plus pg2);
  90.   pm3 = 0.5 * (pf3 plus pg3);
  91.   si(logpm2);  
  92.     logpm2 = faux;
  93.     pm2 = 0.5 * (pf2 plus pg2);
  94.     geo1 = (manu 'QUA4' pf1 pf3 pm3 pg1)
  95.         et (manu 'QUA4' pf3 pf2 pm2 pm3)
  96.         et (manu 'QUA4' pm3 pm2 pg2 pg1);
  97.   sino;
  98.     logpm2 = vrai;
  99.     geo1 = (manu 'QUA4' pf1 pf3 pm3 pm2)
  100.         et (manu 'QUA4' pf3 pf2 pg2 pm3)
  101.         et (manu 'QUA4' pm2 pm3 pg2 pg1);
  102.   finsi;
  103.   si(&bmotif1 ega 1); motif1 = geo1;
  104.   sino;               motif1 = motif1 et geo1;
  105.   fins;
  106. fin bmotif1;
  107. * on recupere motif1
  108. sb1 = motif1;
  109.  
  110. * zone "grossiere"
  111. DENS (dx2);
  112. pc3 = (b) (hgoup);
  113. pc5 = (bg) (hgoup);
  114. dc34 = pc3 droi nb34 pc5;
  115. sc1 = regl db34 dc34;
  116.  
  117. pd3 = ( b) (hb1);
  118. pd5 = (bg) (hb1);
  119. dd34 = pd3 droi nb34 pd5;
  120. sd1 = regl dc34 dd34;
  121.  
  122. * symetrie et assemblages
  123. sa2 sb2 sc2 sd2 = sa1 sb1 sc1 sd1  SYME 'DROIT' pa6 pa2 ;
  124.  
  125. * satot = sa1 et sa2;
  126. * sbtot = (sb1 et sc1 et sd1) et (sb2 et sc2 et sd2);
  127. satot = (sa1 et sb1) et (sa2 et sb2);
  128. sbtot = (sc1 et sd1) et (sc2 et sd2);
  129.  
  130. stot0 = satot et sbtot;
  131. ELIM stot0 (1.E-5*dx1);
  132.  
  133.   mess 'maillage central   ext   total' ;
  134.   mess 'nbno=  ' (nbno satot)  (nbno stot0);
  135.   mess 'nbel=  ' (nbel satot)  (nbel stot0);
  136.  
  137.  
  138. *******************************************************
  139. *** Maillage de la fissure
  140.  
  141. * longueur de la fissure
  142. a0 = 0.025 ;
  143.  
  144. * finesse du maillage de la fissure
  145. na0 = 10;
  146. da0 = a0 / (flot na0) ;
  147.  
  148. * maillage
  149. DENS (da0) ;
  150. ptip0 = (bg) (0.);
  151. ptip1 = (a0) (0.);
  152. lcrack0 = (ptip0 droi na0 ptip1) coul 'ROUG';
  153.  
  154. * table (+ pratique pour gérer la propa)
  155. TXFEM = tabl;
  156. TXFEM . 0 = tabl;
  157. TXFEM . 0 . 'POINTE1' = ptip1;
  158. TXFEM . 0 . 'FISSURE' = lcrack0;
  159.  
  160. * creation des level set 
  161. psi0 phi0 = PSIPHI satot lcrack0 'DEUX'  ptip1;
  162. isolv7 = ( prog (-1.*a0) (-5.*dx1) PAS dx1 (-1.*dx1)
  163.            PAS (0.2*dx1) dx1 PAS dx1 (5.*dx1) a0 );
  164. si(graph); 
  165. *opti isov lign;
  166.   TRAC phi0 (stot0 et lcrack0) isolv7 ;
  167.   TRAC psi0 (stot0 et lcrack0) isolv7;
  168. *opti isov surf;
  169. finsi;
  170.  
  171.  
  172. *******************************************************
  173. *** Modeles et materiaux
  174.  
  175. * definition des models elementaires
  176. mod1a = MODE satot mecanique elastique isotrope plastique xq4R;
  177. mod1b = MODE sbtot mecanique elastique isotrope;
  178.  
  179. * enrichissement
  180. Che1X = tabl;
  181. Che1X . 0 = TRIE mod1a psi0 phi0;
  182. Rel1X = tabl;
  183.  
  184. * constructionsion des blocages des ddl X-fem non actifs dans 
  185. * les éléments de transition.
  186. Rel1x . 0 = RELA mod1a;
  187. * on crée le ou les materiau et on les assemble
  188. Ey1 = 5.520E+8 '/' 2.627E-3 ;
  189. nu1 = 0.3 ;
  190. rho1= 7800. ;
  191.  
  192. Epsprog0 = (PROG
  193.    0.000E+0 2.627E-3  2.967E-3 3.276E-3
  194.    3.585E-3 5.129E-3 8.221E-3 1.132E-2 1.442E-2
  195.    1.754E-2 3.334E-2 5.006E-2 6.914E-2 9.391E-2
  196.    1.314E-1 1.956E-1 3.134E-1 5.336E-1) ;
  197. Sigprog0 = (PROG
  198.    0.000E+0 5.520E+2  5.530E+2 5.540E+2
  199.    5.550E+2 5.600E+2 5.700E+2 5.800E+2 5.900E+2
  200.    6.000E+2 6.500E+2 7.000E+2 7.500E+2 8.000E+2
  201.    8.500E+2 9.000E+2 9.500E+2 1.000E+3) * 1.E6 ;
  202. Tcurve = EVOL 'MANU' 'Epsilon' Epsprog0 'Sigma' Sigprog0 ;
  203. lsm1   = Sigprog0 enle 1 ;
  204. lep1   = (Epsprog0 enle 1) - (lsm1 / ey1) ;
  205. evec1  = evol vert manu eps lep1 sig lsm1 ;
  206. si graph ;
  207.   dess (Tcurve et evec1) titr ' Courbes de traction et d ecrouissage (vert)' ;
  208. fins ;
  209.  
  210. mata = MATE mod1a 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'ECRO' evec1 ;
  211. matb = MATE mod1b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 ;
  212.  
  213. * apres cela on assemble
  214. mod1tot = mod1a et mod1b;
  215. mat1tot = mata et matb ;
  216.  
  217.  
  218. *******************************************************
  219. *** CL et chargement
  220.  
  221. * on bloque les depl de corps rigides
  222. pgoup1 = sbtot POIN 'PROC' (0. hgoup);
  223. pgoup2 = sbtot POIN 'PROC' (0. (-1.*hgoup));
  224. pgouptot = pgoup1 et pgoup2;
  225. cl1x = BLOQ  'UX' (pgouptot);
  226. cl1y = BLOQ  'UY' (pgouptot);
  227.  
  228. cl1tot = cl1x et cl1y;
  229.  
  230. * deplacement impose
  231.  
  232. * evolution
  233. du0 = 0.10E-3 ;
  234. *ufin0 = 0.5E-3  ;
  235. ufin0 = 1.E-3  ;
  236. uval0 = PROG  0.  'PAS' du0 ufin0 ;
  237. tval0 = uval0 * 1.E-2 ;
  238. dtval0 = du0 * 1.E-2 ;
  239. Nmax0 = DIME tval0 ;
  240. uev0  = EVOL 'ROUG' 'MANU' 't' tval0 'u'  uval0  ;
  241.  
  242. * chpoint
  243. uchp0 = (MANU 'CHPO' pgoup1 2 'UX' 0. 'UY'   1. 'NATU' 'DIFFUS')
  244.      ET (MANU 'CHPO' pgoup2 2 'UX' 0. 'UY'  -1. 'NATU' 'DIFFUS');
  245. fchp0 = DEPI cl1tot uchp0;
  246.  
  247. * chargement
  248. uimp1  = CHAR 'DIMP' fchp0 uev0 ;
  249.  
  250. tdess1 = tabl;
  251. tdess1 . 1 = mot 'MARQ PLUS';
  252. si(graph);
  253.   dess uev0 tdess1;
  254.   trac (vect uchp0 'DEPL' 'VERT') (stot0 et lcrack0);
  255. finsi;
  256.  
  257.  
  258. *******************************************************
  259. ***           MISE EN CHARGE AVEC PASAPAS           ***
  260.  
  261. * temps de mise en charge
  262. tcalc0 = PROG  0. 'PAS' dtval0 (3.*dtval0);
  263. Ncalc0 = DIME tcalc0;
  264.  
  265. * Initialisation de la table 
  266. TAB1 = TABL;
  267. TAB1.'MODELE'              = mod1tot;
  268. TAB1.'CARACTERISTIQUES'    = mat1tot;
  269. TAB1.'BLOCAGES_MECANIQUES' = cl1tot et(Rel1x . 0) ;
  270. TAB1.'CHARGEMENT'          = uimp1;
  271. TAB1.'PRECISION' = 1E-05;  
  272. TAB1.'PRECISINTER' = 1E-05;
  273. *TAB1. 'TEMPS0'         = 0.  ;
  274. *TAB1. 'PROCESSEURS' =  'MONO';
  275. TAB1.'TEMPS_CALCULES' = tcalc0 ;
  276. TAB1.'TEMPS_SAUVES'   = tcalc0 ;
  277.  
  278. * Resolution non-lineaire ******
  279.   PASAPAS TAB1;
  280.  
  281. * Mise a Jour de la table TXFEM (a terme devra etre inclus dans PASAPAS)
  282. i = 0 ;
  283. repe BOUX0 (Ncalc0 - 1);
  284.    i = i + 1;
  285.    TXFEM . i = TXFEM . (i - 1);
  286.    Che1X . i = Che1X . (i - 1);
  287. fin BOUX0;
  288.  
  289.  
  290. * calcul de l integrale J *****
  291.  
  292. * * J local
  293. * si (fllong);
  294. *    SUPTAB = TABL ;
  295. *    SUPTAB.'OBJECTIF' = MOT 'J';
  296. *    SUPTAB.'PSI' = psi0;
  297. *    SUPTAB.'PHI' = phi0;
  298. *    SUPTAB.'FRONT_FISSURE' = ptip1 ;
  299. *    SUPTAB.'SOLUTION_PASAPAS' = TAB1;
  300. *    SUPTAB.'COUCHE' = 2;
  301. *    G_THETA SUPTAB;
  302. *    si(graph);
  303. *       q7 = (SUPTAB . 'CHAMP_THET1') ;
  304. *       vq7 = VECT q7 AZUR ;
  305. *       trac vq7 (stot0 et lcrack0);
  306. *    fins;
  307. * fins;
  308.  
  309. * calcul de l integrale J global
  310. JTAB = TABL;
  311. JTAB.'OBJECTIF' = MOT 'J';
  312. JTAB.'PSI' = psi0;
  313. JTAB.'PHI' = phi0;
  314. JTAB.'FRONT_FISSURE' = ptip1 ;
  315. JTAB.'SOLUTION_PASAPAS' = TAB1;
  316. *definition dun champ_theta perso pour calcul de j global
  317. sared = ELEM satot 'APPUYE' 'LARGEMEMENT' (CONT satot);
  318. sared = (satot DIFF sared) coul vert;
  319. monch1 = MANU 'CHPO' sared 'UX' 1. 'UY' 0.;
  320. monvec1 = vect monch1 vert;
  321. si(graph);
  322.   trac monvec1 (stot0 et lcrack0);
  323. *  trac (exco 'UX' monch1) (stot0 et lcrack0);
  324. fins;
  325. JTAB.'CHAMP_THETA' = monch1;
  326. G_THETA JTAB;
  327.  
  328.  
  329. *******************************************************
  330. ***         PROPAGATION élémentaire de fissure      ***
  331. *rem: + tard, prevoir procedure surtout pour le 3d
  332. NBPROPA = 2;
  333. da1 = 1.5 * dx1;
  334.  
  335. repe BPROPA NBPROPA;
  336.  
  337. i = i + 1;
  338.  
  339. *** Maillage de la fissure ***
  340. ptip1 = ptip1 plus (da1 0.);
  341. lcrack1 = (lcrack0 droi ptip1) coul 'ROUG';
  342. TXFEM . i = tabl;
  343. TXFEM . i . 'POINTE1' = ptip1;
  344. TXFEM . i . 'FISSURE' = lcrack1;
  345.  
  346. *** Actualisations ***
  347. * des level set
  348. psi1 phi1 = PSIPHI satot lcrack1 'DEUX' ptip1 ;
  349. si(graph);
  350. *opti isov lign;
  351.   trac (stot0 et lcrack1);
  352.   TRAC phi1 (stot0 et lcrack1) isolv7 ;
  353.   TRAC psi1 (stot0 et lcrack1) isolv7;
  354. *opti isov surf;
  355. finsi;
  356. * du modele, de la rigidite,... (inutile pour le materiau)
  357. Che1X . i = TRIE mod1a psi1 phi1;
  358.  
  359. * des blocages des ddl X-fem non actifs dans 
  360. * les éléments de transition.
  361. Rel1X . i = RELA mod1a;
  362. mod1tot = mod1a et mod1b;
  363. TAB1.'MODELE'              = mod1tot;
  364. TAB1.'BLOCAGES_MECANIQUES' = cl1tot et (Rel1x . i) ;
  365.  
  366. *** resolution non-lineaire PASaPAS ***
  367. tlast1 = extr TAB1.'TEMPS_CALCULES' (dime TAB1.'TEMPS_CALCULES');
  368. tcalc1 = TAB1.'TEMPS_CALCULES' et (PROG (tlast1 + dtval0));
  369. Ncalc1 = dime tcalc1;
  370.  
  371. TAB1.'TEMPS_CALCULES' = tcalc1 ;
  372. TAB1.'TEMPS_SAUVES'   = tcalc1 ;
  373. TAB1.'REEQUILIBRAGE'  = VRAI ;
  374.  
  375. PASAPAS TAB1;
  376.  
  377. *** calcul de l integrale J global (le contour ne change pas) ***
  378. JTAB.'PSI' = psi1;
  379. JTAB.'PHI' = phi1;
  380. JTAB.'FRONT_FISSURE' = ptip1 ;
  381. *JTAB.'CHAMP_THETA' = JTAB.'CHAMP_THET1';
  382. JTAB.'CHAMP_THETA' = monch1;
  383. * on utilise la reprise car le contour ne change pas
  384. G_THETA JTAB;
  385.  
  386. * si (fllong);
  387. * *** calcul de l integrale J local, de K1, de K2 ***
  388. * * a chaque propa il faut une nouvelle table, car le contour change...
  389. * SUPTAB = TABL ;
  390. * SUPTAB.'OBJECTIF' = MOT 'J';
  391. * SUPTAB.'PSI' = psi1;
  392. * SUPTAB.'PHI' = phi1;
  393. * SUPTAB.'FRONT_FISSURE' = ptip1 ;
  394. * SUPTAB.'SOLUTION_PASAPAS' = TAB1;
  395. * SUPTAB.'COUCHE' = 2;
  396. * 
  397. *   G_THETA SUPTAB;
  398. *   q7 = (SUPTAB . 'CHAMP_THET1') ;
  399. * 
  400. *   si(graph);
  401. *     vq7 = VECT q7 AZUR ;
  402. *     trac vq7 (stot0 et lcrack1);
  403. *   fins;
  404. * fins;
  405.  
  406. * i etant l indice de la table de PASAPAS TAB1,
  407. * il faut l'incrementer avant et apres chaque appel a pasapas
  408. i = i + 1;
  409. TXFEM . i = TXFEM . (i - 1);
  410. Che1X . i = Che1X . (i - 1);
  411.  
  412. fin BPROPA;
  413.  
  414.  
  415.  
  416. ******************************************************
  417. ***                POST TRAITEMENT                 ***
  418.  
  419. * courbe J global - t  (le contour ne change pas) ***
  420. si(graph);
  421.    dess JTAB . EVOLUTION_RESULTATS 'TITR' 'Jglobal - t';
  422. finsi;
  423.  
  424. * courbe force ouverture  ****************************
  425. cod1 = prog;
  426. for1 = prog;
  427. ipost = -1;
  428. Npost = i;
  429. repe BPOST1 (Npost + 1);
  430.   ipost = ipost + 1;
  431.   u1i   = TAB1 . 'DEPLACEMENTS' . ipost;
  432.     cod1i = extr u1i 'AY' pa1;
  433.     cod1 = cod1 et (prog cod1i);
  434.   frea1i = TAB1 . 'REACTIONS' . ipost;
  435.     for1i = extr frea1i 'FY' pgoup1;
  436.     for1 = for1 et (prog for1i);
  437. fin BPOST1;
  438. evfcod1 = EVOL 'VERT' 'MANU' 'ouverture' cod1 'force' for1 ;
  439. si(graph);
  440.   dess evfcod1 tdess1 'TITR' 'Force - Ouverture';
  441. finsi;
  442.  
  443. * quelques tracés  ****************************
  444.   u1i = TAB1 . 'DEPLACEMENTS'  . i;
  445.   u1phy = XFEM 'RECO' u1i mod1tot ;
  446.   var1i = TAB1 . 'VARIABLES_INTERNES'  . i;
  447.   sig1i = TAB1 . 'CONTRAINTES'  . i;
  448.   def1i = ELAS mod1tot sig1i mat1tot;
  449. si(graph);
  450.   TRAC var1i mod1tot (DEFO u1phy stot0 20.) 'TITR' 'EPSE';
  451.   TRAC (exco EPYY def1i) mod1tot (DEFO u1phy stot0 20.) 'TITR' 'EPYY';
  452. finsi;
  453.  
  454. *
  455. *opti donn 5 ;
  456. fin;
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales