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.  
  160. * constructionsion des blocages des ddl X-fem non actifs dans
  161. * les éléments de transition et des relation de conformitées dues au
  162. * raffinent.
  163. Rel1X = tabl;
  164. mod1tot = mod1x et mod1f;
  165. Rel1X.0 = RELA mod1tot;
  166. list resu Rel1X.0;
  167. *list resu Che1X . 0;
  168. * apres cela on assemble
  169. mod1tot = mod1x et mod1f;
  170.  
  171. * on crée le ou les materiau et on les assemble
  172. Ey1 = 2.E11 ;
  173. nu1 = 0.3;
  174. rho1 = 7800.;
  175. mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  176.  
  177.  
  178. *******************************************************
  179. *** Matrices + CL
  180.  
  181. Rig1tot = RIGI mod1tot mat1tot;
  182. list resu RIG1TOT;
  183. * on bloque les depl de corps rigides
  184.  
  185. pb1 = surfc POIN 'PROC' (0.0 Lb1);
  186. pb2 = surfc POIN 'PROC' (0.0 (-1.0*Lb1));
  187. pb3 = surfc POIN 'PROC' (Lb1 0.0);
  188. pb4 = surfc POIN 'PROC' ((-1.0*Lb1) 0.0);
  189.  
  190. cl1x = BLOQ 'UX' (pb1 et pb2);
  191. cl1y = BLOQ 'UY' (pb3 et pb4);
  192.  
  193. cl1tot = cl1x et cl1y;
  194.  
  195. K1tot = Rig1tot et cl1tot et Rel1X.0 ;
  196.  
  197. syy0 = 1.E9;
  198.  
  199. pre1 = PRES 'MASSIF' mod1tot (L4 et L2) (-1.*syy0);
  200.  
  201. si (graph);
  202. TRAC (VECT pre1 'FORC' 'BLEU') (stot2 et lcrack0) TITR 'Chargement';
  203. finsi;
  204. list resu cl1tot;
  205.  
  206. list resu Rel1X.0;
  207.  
  208. *******************************************************
  209. *** RESOlution avec RESO
  210.  
  211. u0 = RESO K1tot pre1;
  212.  
  213. *******************************************************
  214. *** Post traitement
  215.  
  216. * verif du fonctionnement des principaux operateurs
  217. defel0 = EPSI mod1tot u0;
  218. sig0 = SIGM mod1tot u0 mat1tot;
  219. svm0 = VMIS sig0 mod1tot;
  220. gru0 = GRAD mod1tot u0;
  221. def0 = ELAS mod1tot sig0 mat1tot;
  222. ddef0 = def0 - defel0;
  223.  
  224. * calcul du residu
  225. Fint0 = BSIG sig0 mod1tot;
  226. Freac0 = REAC cl1tot u0;
  227. Res0 = Fint0 - (Freac0 + pre1);
  228.  
  229. * * calcul energie elastique = energie exterieure
  230. * wdefel0= 0.5 * (ENER sig0 defel0 mod1tot);
  231. * Iwdef0 = INTG wdefel0 mod1tot;
  232. * wext = 0.5 * (PSCA u0 pre1 (mots UY) (mots FY));
  233. * Iwext = maxi (RESU wext);
  234. * * travail des efforts internes de deformations
  235. * wdef0 = WORK mod1tot sig0 gru0;
  236.  
  237. * messages
  238. mess 'erreur relative / Residu =' ((maxi Res0) / (maxi pre1));
  239. * mess 'erreur relative / Energie =' ((Iwdef0-Iwext) / (abs Iwext));
  240. mess 'ecart relatif EPSI / SIGM+ELAS= ' (maxi ddef0 abs) ;
  241.  
  242. * traces
  243. u0phy = XFEM 'RECO' u0 mod1tot ;
  244.  
  245. si (graph);
  246. def0 = DEFO u0phy stot2 20.;
  247. TRAC sig0 mod1tot def0 'TITR' 'Contraintes';
  248.  
  249. ucrk0u ucrk0d = XFEM 'FISS' lcrack0 u0 mod1x ;
  250.  
  251.  
  252. def0 = def0 et (DEFO ucrk0u lcrack0 20. coul rouge)
  253. et (DEFO ucrk0d lcrack0 20. coul vert);
  254. TRAC def0 TITR 'DEFORMEE';
  255. * trac defel0 mod1tot (lcrack0);
  256. * TRAC svm0 mod1tot (stot1 et lcrack0);
  257. * TRAC wdefel0 mod1tot (stot1 et lcrack0);
  258. * TRAC gru0 mod1tot (stot1 et lcrack0);
  259. finsi;
  260.  
  261.  
  262. *******************************************************
  263. * calcul de l integrale J, de K1, de K2
  264.  
  265. *** valeur theorique de J,K1,K2 ***
  266. k0th1 = syy0 * ((a0 * pi) ** 0.5) * (cos alpha0);
  267. Estar1 = (1. - (nu1**2)) / Ey1;
  268. Gth1 = Estar1 * (k0th1 ** 2);
  269.  
  270. * facteur de forme pour a/b=0.1 et h/b=Inf
  271. *facfor1 = 1.006;
  272. * facteur de forme pour a/b=0.1 et h/b=1
  273. facfor1 = 1.014;
  274. Gth1 = (facfor1 ** 2) * Gth1;
  275. K1ref = k0th1 * (cos alpha0) * facfor1;
  276. K2ref = k0th1 * (sin alpha0) * facfor1;
  277.  
  278.  
  279.  
  280. *** G_THETA -> J ***
  281. SUPTAB = TABL ;
  282. SUPTAB.'OBJECTIF' = MOT 'J';
  283. SUPTAB.'PSI' = psi0;
  284. SUPTAB.'PHI' = phi0;
  285. SUPTAB.'FRONT_FISSURE' = ptip1 ;
  286. SUPTAB.'MODELE' = mod1tot;
  287. SUPTAB.'CARACTERISTIQUES' = mat1tot;
  288. SUPTAB.'SOLUTION_RESO' = u0;
  289. SUPTAB.'CHARGEMENTS_MECANIQUES' = pre1;
  290. * on calcule J pour differents contours
  291. prj1 = prog ;
  292. prcou1= prog ;
  293. lecou1 = lect 1 pas 1 5;
  294. ncou1 = dime lecou1;
  295. jmoy1 = 0.;
  296. repe BCOU1 ncou1;
  297. icou1 = extr lecou1 (&BCOU1);
  298. SUPTAB.'COUCHE' = icou1;
  299. G_THETA SUPTAB;
  300. prj1 = prj1 et (prog (SUPTAB . 'RESULTATS'));
  301. prcou1=prcou1 et (prog (flot icou1));
  302. jmoy1=jmoy1 + (SUPTAB . 'RESULTATS');
  303. fin BCOU1;
  304.  
  305. evjcou1 = evol BLEU manu 'CONTOUR' prcou1 'J' prj1 ;
  306.  
  307. si(graph);
  308. dess evjcou1 TITR 'J en fonction du contour';
  309. fins;
  310.  
  311. * extraction d un J moyen
  312. jmoy1 = jmoy1 / ncou1;
  313.  
  314.  
  315. *erreur relative / J
  316. errG1 = abs (1. - (jmoy1 / Gth1));
  317.  
  318.  
  319. *** Calcul des FIC par methode integrale G_THETA ***
  320. KTAB = TABL;
  321. KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
  322. KTAB .'PSI' = psi0;
  323. KTAB .'PHI' = phi0;
  324. KTAB .'FRONT_FISSURE' = ptip1 ;
  325. KTAB .'MODELE' = mod1tot;
  326. KTAB .'CARACTERISTIQUES' = mat1tot;
  327. KTAB .'SOLUTION_RESO' = u0;
  328. KTAB .'CHARGEMENTS_MECANIQUES' = pre1;
  329. KTAB . 'COUCHE' = ncou1;
  330.  
  331. G_THETA KTAB;
  332.  
  333. K1G = KTAB . 'RESULTATS' . 'I';
  334. K2G = KTAB . 'RESULTATS' . 'II';
  335. errK1G = ABS (1 - (ABS ( K1G/K1ref)));
  336. errK2G = ABS (1 - (ABS ( K2G/K2ref)));
  337. GK1K2 = Estar1 * ((K1G**2) + (K2G**2)) ;
  338.  
  339.  
  340. *********************************************************
  341.  
  342. mess '------------------------------------------------';
  343. mess (chai ' G K1 K2 ');
  344. mess (chai 'reference ' Gth1 ' ' K1ref ' ' K2ref );
  345. mess (chai ' G_THETA ' jmoy1 ' ' K1G ' ' K2G );
  346. mess (chai ' ERR_REL ' errG1 ' ' errK1G ' ' errK2G );
  347.  
  348.  
  349. Si (errG1 < 0.01);
  350. 'MESS' ' ----------------------' ;
  351. 'MESS' ' SUCCES DU CAS-TEST !' ;
  352. 'MESS' ' ----------------------' ;
  353. Erre 0;
  354. Sinon;
  355. 'MESS' ' ---------------------' ;
  356. 'MESS' ' ECHEC DU CAS-TEST !' ;
  357. 'MESS' ' ---------------------' ;
  358. Erre 5;
  359. finsi;
  360.  
  361.  
  362. fin;
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  

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