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. * apres cela on assemble
  123. mod1tot = mod1a ;
  124.  
  125. * on crée le ou les materiau et on les assemble
  126. Ey1 = 5.520E+2 '/' 2.627E-3 ;
  127. nu1 = 0.3 ;
  128. rho1= 7.8e-6 ;
  129.  
  130. Epsprog0 = (PROG
  131. 0.000E+0 2.627E-3 2.967E-3 3.276E-3
  132. 3.585E-3 5.129E-3 8.221E-3 1.132E-2 1.442E-2
  133. 1.754E-2 3.334E-2 5.006E-2 6.914E-2 9.391E-2
  134. 1.314E-1 1.956E-1 3.134E-1 5.336E-1) ;
  135. Sigprog0 = (PROG
  136. 0.000E+0 5.520E+2 5.530E+2 5.540E+2
  137. 5.550E+2 5.600E+2 5.700E+2 5.800E+2 5.900E+2
  138. 6.000E+2 6.500E+2 7.000E+2 7.500E+2 8.000E+2
  139. 8.500E+2 9.000E+2 9.500E+2 1.000E+3);
  140. Tcurve = EVOL 'MANU' 'Epsilon' Epsprog0 'Sigma' Sigprog0 ;
  141.  
  142. mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
  143. 'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
  144.  
  145.  
  146. *******************************************************
  147. *** Matrices + CL
  148. *Rig1tot = RIGI mod1tot mat1tot;
  149.  
  150. * on bloque les depl de corps rigides
  151. pgoup1 = stot0 POIN 'PROC' (0. hgoup);
  152. pgoup2 = stot0 POIN 'PROC' (0. (-1.*hgoup));
  153. pgouptot = pgoup1 et pgoup2;
  154. cl1x = BLOQ 'UX' (pgouptot);
  155. cl1y = BLOQ 'UY' (pgouptot);
  156.  
  157. cl1tot = cl1x et cl1y;
  158.  
  159.  
  160. *******************************************************
  161. *** deplacement impose ***
  162.  
  163. * evolution
  164. du0 = 0.04 ;
  165. ufin0 = 0.44 ;
  166.  
  167. si logdech;
  168. uval0 = PROG 0. 'PAS' du0 ufin0 'PAS' (0. - du0) (0.8*ufin0);
  169. tval0 = (PROG 0. 'PAS' du0 ufin0 'PAS' du0 (1.2*ufin0))* 1.E-2 ;
  170. sinon;
  171. uval0 = PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0;
  172. tval0 = (PROG 0. 'PAS' (2.5*du0) (0.7*ufin0) 'PAS' du0 ufin0)* 1.E-2;
  173. finsi;
  174. Nmax0 = DIME tval0 ;
  175. uev0 = EVOL 'ROUG' 'MANU' 't' tval0 'u' uval0 ;
  176.  
  177. * chpoint
  178. uchp0 = (MANU 'CHPO' pgoup1 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS')
  179. ET (MANU 'CHPO' pgoup2 2 'UX' 0. 'UY' -1. 'NATU' 'DIFFUS');
  180. fchp0 = DEPI cl1tot uchp0;
  181.  
  182. * chargement
  183. uimp1 = CHAR 'DIMP' fchp0 uev0 ;
  184.  
  185. tdess1 = tabl;
  186. tdess1 . 1 = mot 'MARQ PLUS';
  187. si(graph);
  188. dess uev0 tdess1;
  189. trac (vect uchp0 'DEPL' 'VERT') (stot0 et lcrack0);
  190. finsi;
  191.  
  192.  
  193. *******************************************************
  194. *** RESOlution avec PASaPAS
  195.  
  196. * Initialisation calcul
  197.  
  198. TAB1 = TABL;
  199. TAB1.'MODELE' = mod1tot;
  200. TAB1.'CARACTERISTIQUES' = mat1tot;
  201. TAB1.'BLOCAGES_MECANIQUES' = cl1tot;
  202. TAB1.'CHARGEMENT' = uimp1;
  203. TAB1.'PRECISION' = 1E-05;
  204. TAB1.'PRECISINTER' = 1E-05;
  205. *TAB1. 'TEMPS0' = 0. ;
  206. *TAB1. 'PROCESSEURS' = 'MONO';
  207. TAB1.'TEMPS_CALCULES' = tval0 ;
  208. TAB1.'TEMPS_SAUVES' = tval0 ;
  209. * Resolution non-lineaire
  210. PASAPAS TAB1;
  211.  
  212. ****************
  213. * Cas fissure maillée :
  214. sbtot = satot;
  215.  
  216. * definition des models elementaires
  217. mod1b = MODE sbtot mecanique elastique isotrope
  218. plastique_endom rousselier ;
  219. mat1b = MATE mod1b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'TRAC' Tcurve
  220. 'SIG1' 450. 'D' 2. 'F' 0.0014 'FC' 0.35;
  221. * on bloque les depl de corps rigides
  222. dcl0 = d11d et d1d;
  223. cl0 = BLOQ 'UY' dcl0;
  224.  
  225.  
  226. pgoup1b = sbtot POIN 'PROC' (0. hgoup);
  227. cl1xb = BLOQ 'UX' (pgoup1b);
  228. cl1yb = BLOQ 'UY' (pgoup1b);
  229.  
  230. cl1totb = cl1xb et cl1yb ;
  231. cltot = cl1totb et cl0;
  232. *** deplacement impose ***
  233. * chpoint
  234. uchp0b = (MANU 'CHPO' pgoup1b 2 'UX' 0. 'UY' 1. 'NATU' 'DIFFUS');
  235. fchp0b = DEPI cl1totb uchp0b;
  236.  
  237. * chargement
  238. uimp1b = CHAR 'DIMP' fchp0b uev0 ;
  239.  
  240. tdess1 = tabl;
  241. tdess1 . 1 = mot 'MARQ PLUS';
  242. si(graph);
  243. dess uev0 tdess1;
  244. trac (vect uchp0b 'DEPL' 'VERT') (sbtot et lcrack0);
  245. finsi;
  246.  
  247.  
  248. *******************************************************
  249. * Initialisation calcul
  250.  
  251. TAB2 = TABL;
  252. TAB2.'MODELE' = mod1b;
  253. TAB2.'CARACTERISTIQUES' = mat1b;
  254. TAB2.'BLOCAGES_MECANIQUES' = cltot;
  255. TAB2.'CHARGEMENT' = uimp1b;
  256. TAB2.'PRECISION' = 1E-05;
  257. TAB2.'PRECISINTER' = 1E-05;
  258. *TAB2. 'PROCESSEURS' = 'MONO';
  259. TAB2.'TEMPS_CALCULES' = tval0 ;
  260. TAB2.'TEMPS_SAUVES' = tval0 ;
  261. * Resolution non-lineaire
  262. PASAPAS TAB2;
  263.  
  264.  
  265.  
  266. ******************************************************
  267. * post traitement de pasapas
  268.  
  269. * courbes force ouverture ****************************
  270. cod1 = prog ; for1 = prog;
  271. cod2 = prog ; for2 = prog;
  272. i = -1;
  273. repe BPOST1 (Nmax0);
  274. i = i + 1;
  275. u1i = TAB1 . 'DEPLACEMENTS' . i;
  276. cod1i = extr u1i 'UY' p4;
  277. cod1 = cod1 et (prog cod1i);
  278. frea1i = TAB1 . 'REACTIONS' . i;
  279. for1i = extr frea1i 'FY' pgoup1;
  280. for1 = for1 et (prog for1i);
  281.  
  282. u2i = TAB2 . 'DEPLACEMENTS' . i;
  283. cod2i = extr u2i 'UY' p4;
  284. cod2 = cod2 et (prog cod2i);
  285. frea2i = TAB2 . 'REACTIONS' . i;
  286. for2i = extr frea2i 'FY' pgoup1b;
  287. for2 = for2 et (prog for2i);
  288. fin BPOST1;
  289. evfcod1 = EVOL 'BLEU' 'MANU' 'ouverture' cod1 'force' for1 ;
  290. evfcod2 = EVOL 'ROUG' 'MANU' 'ouverture' cod2 'force' for2 ;
  291. si(graph);
  292. dess (evfcod1 et evfcod2);
  293. finsi;
  294.  
  295. * tracés a la fin de la mise en charge **************
  296. sig1i = TAB1 . 'CONTRAINTES' . i;
  297. var1i = TAB1 . 'VARIABLES_INTERNES' . i;
  298. svm1 = VMIS sig1i mod1tot;
  299. sig2i = TAB2 . 'CONTRAINTES' . i;
  300. var2i = TAB2 . 'VARIABLES_INTERNES' . i;
  301. svm2 = VMIS sig2i mod1b;
  302. si(graph);
  303. TRAC svm1 mod1tot (DEFO u1i (stot0 et lcrack0) 3.);
  304. TRAC svm2 mod1b (DEFO u2i (sbtot et lcrack0) 3.);
  305. TRAC var1i mod1tot (DEFO u1i stot0 3.);
  306. TRAC var2i mod1b (DEFO u2i sbtot 3.);
  307. finsi;
  308.  
  309. * Reconstruction du déplacement réel pour avoir la déformée :
  310. chp2UX = XFEM 'RECO' tab1.deplacements. i mod1tot ;
  311. si graph;
  312. stot0d = (defo (coul stot0 roug) chp2UX) ;
  313. sbtotd = defo (coul sbtot jaun) tab2.deplacements. i;
  314. trac stot0d;
  315. trac (stot0d et (defo stot0 tab1.deplacements. i));
  316. finsi;
  317.  
  318. *******************************************************
  319. * propagation élémentaire de fissure
  320. * ... a venir
  321. *==============================================================*
  322. *
  323. err6 = (ABS (((maxi for2) - (maxi for1))/(maxi for2)));
  324. si logmess;
  325. mess 'erreur sur l effort : ' err6;
  326. finsi;
  327.  
  328.  
  329. SI (err6 >EG 1.);
  330. ERRE 5;
  331. FINSI;
  332.  
  333.  
  334.  
  335. *opti donn 5 ;
  336. fin;
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  

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