Télécharger raff03.dgibi

Retour à la liste

Numérotation des lignes :

  1. *******************************************************************
  2. ** raff03.dgibi
  3. *******************************************************************
  4. *
  5. * Calcul en mecanique de la rupture avec un maillage raffine par
  6. * l'operateur RAFF et des elements X-FEM.
  7. *
  8. * d'une plaque elastique en traction en 2D déformations planes
  9. * avec fissure horizontale
  10. * test avec des éléments QUA4 QUA6 TRI3 et TRI6
  11. *
  12. * Dimentions du quart de plaque:
  13. * l1= 0.5 m
  14. * l2= 0.5 m
  15. * Longueur de la fissure :
  16. * a=0.005 m
  17. * Inclinaison :
  18. * alpha=20 degre
  19. *
  20. *création : gg, le 15.03.2017
  21. *
  22. * Comparaison de la simulations avec une solution de reference
  23. * pour g,k1 et k2 donne par :
  24. * [TADA, STRESS ANALYSIS HANDBOOK,1973]
  25. * [ISIDA,IJF(7),301-316,1971]
  26. *
  27. ************************************************************************
  28.  
  29. *** Options de calcul : 2D
  30. opti dime 2 elem qua4 mode plan defo;
  31. *opti trac psc;
  32. *graph = vrai;
  33. graph = faux;
  34.  
  35. *******************************************************
  36. *** Maillage : Plaque rectangulaire saine
  37.  
  38. * dimension et finesse
  39.  
  40.  
  41. Lb1 = 0.50 ;
  42. nb1 = 20;
  43.  
  44. dx2 = (2.*Lb1) / (flot nb1) ;
  45. dx1 = dx2/16;
  46.  
  47.  
  48. P1 = (-1*Lb1) (-1*Lb1) ;
  49. P2 = (-1*Lb1) Lb1 ;
  50. P3 = Lb1 Lb1 ;
  51. P4 = Lb1 (-1*Lb1) ;
  52.  
  53. L1 = DROITE nb1 P1 P2 ;
  54. L2 = DROITE nb1 P2 P3 ;
  55. L3 = DROITE nb1 P3 P4 ;
  56. L4 = DROITE nb1 P4 P1 ;
  57.  
  58. stot1 = DALL L1 L2 L3 L4 ;
  59.  
  60. *******************************************************
  61. *** Maillage de la fissure
  62.  
  63. * longueur et angle de la fissure
  64. *a0 = 0.1;
  65. a0 = 0.05;
  66. *alpha0 = 0.;
  67. alpha0 = 20.;
  68. *alpha0 = 40.;
  69. *alpha0 = 45.;
  70. *alpha0 = 50.;
  71.  
  72. * finesse du maillage de la fissure
  73. *na0 = 10;
  74. *da0 = a0 / (flot na0) ;
  75. da0 = dx1 / 10.;
  76. na0 = enti (a0 / da0);
  77.  
  78. * maillage
  79. DENS (da0) ;
  80. ptip1 = (a0 * (cos alpha0)) (a0 * (sin alpha0));
  81. ptip2 = (-1.*a0 * (cos alpha0)) (-1.*a0 * (sin alpha0));
  82. lcrack0 = (ptip2 droi na0 ptip1) coul 'ROUG';
  83.  
  84.  
  85. * table (+ pratique pour gérer la propa)
  86. TXFEM = tabl;
  87. TXFEM . 0 = tabl;
  88. TXFEM . 0 . 'POINTE1' = ptip1;
  89. TXFEM . 0 . 'POINTE2' = ptip2;
  90. TXFEM . 0 . 'FISSURE' = lcrack0;
  91.  
  92.  
  93.  
  94. *******************************************************
  95. * raffinement du maillage proche des fronts
  96.  
  97. * rayon min de la zone raffiné
  98. dmin= 1.5*dx2;
  99. dmax=dmin+(dx2);
  100.  
  101.  
  102.  
  103. x1 y1 = coor stot1 ;
  104.  
  105. DF1= (((x1 - (a0 * (cos alpha0)))**2)+
  106. ((y1 - (a0 * (sin alpha0)))**2))**0.5;
  107.  
  108. DF2= (((x1 + (a0 * (cos alpha0)))**2)+
  109. ((y1 + (a0 * (sin alpha0)))**2))**0.5;
  110.  
  111. *DF3= (0.5*DF1)+(0.5*DF2);
  112. DF3= mini DF1 DF2;
  113.  
  114.  
  115. aff = (((dx2 - dx1)*(DF3 - dmax))/(dmax - dmin)) + (dx2) ;
  116. dens1 = born aff 'SCAL' 'COMPRIS' (dx1*1.000001) (dx2*1.000001) ;
  117.  
  118. si (graph);
  119. TRAC dens1 stot1 TITR 'Fonction de densite pour RAFF';
  120. finsi;
  121.  
  122. stot2 = raff stot1 dens1 ;
  123.  
  124.  
  125. PTX = STOT2 POIN 'DROI' ptip1 ptip2 (dx2*0.5);
  126. PTX = PTX POIN 'DROI' (0.0 0.0) ((-1.0*(sin alpha0)) (cos alpha0))
  127. (a0+(dx2*0.5));
  128. SURFX = STOT2 ELEM 'APPUYE' 'LARGEMENT' PTX ;
  129. SURFC = STOT2 DIFF SURFX ;
  130.  
  131. si (graph);
  132. TRAC 'FACE' ((SURFC coul jaune) ET (SURFX COUL BLEU) ET lcrack0 )
  133. TITR 'Surface FEM Jaune et X-FEM BLEU';
  134. finsi ;
  135.  
  136. * creation des level set
  137. psi0 phi0 = PSIPHI SURFX lcrack0 'DEUX' ptip1 ptip2;
  138. isolv7 = prog (-2.*a0) PAS (a0 / 4.) (2.*a0);
  139.  
  140. si (graph);
  141. TRAC (SURFX et lcrack0) TITR 'ELEMENTS X-FEM';
  142. TRAC phi0 (SURFX et SURFC ) isolv7 ;
  143. TRAC psi0 (SURFX et SURFC ) isolv7;
  144. finsi ;
  145.  
  146.  
  147. *opti donn 5 ;
  148. *******************************************************
  149. *** Modeles et materiaux
  150.  
  151. * definition des models elementaires
  152. mod1x = MODE surfx mecanique elastique EPSI LINEAIRE xq4R;
  153. mod1f = MODE surfc mecanique elastique EPSI LINEAIRE;
  154.  
  155. *mod1tot = MODE stot2 mecanique elastique isotrope xq4R;
  156. * enrichissement
  157. Che1X = tabl;
  158. Che1X . 0 = TRIE mod1x psi0 phi0;
  159. *list resu Che1X . 0;
  160. * apres cela on assemble
  161. mod1tot = mod1x et mod1f;
  162.  
  163. * on crée le ou les materiau et on les assemble
  164. Ey1 = 2.E11 ;
  165. nu1 = 0.3;
  166. rho1 = 7800.;
  167. mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  168.  
  169.  
  170. *******************************************************
  171. *** Matrices + CL
  172.  
  173. Rig1tot = RIGI mod1tot mat1tot;
  174.  
  175. * on bloque les depl de corps rigides
  176.  
  177. pb1 = surfc POIN 'PROC' (0.0 Lb1);
  178. pb2 = surfc POIN 'PROC' (0.0 (-1.0*Lb1));
  179. pb3 = surfc POIN 'PROC' (Lb1 0.0);
  180. pb4 = surfc POIN 'PROC' ((-1.0*Lb1) 0.0);
  181.  
  182. cl1x = BLOQ 'UX' (pb1 et pb2);
  183. cl1y = BLOQ 'UY' (pb3 et pb4);
  184.  
  185. cl1tot = cl1x et cl1y;
  186.  
  187. K1tot = Rig1tot et cl1tot ;
  188.  
  189. syy0 = 1.E9;
  190.  
  191. pre1 = PRES 'MASSIF' mod1tot (L4 et L2) (-1.*syy0);
  192.  
  193. si (graph);
  194. TRAC (VECT pre1 'FORC' 'BLEU') (stot2 et lcrack0) TITR 'Chargement';
  195. finsi;
  196.  
  197.  
  198. *******************************************************
  199. *** RESOlution avec RESO
  200.  
  201. u0 = RESO K1tot pre1;
  202.  
  203. *******************************************************
  204. *** Post traitement
  205.  
  206. * verif du fonctionnement des principaux operateurs
  207. defel0 = EPSI mod1tot u0;
  208. sig0 = SIGM mod1tot u0 mat1tot;
  209. svm0 = VMIS sig0 mod1tot;
  210. gru0 = GRAD mod1tot u0;
  211. def0 = ELAS mod1tot sig0 mat1tot;
  212. ddef0 = def0 - defel0;
  213.  
  214. * calcul du residu
  215. Fint0 = BSIG sig0 mod1tot;
  216. Freac0 = REAC cl1tot u0;
  217. Res0 = Fint0 - (Freac0 + pre1);
  218.  
  219. * * calcul energie elastique = energie exterieure
  220. * wdefel0= 0.5 * (ENER sig0 defel0 mod1tot);
  221. * Iwdef0 = INTG wdefel0 mod1tot;
  222. * wext = 0.5 * (PSCA u0 pre1 (mots UY) (mots FY));
  223. * Iwext = maxi (RESU wext);
  224. * * travail des efforts internes de deformations
  225. * wdef0 = WORK mod1tot sig0 gru0;
  226.  
  227. * messages
  228. mess 'erreur relative / Residu =' ((maxi Res0) / (maxi pre1));
  229. * mess 'erreur relative / Energie =' ((Iwdef0-Iwext) / (abs Iwext));
  230. mess 'ecart relatif EPSI / SIGM+ELAS= ' (maxi ddef0 abs) ;
  231.  
  232. * traces
  233. u0phy = XFEM 'RECO' u0 mod1tot ;
  234.  
  235. si (graph);
  236. def0 = DEFO u0phy stot2 20.;
  237. TRAC sig0 mod1tot def0 'TITR' 'Contraintes';
  238.  
  239. ucrk0u ucrk0d = XFEM 'FISS' lcrack0 u0 mod1x ;
  240.  
  241.  
  242. def0 = def0 et (DEFO ucrk0u lcrack0 20. coul rouge)
  243. et (DEFO ucrk0d lcrack0 20. coul vert);
  244. TRAC def0 TITR 'DEFORMEE';
  245. * trac defel0 mod1tot (lcrack0);
  246. * TRAC svm0 mod1tot (stot1 et lcrack0);
  247. * TRAC wdefel0 mod1tot (stot1 et lcrack0);
  248. * TRAC gru0 mod1tot (stot1 et lcrack0);
  249. finsi;
  250.  
  251.  
  252. *******************************************************
  253. * calcul de l integrale J, de K1, de K2
  254.  
  255. *** valeur theorique de J,K1,K2 ***
  256. k0th1 = syy0 * ((a0 * pi) ** 0.5) * (cos alpha0);
  257. Estar1 = (1. - (nu1**2)) / Ey1;
  258. Gth1 = Estar1 * (k0th1 ** 2);
  259.  
  260. * facteur de forme pour a/b=0.1 et h/b=Inf
  261. *facfor1 = 1.006;
  262. * facteur de forme pour a/b=0.1 et h/b=1
  263. facfor1 = 1.014;
  264. Gth1 = (facfor1 ** 2) * Gth1;
  265. K1ref = k0th1 * (cos alpha0) * facfor1;
  266. K2ref = k0th1 * (sin alpha0) * facfor1;
  267.  
  268.  
  269.  
  270. *** G_THETA -> J ***
  271. SUPTAB = TABL ;
  272. SUPTAB.'OBJECTIF' = MOT 'J';
  273. SUPTAB.'PSI' = psi0;
  274. SUPTAB.'PHI' = phi0;
  275. SUPTAB.'FRONT_FISSURE' = ptip1 ;
  276. SUPTAB.'MODELE' = mod1tot;
  277. SUPTAB.'CARACTERISTIQUES' = mat1tot;
  278. SUPTAB.'SOLUTION_RESO' = u0;
  279. SUPTAB.'CHARGEMENTS_MECANIQUES' = pre1;
  280. * on calcule J pour differents contours
  281. prj1 = prog ;
  282. prcou1= prog ;
  283. lecou1 = lect 1 pas 1 5;
  284. ncou1 = dime lecou1;
  285. jmoy1 = 0.;
  286. repe BCOU1 ncou1;
  287. icou1 = extr lecou1 (&BCOU1);
  288. SUPTAB.'COUCHE' = icou1;
  289. G_THETA SUPTAB;
  290. prj1 = prj1 et (prog (SUPTAB . 'RESULTATS'));
  291. prcou1=prcou1 et (prog (flot icou1));
  292. jmoy1=jmoy1 + (SUPTAB . 'RESULTATS');
  293. fin BCOU1;
  294.  
  295. evjcou1 = evol BLEU manu 'CONTOUR' prcou1 'J' prj1 ;
  296.  
  297. si(graph);
  298. dess evjcou1 TITR 'J en fonction du contour';
  299. fins;
  300.  
  301. * extraction d un J moyen
  302. jmoy1 = jmoy1 / ncou1;
  303.  
  304.  
  305. *erreur relative / J
  306. errG1 = abs (1. - (jmoy1 / Gth1));
  307.  
  308.  
  309. *** Calcul des FIC par methode integrale G_THETA ***
  310. KTAB = TABL;
  311. KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
  312. KTAB .'PSI' = psi0;
  313. KTAB .'PHI' = phi0;
  314. KTAB .'FRONT_FISSURE' = ptip1 ;
  315. KTAB .'MODELE' = mod1tot;
  316. KTAB .'CARACTERISTIQUES' = mat1tot;
  317. KTAB .'SOLUTION_RESO' = u0;
  318. KTAB .'CHARGEMENTS_MECANIQUES' = pre1;
  319. KTAB . 'COUCHE' = ncou1;
  320.  
  321. G_THETA KTAB;
  322.  
  323. K1G = KTAB . 'RESULTATS' . 'I';
  324. K2G = KTAB . 'RESULTATS' . 'II';
  325. errK1G = ABS (1 - (ABS ( K1G/K1ref)));
  326. errK2G = ABS (1 - (ABS ( K2G/K2ref)));
  327. GK1K2 = Estar1 * ((K1G**2) + (K2G**2)) ;
  328.  
  329.  
  330. *********************************************************
  331.  
  332. mess '------------------------------------------------';
  333. mess (chai ' G K1 K2 ');
  334. mess (chai 'reference ' Gth1 ' ' K1ref ' ' K2ref );
  335. mess (chai ' G_THETA ' jmoy1 ' ' K1G ' ' K2G );
  336. mess (chai ' ERR_REL ' errG1 ' ' errK1G ' ' errK2G );
  337.  
  338.  
  339. Si (errG1 < 0.01);
  340. 'MESS' ' ----------------------' ;
  341. 'MESS' ' SUCCES DU CAS-TEST !' ;
  342. 'MESS' ' ----------------------' ;
  343. Erre 0;
  344. Sinon;
  345. 'MESS' ' ---------------------' ;
  346. 'MESS' ' ECHEC DU CAS-TEST !' ;
  347. 'MESS' ' ---------------------' ;
  348. Erre 5;
  349. finsi;
  350.  
  351.  
  352. fin;
  353.  
  354.  
  355.  
  356.  

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