Télécharger xfem02.dgibi

Retour à la liste

Numérotation des lignes :

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

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