Télécharger la plaquette Cast3M


fr en

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

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