Télécharger xfem03.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : xfem03.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ********************************************************************
  5. *
  6. * xfem03.dgibi
  7. *
  8. * calcul avec elements xfem (XQ4R)
  9. * d'une plaque en traction modele de Rousselier avec fissure droite
  10. *
  11. * création : as, le 24.11.2009 : Résultats moyens à cause des fonctions
  12. * de pointes de fissure adaptées à la rupture fragile sans prise en compte
  13. * de la plasticité étendue !!!
  14. *
  15. ********************************************************************
  16. *opti echo 1 ;
  17. graph = faux;
  18. logmess = faux;
  19. * si on veut tester la décharge:
  20. logdech = faux;
  21.  
  22. *******************************************************
  23. *** Options de calcul : 2D
  24. opti dime 2 elem qua4 mode plan defo;
  25.  
  26.  
  27. *******************************************************
  28. *** Maillage : Plaque rectangulaire saine
  29.  
  30. * dimensions
  31. * bord droit, gauche, hauteur
  32. b= 25. ; bg = -0. * b ; hb1= 0.5 * b ;
  33. hgoup= 0.35 * b;
  34.  
  35. * finesse, pour la comparaison FEM/XFEM, un nombre d'éléments
  36. * pair est nécessaire
  37. nx1 = 6; dx1 = (b-bg) / (flot nx1) ;
  38. ha1 = (flot nx1) * dx1;
  39.  
  40. * demi-zone centrale
  41. DENS (dx1) ;
  42. p1 = bg 0. ; p2 = b 0. ;
  43. p3 = b hb1 ; p4 = bg hb1 ;
  44.  
  45. p11 = (bg + ((b-bg)/2.5)) 0. ; p33 = (b - ((b-bg)/2.5)) (hb1/5.);
  46. p44 = p11 plus (0. (coor p33 2)); p22 = p33 moin (0. (coor p33 2));
  47.  
  48. pa1 = p1 plus (((b-bg)/2.) 0.);
  49.  
  50. d1g = droi (nx1/2) p1 p11; d11g = droi (nx1/2) p11 pa1;
  51. d11d = droi (nx1/2) pa1 p22; d1d = droi (nx1/2) p22 p2;
  52. d11 = d11d et d11g; d22 = droi (nx1/2) p22 p33;
  53. d33 = droi nx1 p33 p44; d44 = droi (nx1/2) p44 p11;
  54.  
  55. d2 = droi (nx1/2) p2 p3; d4 = droi (nx1/2) p4 p1;
  56. d3 = droi nx1 p3 p4;
  57.  
  58. d333 = droi (nx1/2) p33 p3;
  59.  
  60. saco = dall d11 d22 d33 d44;
  61.  
  62. sah = regl d33 (nx1/2) d3; sag = regl d44 (nx1/2) d4;
  63. sad = dall d22 d1d d2 d333;
  64.  
  65. satot = saco et sah et sad et sag;
  66.  
  67. stot0 = satot et (syme droi p1 p2 satot);
  68. ELIM stot0 (1.E-5*dx1);
  69.  
  70. si(logmess);
  71. mess 'maillage central ext total' ;
  72. mess 'nbno= ' (nbno satot) (nbno stot0);
  73. mess 'nbel= ' (nbel satot) (nbel stot0);
  74. finsi;
  75.  
  76.  
  77. *******************************************************
  78. *** Maillage de la fissure
  79.  
  80. * longueur de la fissure
  81. a0 = (b - bg) /2. ;
  82.  
  83. * finesse du maillage de la fissure
  84. na0 = 10;
  85. da0 = a0 / (flot na0) ;
  86.  
  87. * maillage
  88. DENS (da0) ;
  89. ptip0 = (bg) (0.);
  90. ptip1 = (a0 + bg) (0.);
  91. lcrack0 = (ptip0 droi na0 ptip1) coul 'ROUG';
  92.  
  93. * table (+ pratique pour gérer la propa)
  94. TXFEM = tabl;
  95. TXFEM . 0 = tabl;
  96. TXFEM . 0 . 'POINTE1' = ptip1;
  97. TXFEM . 0 . 'FISSURE' = lcrack0;
  98.  
  99. * creation des level set
  100. psi0 phi0 = PSIPHI stot0 lcrack0 'DEUX' ptip1;
  101. isolv7 = ( prog (-1.*a0) (-5.*dx1) PAS dx1 (-1.*dx1)
  102. PAS (0.2*dx1) dx1 PAS dx1 (5.*dx1) a0 );
  103. si(graph);
  104. *opti isov lign;
  105. TRAC phi0 (stot0 et lcrack0) isolv7 ;
  106. TRAC psi0 (stot0 et lcrack0) isolv7;
  107. *opti isov surf;
  108. finsi;
  109.  
  110.  
  111. *******************************************************
  112. *** Modeles et materiaux
  113.  
  114. * definition des models elementaires
  115. mod1a = MODE stot0 mecanique elastique isotrope
  116. plastique_endom rousselier xq4R;
  117.  
  118. * enrichissement
  119. Che1X = tabl;
  120. Che1X . 0 = TRIE mod1a psi0 phi0;
  121.  
  122. * constructionsion des blocages des ddl X-fem non actifs dans
  123. * les éléments de transition.
  124. Rel1X = tabl;
  125. Rel1x . 0 = RELA mod1a;
  126.  
  127. * apres cela on assemble
  128. mod1tot = mod1a ;
  129.  
  130. * on crée le ou les materiau et on les assemble
  131. Ey1 = 5.520E+2 '/' 2.627E-3 ;
  132. nu1 = 0.3 ;
  133. rho1= 7.8e-6 ;
  134.  
  135. Epsprog0 = (PROG
  136. 0.000E+0 2.627E-3 2.967E-3 3.276E-3
  137. 3.585E-3 5.129E-3 8.221E-3 1.132E-2 1.442E-2
  138. 1.754E-2 3.334E-2 5.006E-2 6.914E-2 9.391E-2
  139. 1.314E-1 1.956E-1 3.134E-1 5.336E-1) ;
  140. Sigprog0 = (PROG
  141. 0.000E+0 5.520E+2 5.530E+2 5.540E+2
  142. 5.550E+2 5.600E+2 5.700E+2 5.800E+2 5.900E+2
  143. 6.000E+2 6.500E+2 7.000E+2 7.500E+2 8.000E+2
  144. 8.500E+2 9.000E+2 9.500E+2 1.000E+3);
  145. Tcurve = EVOL 'MANU' 'Epsilon' Epsprog0 'Sigma' Sigprog0 ;
  146.  
  147. mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
  148. 'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
  149.  
  150.  
  151. *******************************************************
  152. *** Matrices + CL
  153. *Rig1tot = RIGI mod1tot mat1tot;
  154.  
  155. * on bloque les depl de corps rigides
  156. pgoup1 = stot0 POIN 'PROC' (0. hgoup);
  157. pgoup2 = stot0 POIN 'PROC' (0. (-1.*hgoup));
  158. pgouptot = pgoup1 et pgoup2;
  159. cl1x = BLOQ 'UX' (pgouptot);
  160. cl1y = BLOQ 'UY' (pgouptot);
  161.  
  162. cl1tot = cl1x et cl1y;
  163.  
  164.  
  165. *******************************************************
  166. *** deplacement impose ***
  167.  
  168. * evolution
  169. du0 = 0.04 ;
  170. ufin0 = 0.44 ;
  171.  
  172. si logdech;
  173. uval0 = PROG 0. 'PAS' du0 ufin0 'PAS' (0. - du0) (0.8*ufin0);
  174. tval0 = (PROG 0. 'PAS' du0 ufin0 'PAS' du0 (1.2*ufin0))* 1.E-2 ;
  175. sinon;
  176. uval0 = PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0;
  177. tval0 = (PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0)* 1.E-2;
  178. finsi;
  179. Nmax0 = DIME tval0 ;
  180. uev0 = EVOL 'ROUG' 'MANU' 't' tval0 'u' uval0 ;
  181.  
  182. * chpoint
  183. uchp0 = (MANU 'CHPO' pgoup1 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS')
  184. ET (MANU 'CHPO' pgoup2 2 'UX' 0. 'UY' -1. 'NATU' 'DIFFUS');
  185. fchp0 = DEPI cl1tot uchp0;
  186.  
  187. * chargement
  188. uimp1 = CHAR 'DIMP' fchp0 uev0 ;
  189.  
  190. tdess1 = tabl;
  191. tdess1 . 1 = mot 'MARQ PLUS';
  192. si(graph);
  193. dess uev0 tdess1;
  194. trac (vect uchp0 'DEPL' 'VERT') (stot0 et lcrack0);
  195. finsi;
  196.  
  197.  
  198. *******************************************************
  199. *** RESOlution avec PASaPAS
  200.  
  201. * Initialisation calcul
  202.  
  203. TAB1 = TABL;
  204. TAB1.'MODELE' = mod1tot;
  205. TAB1.'CARACTERISTIQUES' = mat1tot;
  206. TAB1.'BLOCAGES_MECANIQUES' = cl1tot et (Rel1X.0);
  207. TAB1.'CHARGEMENT' = uimp1;
  208. TAB1.'PRECISION' = 1E-05;
  209. TAB1.'PRECISINTER' = 1E-05;
  210. *TAB1. 'TEMPS0' = 0. ;
  211. *TAB1. 'PROCESSEURS' = 'MONO';
  212. TAB1.'TEMPS_CALCULES' = tval0 ;
  213. TAB1.'TEMPS_SAUVES' = tval0 ;
  214. * Resolution non-lineaire
  215. PASAPAS TAB1;
  216.  
  217. ****************
  218. * Cas fissure maillée :
  219. sbtot = satot;
  220.  
  221. * definition des models elementaires
  222. mod1b = MODE sbtot mecanique elastique isotrope
  223. plastique_endom rousselier ;
  224. mat1b = MATE mod1b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
  225. 'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
  226. * on bloque les depl de corps rigides
  227. dcl0 = d11d et d1d;
  228. cl0 = BLOQ 'UY' dcl0;
  229.  
  230.  
  231. pgoup1b = sbtot POIN 'PROC' (0. hgoup);
  232. cl1xb = BLOQ 'UX' (pgoup1b);
  233. cl1yb = BLOQ 'UY' (pgoup1b);
  234.  
  235. cl1totb = cl1xb et cl1yb ;
  236. cltot = cl1totb et cl0;
  237. *** deplacement impose ***
  238. * chpoint
  239. uchp0b = (MANU 'CHPO' pgoup1b 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS');
  240. fchp0b = DEPI cl1totb uchp0b;
  241.  
  242. * chargement
  243. uimp1b = CHAR 'DIMP' fchp0b uev0 ;
  244.  
  245. tdess1 = tabl;
  246. tdess1 . 1 = mot 'MARQ PLUS';
  247. si(graph);
  248. dess uev0 tdess1;
  249. trac (vect uchp0b 'DEPL' 'VERT') (sbtot et lcrack0);
  250. finsi;
  251.  
  252.  
  253. *******************************************************
  254. * Initialisation calcul
  255.  
  256. TAB2 = TABL;
  257. TAB2.'MODELE' = mod1b;
  258. TAB2.'CARACTERISTIQUES' = mat1b;
  259. TAB2.'BLOCAGES_MECANIQUES' = cltot;
  260. TAB2.'CHARGEMENT' = uimp1b;
  261. TAB2.'PRECISION' = 1E-05;
  262. TAB2.'PRECISINTER' = 1E-05;
  263. *TAB2. 'PROCESSEURS' = 'MONO';
  264. TAB2.'TEMPS_CALCULES' = tval0 ;
  265. TAB2.'TEMPS_SAUVES' = tval0 ;
  266. * Resolution non-lineaire
  267. PASAPAS TAB2;
  268.  
  269.  
  270.  
  271. ******************************************************
  272. * post traitement de pasapas
  273.  
  274. * courbes force ouverture ****************************
  275. cod1 = prog ; for1 = prog;
  276. cod2 = prog ; for2 = prog;
  277. i = -1;
  278. repe BPOST1 (Nmax0);
  279. i = i + 1;
  280. u1i = TAB1 . 'DEPLACEMENTS' . i;
  281. cod1i = extr u1i 'UY' p4;
  282. cod1 = cod1 et (prog cod1i);
  283. frea1i = TAB1 . 'REACTIONS' . i;
  284. for1i = extr frea1i 'FY' pgoup1;
  285. for1 = for1 et (prog for1i);
  286.  
  287. u2i = TAB2 . 'DEPLACEMENTS' . i;
  288. cod2i = extr u2i 'UY' p4;
  289. cod2 = cod2 et (prog cod2i);
  290. frea2i = TAB2 . 'REACTIONS' . i;
  291. for2i = extr frea2i 'FY' pgoup1b;
  292. for2 = for2 et (prog for2i);
  293. fin BPOST1;
  294. evfcod1 = EVOL 'BLEU' 'MANU' 'ouverture' cod1 'force' for1 ;
  295. evfcod2 = EVOL 'ROUG' 'MANU' 'ouverture' cod2 'force' for2 ;
  296. si(graph);
  297. dess (evfcod1 et evfcod2);
  298. finsi;
  299.  
  300. * tracés a la fin de la mise en charge **************
  301. sig1i = TAB1 . 'CONTRAINTES' . i;
  302. var1i = TAB1 . 'VARIABLES_INTERNES' . i;
  303. svm1 = VMIS sig1i mod1tot;
  304. sig2i = TAB2 . 'CONTRAINTES' . i;
  305. var2i = TAB2 . 'VARIABLES_INTERNES' . i;
  306. svm2 = VMIS sig2i mod1b;
  307. si(graph);
  308. TRAC svm1 mod1tot (DEFO u1i (stot0 et lcrack0) 3.);
  309. TRAC svm2 mod1b (DEFO u2i (sbtot et lcrack0) 3.);
  310. TRAC var1i mod1tot (DEFO u1i stot0 3.);
  311. TRAC var2i mod1b (DEFO u2i sbtot 3.);
  312. finsi;
  313.  
  314. * Reconstruction du déplacement réel pour avoir la déformée :
  315. chp2UX = XFEM 'RECO' tab1.deplacements. i mod1tot ;
  316. si graph;
  317. stot0d = (defo (coul stot0 roug) chp2UX) ;
  318. sbtotd = defo (coul sbtot jaun) tab2.deplacements. i;
  319. trac stot0d;
  320. trac (stot0d et (defo stot0 tab1.deplacements. i));
  321. finsi;
  322.  
  323. *******************************************************
  324. * propagation élémentaire de fissure
  325. * ... a venir
  326. *==============================================================*
  327. *
  328. err6 = (ABS (((maxi for2) - (maxi for1))/(maxi for2)));
  329. si logmess;
  330. mess 'erreur sur l effort : ' err6;
  331. finsi;
  332.  
  333.  
  334. SI (err6 >EG 1.);
  335. ERRE 5;
  336. FINSI;
  337.  
  338.  
  339.  
  340. *opti donn 5 ;
  341. fin;
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  

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