Télécharger la plaquette Cast3M


fr en

Télécharger xfem_gd.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : xfem_gd.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. opti DIME 2 ELEM QUA4 mode plan defo;
  5.  
  6. ******************************************************************** 
  7. *              xfemcapipica.dgibi
  8. *
  9. * Test de l'opérateur de passage des contraintes(déformations) 
  10. *    PK2 aux contraintes de Cauchy pour la XFEM
  11. * d'une plaque elastique en traction avec fissure droite 
  12. * 
  13. * création : as, le 14.01.2010 
  14. ********************************************************************
  15. opti echo 0;
  16. * Options :
  17. *########
  18. * pour afficher les messages et les graphes + choix de la sortie en .ps:
  19. lmess = faux ;              lgraph = faux ;
  20. si lgraph ; lps=vrai; sinon; lps= faux; finsi;
  21. si lps;   opti trac psc; opti ftra 'xfem_gd.ps'; finsi;
  22.  
  23. * Maillage :
  24. * ########
  25. dx1 = 0.06 ;          
  26. blim = 1. + 5 * dx1 ;        nbr1  = enti (blim / dx1);
  27. * densité sur la hauteur de la bande :
  28. nbh1 = 3 ;  hh1 = ((flot nbh1)*dx1);
  29.  
  30. p1 =  'POIN' 0.  hh1  ;               p2 = 'POIN' blim hh1 ; 
  31. pm3 =  'POIN' blim  (0.) ;        pm4 = 'POIN' 0. (0.) ;    
  32. p3 = 'POIN' blim  (-1.*hh1) ;     p4 = 'POIN' 0.  (-1.*hh1) ; 
  33.  
  34. lh =  droi nbr1 p1 p2 ; 
  35. ldh =  droi p2 pm3 DINI dx1 DFIN dx1;
  36. lgh =  droi p1 pm4 DINI dx1 DFIN dx1;
  37. lm = droi nbr1 pm4 pm3;
  38. ldb =  droi pm3 p3 DINI dx1 DFIN dx1;
  39. lgb =  droi pm4 p4 DINI dx1 DFIN dx1;
  40. lb = droi nbr1 p4 p3;
  41.  
  42. s1h = 'DALLER' lm ldh lh lgh ;
  43. s1b = 'DALLER' lb ldb lm lgb ;
  44.  
  45. * fissure initiale :
  46. * #################
  47. a0 = 3.*dx1;
  48. pfis0 = -0.01 0.;    pfisi = a0 0.; 
  49. fis0 = coul (droi 1 pfis0 pfisi) roug;
  50.  
  51. * Modèle et materiau :
  52. * ####################
  53. Ey1 = 2.e5 ; nu1=0.3; rho1 = 7.e-6 ;
  54.  
  55. * CHARGEMENT : déplacement imposé :
  56. * #################################
  57. tf = 1.;   NT = ENTI (3.*tf);  
  58. dt0 = tf / (flot NT) ;        l_t = PROG 0. pas dt0 tf; 
  59. ufin = tf ;  l_uimp = PROG 0. pas (ufin / (flot NT)) ufin;
  60. evou = evol 'MANU' 'temps' l_t 'U_imp' l_uimp ;
  61. coeff0 = 0.2;
  62.  
  63.  
  64. **###################
  65. **   Calcul MASSIF  **
  66. **###################
  67. s1q = s1h ;
  68. elim s1q 0.0001;
  69. * Modèle et materiau :
  70. * ####################
  71. Mod1q = s1q MODE mecanique elastique isotrope ;
  72.  
  73. mat1q = MATE Mod1q 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  74.  
  75.  
  76. *   Conditions aux limites, Déplacements imposés :
  77. * ################################################
  78. ptrach = poin s1q proc (p1 plus ((2.*dx1) (-2.*dx1) ));
  79. cl1hx = 'BLOQUE' ptrach 'UX' ;          cl1hy = 'BLOQUE' ptrach 'UY' ;
  80. psyms = elem lm comp (poin s1q proc pfisi) pm3;
  81. cl1sym = 'BLOQUE' psyms 'UY' ;
  82. si lgraph ; trac (s1q et (coul psyms roug));  finsi;
  83.  
  84. cl0 =  (cl1hx 'ET' cl1hy 'ET' cl1sym) ;  
  85.  
  86. dep1  = DEPI cl1hy coeff0;
  87. uimp0 = CHAR 'DIMP' dep1 evou ;
  88.  
  89. *        Resolution iterative du probleme
  90. * ################################################
  91. TAB0 = TABL;   TAB0.TEMPS0 = 0.  ;
  92. TAB0.'MODELE'              = Mod1q  ;
  93. TAB0.'CARACTERISTIQUES'    = mat1q  ;
  94. TAB0.'CHARGEMENT'          = uimp0;
  95. TAB0.'BLOCAGES_MECANIQUES' = cl0   ;
  96. TAB0.'TEMPS_CALCULES' = L_t ;     TAB0.'TEMPS_SAUVES' = L_t  ;
  97. TAB0.'CONVERGENCE_FORCEE' = faux;
  98. TAB0.'DELTAITER' = 190;
  99. TAB0.'MAXITERATION' = 199;
  100. TAB0.'HYPOTHESE_DEFORMATIONS' = 'LINEAIRE';
  101.  
  102. PASAPAS tab0;
  103.  
  104. *opti donn 5;
  105.  
  106. **##################
  107. **   Calcul XFEM  **
  108. **##################
  109. s1 = s1h et s1b;
  110. elim s1 0.0001;
  111.  
  112. si lgraph; trac (s1 et fis0); finsi;
  113.  
  114. * Modèle et materiau :
  115. * ####################
  116. Mod1a = s1 MODE mecanique elastique isotrope xq4R;
  117. Mod1 = Mod1a ;
  118.  
  119. mat1 = MATE Mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  120.  
  121.  
  122. ** Table XFEM : fisssure initiale, critere d'initiation, modele zone cohesive
  123. ** #################################################################
  124. tab_xfem = tabl;
  125. TAB_xfem.fissure = fis0 ;
  126. * creation des level set 
  127. psi1 phi1 = PSIPHI s1 fis0 'DEUX'  pfisi;
  128. * Appel a TRIELE :
  129. cheX1 = trie mod1 psi1 phi1 'SAUT';
  130.  
  131.  
  132. *   Conditions aux limites, Déplacements imposés :
  133. * ################################################
  134. ptrach = poin s1 proc (p1 plus ((2.*dx1) (-2.*dx1) ));
  135. ptracb = poin s1 proc (p4 plus ((2.*dx1) (2.*dx1) ));
  136. cl1hx = 'BLOQUE' ptrach 'UX' ;          cl1bx = 'BLOQUE' ptracb 'UX' ;
  137. cl1hy = 'BLOQUE' ptrach 'UY' ;          cl1by = 'BLOQUE' ptracb 'UY' ;
  138.  
  139. cl1 =  (cl1hx 'ET' cl1bx 'ET'cl1hy 'ET' cl1by) ;  
  140.  
  141. * CHARGEMENT : déplacement imposé :
  142. dep1  = (DEPI cl1hy coeff0) et (DEPI cl1by (-1.*coeff0))   ;
  143. uimp0 = CHAR 'DIMP' dep1 evou ;
  144.  
  145.  
  146. *        Resolution iterative du probleme
  147. * ################################################
  148. TAB1 = TABL;   TAB1.TEMPS0 = 0.  ;
  149. TAB1.'MODELE'              = Mod1  ;
  150. TAB1.'CARACTERISTIQUES'    = mat1  ;
  151. TAB1.'CHARGEMENT'          = uimp0;
  152. TAB1.'BLOCAGES_MECANIQUES' = cl1   ;
  153. TAB1.'TEMPS_CALCULES' = L_t ;     TAB1.'TEMPS_SAUVES' = L_t  ;
  154. TAB1.'CONVERGENCE_FORCEE' = faux;
  155. TAB1.'DELTAITER' = 190;
  156. TAB1.'MAXITERATION' = 199;
  157. TAB1.'HYPOTHESE_DEFORMATIONS' = 'LINEAIRE';
  158.  
  159. *opti donn 5 ;
  160.  
  161. pasapas tab1;
  162.  
  163.  
  164. **######################
  165. **   Post-traitement  **
  166. *########################
  167. naff = mini (lect ((dime tab1.deplacements)  - 1) 
  168.                   ((dime tab0.deplacements)  - 1));
  169.  
  170. si lgraph; 
  171.     trac ((defo s1q tab1.deplacements.naff 1.)
  172.        et (defo (coul s1h roug) tab0.deplacements.naff 1.));
  173. finsi;
  174. conf0=form;
  175.  
  176.  
  177. **#####################
  178. * TEST de PICA et CAPI
  179. **######################
  180.  
  181. CAU1 = pica tab1.contraintes.naff tab1.deplacements.naff 
  182.             tab1.deplacements. 0 tab1.modele;
  183. PIC1 = capi CAU1 tab1.deplacements.naff tab1.deplacements. 0
  184.                                  tab1.modele;
  185. CAU0 = pica tab0.contraintes.naff tab0.deplacements.naff 
  186.             tab0.modele;
  187. PIC0 = capi CAU0 tab0.deplacements.naff tab0.modele;
  188.  
  189. si lmess ;
  190.     mess '  ';
  191.     mess 'Test de CAPI et PICA en grandes defs ';
  192.     mess '  ';
  193.     mess 'Configuration de départ = configuration initiale'; 
  194.     mess '----------------------';
  195.     mess '  ';
  196.     mess 'Elements enrichis :';
  197. finsi;
  198. el1 = 3;  pg1 = 61;
  199. sig01 = extr tab1.contraintes.naff 'SMYY' 1 el1 pg1;
  200. cau11 = extr CAU1 'SMYY' 1 el1 pg1;
  201. pic11 = extr PIC1 'SMYY' 1 el1 pg1;
  202. si lmess; mess 'XFEM:' 'ini' sig01 'pica1' cau11 'inicapi' pic11; fins;
  203.  
  204. pg0 = 3;
  205. sig00 = extr tab0.contraintes.naff 'SMYY' 1 el1 pg0;
  206. cau01 = extr CAU0 'SMYY' 1 el1 pg0;
  207. pic01 = extr PIC0 'SMYY' 1 el1 pg0;
  208. si lmess; mess 'QUA4:' 'ini' sig00 'pica1' cau01 'inicapi' pic01; fins;
  209.  
  210. err6001 = abs ((sig01-sig00)/sig00);
  211. err6002 = abs ((cau11-cau01)/cau01);
  212. err6003 = abs ((pic11-pic01)/pic01);
  213. si lmess; mess '   '; fins;
  214.  
  215. el3 = 2;  
  216. sig013 = extr tab1.contraintes.naff 'SMYY' 1 el3 pg1;
  217. cau13 = extr CAU1 'SMYY' 1 el3 pg1;
  218. pic13 = extr PIC1 'SMYY' 1 el3 pg1;
  219. si lmess; mess 'XFEM:' 'ini' sig013 'pica1' cau13 'inicapi' pic13;fins;
  220. sig003 = extr tab0.contraintes.naff 'SMYY' 1 el3 pg0;
  221. cau03 = extr CAU0 'SMYY' 1 el3 pg0;
  222. pic03 = extr PIC0 'SMYY' 1 el3 pg0;
  223. si lmess; mess 'QUA4:' 'ini' sig003 'pica1' cau03 'inicapi' pic03;fins;
  224. err5001 = abs ((sig013-sig003)/sig003);
  225. err5002 = abs ((cau13-cau03)/cau03);
  226. err5003 = abs ((pic13-pic03)/pic03);
  227.  
  228. si lmess; mess '   '; fins;
  229.  
  230. si lmess; mess  'Element non enrichi:'; fins;
  231. el2 = 11; 
  232. sig012 = extr tab1.contraintes.naff 'SMYY' 1 el2 pg1;
  233. cau12 = extr CAU1 'SMYY' 1 el2 pg1;
  234. pic12 = extr PIC1 'SMYY' 1 el2 pg1;
  235. si lmess; mess 'XFEM:' 'ini' sig012 'pica1' cau12 'inicapi' pic12;fins;
  236. sig002 = extr tab0.contraintes.naff 'SMYY' 1 el2 pg0;
  237. cau02 = extr CAU0 'SMYY' 1 el2 pg0;
  238. pic02 = extr PIC0 'SMYY' 1 el2 pg0;
  239. si lmess; mess 'QUA4:' 'ini' sig002 'pica1' cau02 'inicapi' pic02;fins;
  240. err26001 = abs ((sig012-sig002)/sig002);
  241. err26002 = abs ((cau12-cau02)/cau02);
  242. err26003 = abs ((pic12-pic02)/pic02);
  243.  
  244. errmax1 = maxi (prog err6001 err6002 err6003 err5001 err5002 
  245.            err5003 err26001 err26002 err26003 );
  246.  
  247. * Test de CAPI et PICA en grandes defs
  248. *
  249. * Configuration de départ = configuration initiale'; 
  250. * -------------------------------------------------
  251.  
  252. * XFEM :
  253. uvr1 = XFEM 'RECO' tab1.deplacements. naff mod1 ; 
  254. FORM uvr1;
  255.  
  256. PIC111 = capi CAU1
  257.   (tab1.deplacements. 0 - tab1.deplacements.naff) 
  258. tab1.deplacements. naff  tab1.modele;
  259. CAU111 = pica pic111 (tab1.deplacements. 0 - tab1.deplacements.naff)
  260.             tab1.deplacements. naff tab1.modele;
  261.  
  262. * Standard :
  263. form conf0;
  264. form tab0.deplacements.naff;
  265.  
  266. PIC011 = capi CAU0 (tab0.deplacements. 0 - tab0.deplacements.naff)
  267.                 tab0.modele;
  268. CAU011 = pica pic011 (tab0.deplacements. 0 - tab0.deplacements.naff) 
  269.             tab0.modele;
  270.  
  271. si lmess; 
  272.     mess '  ';
  273.     mess 'Configuration de départ = configuration déformée'; 
  274.     mess '----------------------';
  275.     mess '  ';
  276.     mess 'Elements enrichis :';
  277. fins;
  278. el1 = 3;  pg1 = 61;
  279. sig01 = extr CAU1 'SMYY' 1 el1 pg1;
  280. cau11 = extr CAU111 'SMYY' 1 el1 pg1;
  281. pic11 = extr PIC111 'SMYY' 1 el1 pg1;
  282. si lmess; mess 'XFEM:'  'ini' sig01 'pica1' cau11 'inicapi' pic11;fins;
  283.  
  284. pg0 = 3;
  285. sig00 = extr CAU0 'SMYY' 1 el1 pg0;
  286. cau01 = extr CAU011 'SMYY' 1 el1 pg0;
  287. pic01 = extr PIC011 'SMYY' 1 el1 pg0;
  288. si lmess; mess 'QUA4:'  'ini' sig00 'pica1' cau01 'inicapi' pic01;fins;
  289. err6001 = abs ((sig01-sig00)/sig00);
  290. err6002 = abs ((cau11-cau01)/cau01);
  291. err6003 = abs ((pic11-pic01)/pic01);
  292.  
  293. si lmess; mess '   ';fins;
  294.  
  295. el3 = 2;  
  296. sig013 = extr CAU1 'SMYY' 1 el3 pg1;
  297. cau13 = extr CAU111 'SMYY' 1 el3 pg1;
  298. pic13 = extr PIC111 'SMYY' 1 el3 pg1;
  299. si lmess; mess 'XFEM:' 'ini' sig013 'pica1' cau13 'inicapi' pic13;fins;
  300. sig003 = extr CAU0 'SMYY' 1 el3 pg0;
  301. cau03 = extr CAU011 'SMYY' 1 el3 pg0;
  302. pic03 = extr PIC011 'SMYY' 1 el3 pg0;
  303. si lmess; mess 'QUA4:' 'ini' sig003 'pica1' cau03 'inicapi' pic03;fins;
  304. err5001 = abs ((sig013-sig003)/sig003);
  305. err5002 = abs ((cau13-cau03)/cau03);
  306. err5003 = abs ((pic13-pic03)/pic03);
  307.  
  308. si lmess; mess '   '; fins;
  309.  
  310. si lmess; mess  'Element non enrichi:'; fins;
  311. el2 = 11; 
  312. sig012 = extr CAU1 'SMYY' 1 el2 pg1;
  313. cau12 = extr CAU111 'SMYY' 1 el2 pg1;
  314. pic12 = extr PIC111 'SMYY' 1 el2 pg1;
  315. si lmess; mess 'XFEM:' 'ini' sig012 'pica1' cau12 'inicapi' pic12;fins;
  316. sig002 = extr CAU0 'SMYY' 1 el2 pg0;
  317. cau02 = extr CAU011 'SMYY' 1 el2 pg0;
  318. pic02 = extr PIC011 'SMYY' 1 el2 pg0;
  319. si lmess; mess 'QUA4:' 'ini' sig002 'pica1' cau02 'inicapi' pic02;fins;
  320.  
  321. err26001 = abs ((sig012-sig002)/sig002);
  322. err26002 = abs ((cau12-cau02)/cau02);
  323. err26003 = abs ((pic12-pic02)/pic02);
  324.  
  325. errmax2 = maxi (prog err6001 err6002 err6003 err5001 err5002 
  326.            err5003 err26001 err26002 err26003 );
  327.  
  328. errmax = maxi (prog errmax1 errmax2);
  329.  
  330. si lmess; 
  331.     mess 'Erreur entre calcul standard et XFEM:' errmax;
  332. fins;
  333. * Erreur initiale entre calcul standard et XFEM: 0.10829
  334. * Attention, on fait des erreurs assez importantes car le nombre de 
  335. * pts de Gauss n'est pas le même dans les deux cas :
  336. * -> différences sur le calcul des contraintes initiales
  337. * -> le pt de gauss xq4r est le plus proche possible du pt QUA4
  338. *    mais reste different, 'où une erreur importante pour les éléments
  339. *    en pointe de fissure. 
  340. * pr etre sur, comparer avec éléments xq4r à 4pts de Gauss ->erreur~e-15
  341.  
  342. SI (errmax >EG 20.);
  343.    ERRE 5;
  344. FINSI;
  345.  
  346. fin;
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  

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