Télécharger la plaquette Cast3M


fr en

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

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