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. * ajout de option epsilon lineaire pour la precision des test!
  17. OPTION epsilon lineaire;
  18. * Options :
  19. *########
  20. * pour afficher les messages et les graphes + choix de la sortie en .ps:
  21. lmess = faux ; lgraph = faux ;
  22. si lgraph ; lps=vrai; sinon; lps= faux; finsi;
  23. si lps; opti trac psc; opti ftra 'xfem_gd.ps'; finsi;
  24.  
  25. * Maillage :
  26. * ########
  27. dx1 = 0.06 ;
  28. blim = 1. + 5 * dx1 ; nbr1 = enti (blim / dx1);
  29. * densité sur la hauteur de la bande :
  30. nbh1 = 3 ; hh1 = ((flot nbh1)*dx1);
  31.  
  32. p1 = 'POIN' 0. hh1 ; p2 = 'POIN' blim hh1 ;
  33. pm3 = 'POIN' blim (0.) ; pm4 = 'POIN' 0. (0.) ;
  34. p3 = 'POIN' blim (-1.*hh1) ; p4 = 'POIN' 0. (-1.*hh1) ;
  35.  
  36. lh = droi nbr1 p1 p2 ;
  37. ldh = droi p2 pm3 DINI dx1 DFIN dx1;
  38. lgh = droi p1 pm4 DINI dx1 DFIN dx1;
  39. lm = droi nbr1 pm4 pm3;
  40. ldb = droi pm3 p3 DINI dx1 DFIN dx1;
  41. lgb = droi pm4 p4 DINI dx1 DFIN dx1;
  42. lb = droi nbr1 p4 p3;
  43.  
  44. s1h = 'DALLER' lm ldh lh lgh ;
  45. s1b = 'DALLER' lb ldb lm lgb ;
  46.  
  47. * fissure initiale :
  48. * #################
  49. a0 = 3.*dx1;
  50. pfis0 = -0.01 0.; pfisi = a0 0.;
  51. fis0 = coul (droi 1 pfis0 pfisi) roug;
  52.  
  53. * Modèle et materiau :
  54. * ####################
  55. Ey1 = 2.e5 ; nu1=0.3; rho1 = 7.e-6 ;
  56.  
  57. * CHARGEMENT : déplacement imposé :
  58. * #################################
  59. tf = 1.; NT = ENTI (3.*tf);
  60. dt0 = tf / (flot NT) ; l_t = PROG 0. pas dt0 tf;
  61. ufin = tf ; l_uimp = PROG 0. pas (ufin / (flot NT)) ufin;
  62. evou = evol 'MANU' 'temps' l_t 'U_imp' l_uimp ;
  63. coeff0 = 0.2;
  64.  
  65.  
  66. **###################
  67. ** Calcul MASSIF **
  68. **###################
  69. s1q = s1h ;
  70. elim s1q 0.0001;
  71. * Modèle et materiau :
  72. * ####################
  73. Mod1q = s1q MODE mecanique elastique isotrope ;
  74.  
  75. mat1q = MATE Mod1q 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  76.  
  77.  
  78. * Conditions aux limites, Déplacements imposés :
  79. * ################################################
  80. ptrach = poin s1q proc (p1 plus ((2.*dx1) (-2.*dx1) ));
  81. cl1hx = 'BLOQUE' ptrach 'UX' ; cl1hy = 'BLOQUE' ptrach 'UY' ;
  82. psyms = elem lm comp (poin s1q proc pfisi) pm3;
  83. cl1sym = 'BLOQUE' psyms 'UY' ;
  84. si lgraph ; trac (s1q et (coul psyms roug)); finsi;
  85.  
  86. cl0 = (cl1hx 'ET' cl1hy 'ET' cl1sym) ;
  87.  
  88. dep1 = DEPI cl1hy coeff0;
  89. uimp0 = CHAR 'DIMP' dep1 evou ;
  90.  
  91. * Resolution iterative du probleme
  92. * ################################################
  93. TAB0 = TABL; TAB0.TEMPS0 = 0. ;
  94. TAB0.'MODELE' = Mod1q ;
  95. TAB0.'CARACTERISTIQUES' = mat1q ;
  96. TAB0.'CHARGEMENT' = uimp0;
  97. TAB0.'BLOCAGES_MECANIQUES' = cl0 ;
  98. TAB0.'TEMPS_CALCULES' = L_t ; TAB0.'TEMPS_SAUVES' = L_t ;
  99. TAB0.'CONVERGENCE_FORCEE' = faux;
  100. TAB0.'DELTAITER' = 190;
  101. TAB0.'MAXITERATION' = 199;
  102.  
  103. PASAPAS tab0;
  104.  
  105. *opti donn 5;
  106.  
  107. **##################
  108. ** Calcul XFEM **
  109. **##################
  110. s1 = s1h et s1b;
  111. elim s1 0.0001;
  112.  
  113. si lgraph; trac (s1 et fis0); finsi;
  114.  
  115. * Modèle et materiau :
  116. * ####################
  117. Mod1a = s1 MODE mecanique elastique isotrope xq4R;
  118. Mod1 = Mod1a ;
  119.  
  120. mat1 = MATE Mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  121.  
  122.  
  123. ** Table XFEM : fisssure initiale, critere d'initiation, modele zone cohesive
  124. ** #################################################################
  125. tab_xfem = tabl;
  126. TAB_xfem.fissure = fis0 ;
  127. * creation des level set
  128. psi1 phi1 = PSIPHI s1 fis0 'DEUX' pfisi;
  129. * Appel a TRIELE :
  130. cheX1 = trie mod1 psi1 phi1 'SAUT';
  131.  
  132.  
  133. * Conditions aux limites, Déplacements imposés :
  134. * ################################################
  135. ptrach = poin s1 proc (p1 plus ((2.*dx1) (-2.*dx1) ));
  136. ptracb = poin s1 proc (p4 plus ((2.*dx1) (2.*dx1) ));
  137. cl1hx = 'BLOQUE' ptrach 'UX' ; cl1bx = 'BLOQUE' ptracb 'UX' ;
  138. cl1hy = 'BLOQUE' ptrach 'UY' ; cl1by = 'BLOQUE' ptracb 'UY' ;
  139.  
  140. cl1 = (cl1hx 'ET' cl1bx 'ET'cl1hy 'ET' cl1by) ;
  141.  
  142. * CHARGEMENT : déplacement imposé :
  143. dep1 = (DEPI cl1hy coeff0) et (DEPI cl1by (-1.*coeff0)) ;
  144. uimp0 = CHAR 'DIMP' dep1 evou ;
  145.  
  146.  
  147. * Resolution iterative du probleme
  148. * ################################################
  149. TAB1 = TABL; TAB1.TEMPS0 = 0. ;
  150. TAB1.'MODELE' = Mod1 ;
  151. TAB1.'CARACTERISTIQUES' = mat1 ;
  152. TAB1.'CHARGEMENT' = uimp0;
  153. TAB1.'BLOCAGES_MECANIQUES' = cl1 ;
  154. TAB1.'TEMPS_CALCULES' = L_t ; TAB1.'TEMPS_SAUVES' = L_t ;
  155. TAB1.'CONVERGENCE_FORCEE' = faux;
  156. TAB1.'DELTAITER' = 190;
  157. TAB1.'MAXITERATION' = 199;
  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.  

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