Télécharger la plaquette Cast3M


fr en

Télécharger xfem01.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : xfem01.dgibi
  2. ************************************************************************
  3.  
  4. ************************************************************************
  5.  
  6. *******************************************************************
  7. * 
  8. *              xfem01.dgibi
  9. *
  10. * calcul avec elements xfem (XQ4R) 
  11. * d'une plaque elastique en traction
  12. * avec fissure inclinée 
  13. * 
  14. * création : bp, le 24.09.2009
  15. * ajout g_theta : bp, 08.06.2010
  16. *   solution de reference pour g,k1 et k2 donne par :    
  17. *   [TADA, STRESS ANALYSIS HANDBOOK,1973]
  18. *   [ISIDA,IJF(7),301-316,1971] -> facteur de forme
  19. * propagation elementaire
  20. ********************************************************************
  21. opti echo 1 ;
  22. *** Traces 
  23. *graph = faux;
  24. *graph = vrai;
  25. graph = vrai; opti trac PSC ; opti EPTR 5;
  26. *** finesse du maillage (=flfin) et cas test complet (=fllong)
  27. *flfin = faux;
  28. flfin = vrai;
  29. fllong = faux;
  30. *fllong = vrai;
  31.  
  32.  
  33. *******************************************************
  34. *** Options de calcul : 2D
  35. opti dime 2 elem qua4 mode plan defo;
  36.  
  37.  
  38. *******************************************************
  39. *** Maillage : Plaque rectangulaire saine
  40.  
  41. * dimension et finesse
  42. *La1 = 0.1 ;
  43. La1 = 0.15 ;
  44. Lb1 = 0.5 ; 
  45. na1 = 20;
  46. si(flfin);
  47. La1 = 0.2 ;
  48.  na1 = 40;
  49. *na1 = 100;  
  50. * na1 = 200;  
  51. * na1 = 400;
  52. fins;
  53. dx1 = (2.*La1) / (flot na1) ;
  54. nb1 = na1 / 2;
  55. dx2 = (2.*Lb1) / (flot nb1) ;
  56.  
  57. * zone centrale
  58. DENS (dx1) ;
  59. pa1 = (-1.*La1) (-1.*La1) ;
  60. pa2 = (    La1) (-1.*La1) ;
  61. pa3 = (    La1) (    La1) ;
  62. pa4 = (-1.*La1) (    La1) ;
  63. da1 = pa1 droi na1 pa2;
  64. da2 = pa2 droi na1 pa3;
  65. da3 = pa3 droi na1 pa4;
  66. da4 = pa4 droi na1 pa1;
  67. sa1 = dall da1 da2 da3 da4;
  68.  
  69. * zone exterieure (pour avoir un milieu quasi-infini)
  70. nb1 = na1 / 2;
  71. dx2 = (2.*Lb1) / (flot nb1) ;
  72.  
  73. DENS (dx2) ;
  74. pb1 = (-1.*Lb1) (-1.*Lb1) ;
  75. pb2 = (    Lb1) (-1.*Lb1) ;
  76. pb3 = (    Lb1) (    Lb1) ;
  77. pb4 = (-1.*Lb1) (    Lb1) ;
  78.  
  79. pb11 = (-1.*La1) (-1.*Lb1) ;
  80. pb12 = ( 1.*La1) (-1.*Lb1) ;
  81. pb21 = (    Lb1) (-1.*La1) ;
  82. pb22 = (    Lb1) ( 1.*La1) ;
  83. pb31 = (    La1) (    Lb1) ;
  84. pb32 = (-1.*La1) (    Lb1) ;
  85. pb41 = (-1.*Lb1) (    La1) ;
  86. pb42 = (-1.*Lb1) (-1.*La1) ;
  87.  
  88. db1 = pb11 droi na1 pb12;
  89. db2 = pb21 droi na1 pb22;
  90. db3 = pb31 droi na1 pb32;
  91. db4 = pb41 droi na1 pb42;
  92.  
  93. sa1b1 = (regl da1 db1 'DINI' dx1 'DFIN' dx2);
  94. sa2b2 = (regl da2 db2 'DINI' dx1 'DFIN' dx2);
  95. sa3b3 = (regl da3 db3 'DINI' dx1 'DFIN' dx2);
  96. sa4b4 = (regl da4 db4 'DINI' dx1 'DFIN' dx2);
  97. sabtot= sa1b1 et sa2b2 et sa3b3 et sa4b4;
  98.  
  99. dc1 = (cote sa1b1 4) coul bleu;
  100. dd1 = (cote sa1b1 2) coul vert;
  101. dc2 = (cote sa2b2 4) coul bleu;
  102. dd2 = (cote sa2b2 2) coul vert;
  103. dc3 = (cote sa3b3 4) coul bleu;
  104. dd3 = (cote sa3b3 2) coul vert;
  105. dc4 = (cote sa4b4 4) coul bleu;
  106. dd4 = (cote sa4b4 2) coul vert;
  107.  
  108. *1er coin
  109. dc2ext = dc2 plus (pb2 moin pb21);
  110. dd1ext = dd1 plus (pb2 moin pb12);
  111. dsb2 = dc2ext et dd1ext et dc2 et dd1;
  112. elim dsb2 (1.E-4*dx1);
  113. sb2 = dall dc2ext dd1ext dc2 dd1;
  114.  
  115. *2eme coin
  116. dc3ext = dc3 plus (pb3 moin pb31);
  117. dd2ext = dd2 plus (pb3 moin pb22);
  118. dsb3 = dc3ext et dd2ext et dc3 et dd2;
  119. elim dsb3 (1.E-4*dx1);
  120. sb3 = dall dc3ext dd2ext dc3 dd2;
  121.  
  122. *3eme coin
  123. dc4ext = dc4 plus (pb4 moin pb41);
  124. dd3ext = dd3 plus (pb4 moin pb32);
  125. dsb4 = dc4ext et dd3ext et dc4 et dd3;
  126. elim dsb4 (1.E-4*dx1);
  127. sb4 = dall dc4ext dd3ext dc4 dd3;
  128.  
  129. *4eme coin
  130. dc1ext = dc1 plus (pb1 moin pb11);
  131. dd4ext = dd4 plus (pb1 moin pb42);
  132. dsb1 = dc1ext et dd4ext et dc1 et dd4;
  133. elim dsb1 (1.E-4*dx1);
  134. sb1 = dall dc1ext dd4ext dc1 dd4;
  135.  
  136. sbtot = sb1 et sb2 et sb3 et sb4;
  137.  
  138.  
  139. *petite modif pour propager + loin
  140. si(flfin);
  141.   sa1 = sa1 et sa2b2 et sa4b4;
  142.   sb1 = sa1b1 et sa3b3 et sbtot;
  143. sino;
  144. * attention on reutilise sb1 !!!
  145.   sb1 = sabtot et sbtot ; 
  146. fins;
  147.  
  148. stot1 = sa1 et sb1;
  149.  
  150.   mess 'maillage central   ext   total' ;
  151.   mess 'nbno=  ' (nbno sa1) (nbno sb1) (nbno stot1);
  152.   mess 'nbel=  ' (nbel sa1) (nbel sb1) (nbel stot1);
  153.  
  154. *opti donn 5;
  155. *******************************************************
  156. *** Maillage de la fissure
  157.  
  158. * longueur et angle de la fissure
  159. *a0 = 0.1;
  160. a0 = 0.05;
  161. *alpha0 = 0.;
  162. alpha0 = 20.;
  163. *alpha0 = 40.;
  164. *alpha0 = 45.;
  165. *alpha0 = 50.;
  166.  
  167. * finesse du maillage de la fissure
  168. *na0 = 10;
  169. *da0 = a0 / (flot na0) ;
  170. da0 = dx1 / 10.; 
  171. na0 = enti (a0 / da0);
  172.  
  173. * maillage
  174. DENS (da0) ;
  175. ptip1 = (a0 * (cos alpha0)) (a0 * (sin alpha0));
  176. ptip2 = (-1.*a0 * (cos alpha0)) (-1.*a0 * (sin alpha0));
  177. lcrack0 = (ptip2 droi na0 ptip1) coul 'ROUG';
  178.  
  179.  
  180. * table (+ pratique pour gérer la propa)
  181. TXFEM = tabl;
  182. TXFEM . 0 = tabl;
  183. TXFEM . 0 . 'POINTE1' = ptip1;
  184. TXFEM . 0 . 'POINTE2' = ptip2;
  185. TXFEM . 0 . 'FISSURE' = lcrack0;
  186.  
  187. * creation des level set 
  188. psi0 phi0 = PSIPHI sa1 lcrack0 'DEUX' ptip1 ptip2;
  189. isolv7 = prog (-2.*a0) PAS (a0 / 4.) (2.*a0);
  190. si(graph); 
  191. *opti isov lign;
  192.   TRAC  (stot1 et lcrack0)  ;
  193.   TRAC phi0 (stot1 et lcrack0) isolv7 ;
  194.   TRAC psi0 (stot1 et lcrack0) isolv7;
  195. *opti isov surf;
  196. finsi;
  197.  
  198.  
  199.  
  200. *opti donn 5 ;
  201. *******************************************************
  202. *** Modeles et materiaux
  203.  
  204. * definition des models elementaires
  205. mod1a = MODE sa1 mecanique elastique isotrope xq4R 'CONSTITUANT' 'A';
  206. mod1b = MODE sb1 mecanique elastique isotrope 'CONSTITUANT' 'B';
  207.  
  208. * apres cela on assemble
  209. mod1tot = mod1a et mod1b;
  210.  
  211. * enrichissement
  212. Che1X = tabl;
  213. Che1X . 0 = TRIE mod1tot psi0 phi0;
  214. Rel1X = tabl;
  215.  
  216. * constructionsion des blocages des ddl X-fem non actifs dans 
  217. * les éléments de transition.
  218. Rel1X . 0 = RELA mod1tot;
  219.  
  220. * on crée le ou les materiau et on les assemble
  221. Ey1  = 2.E11 ;
  222. nu1  = 0.3;
  223. rho1 = 7800.;
  224. mat1tot = MATE mod1tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  225.  
  226.  
  227. *******************************************************
  228. *** Matrices + CL
  229.  
  230. Rig1tot = RIGI mod1tot mat1tot;
  231.  
  232. * on bloque les depl de corps rigides
  233. pb12 = sb1 POIN 'PROC' (0.5 * (pb1 PLUS pb2));
  234. pb23 = sb1 POIN 'PROC' (0.5 * (pb2 PLUS pb3));
  235. pb34 = sb1 POIN 'PROC' (0.5 * (pb3 PLUS pb4));
  236. pb41 = sb1 POIN 'PROC' (0.5 * (pb4 PLUS pb1));
  237.  
  238. cl1x = BLOQ  'UX' (pb12 et pb34);
  239. cl1y = BLOQ  'UY' (pb41 et pb23);
  240.  
  241. cl1tot = cl1x et cl1y;
  242.  
  243. K1tot = Rig1tot et cl1tot et (Rel1X . 0);
  244.  
  245. syy0 = 1.E9;
  246. dbup = dc4ext et db3 et dd2ext;
  247. dbdo = dd4ext et db1 et dc2ext;
  248. pre1 = PRES 'MASSIF' mod1tot (dbup et dbdo) (-1.*syy0);
  249.  
  250. si(graph); 
  251.   TRAC (VECT pre1 'FORC' 'BLEU') (stot1 et lcrack0);
  252. finsi;
  253.  
  254.  
  255. *******************************************************
  256. *** RESOlution avec RESO
  257.  
  258. u0 = RESO K1tot pre1;
  259.  
  260.  
  261. *******************************************************
  262. *** Post traitement
  263.  
  264. * verif du fonctionnement des principaux operateurs
  265. defel0 = EPSI mod1tot u0;
  266. sig0   = SIGM mod1tot u0 mat1tot;
  267. svm0   = VMIS sig0 mod1tot;
  268. gru0 = GRAD mod1tot u0;
  269. def0 = ELAS mod1tot sig0 mat1tot;
  270. ddef0 = def0 - defel0;
  271.  
  272. * calcul du residu
  273. Fint0 = BSIG sig0 mod1tot;
  274. Freac0 = REAC cl1tot u0;
  275. Res0 = Fint0 - (Freac0 + pre1);
  276.  
  277. * * calcul energie elastique = energie exterieure
  278. * wdefel0= 0.5 * (ENER sig0 defel0 mod1tot);
  279. * Iwdef0 = INTG wdefel0 mod1tot;
  280. * wext = 0.5 * (PSCA u0 pre1 (mots UY) (mots FY));
  281. * Iwext = maxi (RESU wext);
  282. * * travail des efforts internes de deformations
  283. * wdef0 = WORK mod1tot sig0 gru0;
  284.  
  285. * messages
  286.   mess 'erreur relative / Residu =' ((maxi Res0) / (maxi pre1));
  287. *   mess 'erreur relative / Energie =' ((Iwdef0-Iwext) / (abs Iwext));
  288.   mess 'ecart relatif EPSI / SIGM+ELAS= ' (maxi ddef0 abs) ;
  289.  
  290. * traces
  291. u0phy = XFEM 'RECO' u0 mod1tot ;
  292. ucrk0u ucrk0d = XFEM 'FISS' lcrack0 u0 mod1tot ;
  293. si(graph); 
  294.   def0 = DEFO u0phy stot1 20.;
  295.   TRAC sig0 mod1tot def0;
  296.   def0 = def0 et (DEFO ucrk0u lcrack0 20.)
  297.               et (DEFO ucrk0d lcrack0 20.);
  298.   TRAC def0;
  299. *  trac defel0 mod1tot (lcrack0);
  300. *  TRAC svm0 mod1tot  (stot1 et lcrack0);
  301. *  TRAC wdefel0 mod1tot (stot1 et lcrack0);
  302. *  TRAC gru0 mod1tot (stot1 et lcrack0);
  303. * test de l'option CHAN CHPO SUPP :
  304.   sig0_p = CHAN 'CHPO' sig0 mod1a 'SUPP';
  305.   geo1 = EXTR sig0_p 'MAILLAGE' ;
  306.   TRAC (EXCO sig0_p 'SMYY') geo1 (sa1 COUL gris);
  307.   TRAC (EXCO sig0 'SMYY') mod1a;
  308. finsi;
  309.  
  310.  
  311.  
  312. * OPTI DONN 5 echo 1 trac X;
  313. *******************************************************
  314. * calcul de l integrale J, de K1, de K2
  315.  
  316. ***      valeur theorique de J,K1,K2           ***
  317. k0th1 = syy0 * ((a0 * pi) ** 0.5) * (cos alpha0);
  318. Estar1 = (1. - (nu1**2)) / Ey1;
  319. Gth1  = Estar1 * (k0th1 ** 2);
  320. * facteur de forme pour a/b=0.1 et h/b=Inf
  321. *facfor1 = 1.006;
  322. * facteur de forme pour a/b=0.1 et h/b=1
  323. facfor1 = 1.014;
  324. Gth1 = (facfor1 ** 2) * Gth1;
  325. K1ref = k0th1 * (cos alpha0) * facfor1;
  326. K2ref = k0th1 * (sin alpha0) * facfor1;
  327.  
  328. ***      G_THETA   ->   J         ***
  329. SUPTAB = TABL ;
  330. SUPTAB.'OBJECTIF' = MOT 'J';
  331. SUPTAB.'PSI' = psi0;
  332. SUPTAB.'PHI' = phi0;
  333. SUPTAB.'FRONT_FISSURE' = ptip1 ; 
  334. SUPTAB.'MODELE' = mod1tot;
  335. SUPTAB.'CARACTERISTIQUES' = mat1tot;
  336. SUPTAB.'SOLUTION_RESO' = u0;
  337. SUPTAB.'CHARGEMENTS_MECANIQUES' = pre1;
  338.  
  339. *on teste toutes les options
  340. si(fllong);
  341.  
  342.   SUPTAB.'COUCHE' = 0;
  343.   G_THETA SUPTAB;
  344. * on teste l'option ou l on choisit soi meme le champ theta
  345.   q7 = (SUPTAB . 'CHAMP_THET1') ;
  346.   si(graph);
  347.     vq7 = VECT q7 BLEU ;
  348.     trac vq7 (stot1 et lcrack0);
  349.   finsi;
  350.   OTER SUPTAB 'COUCHE';
  351.   SUPTAB . 'CHAMP_THETA' = q7;
  352.   G_THETA SUPTAB;
  353. * on calcule J pour differents contours
  354.   prj1 = prog (SUPTAB . 'RESULTATS');
  355.   prcou1= prog 0.;
  356.   si(flfin);
  357.     si (na1 > 10);  lecou1 = lect 1 pas 1 15;
  358.     sino;           lecou1 = lect 1 pas 1 (na1/4);
  359.     fins;
  360.   sino;
  361.       lecou1 = lect 1 pas 1 3;
  362.   fins;
  363.   ncou1 = dime lecou1;
  364.   jmoy1 = 0.;
  365.   repe BCOU1 ncou1;
  366.     icou1 = extr lecou1 (&BCOU1);
  367.     SUPTAB.'COUCHE' = icou1;
  368.     G_THETA SUPTAB; 
  369.     prj1 = prj1 et (prog (SUPTAB . 'RESULTATS'));
  370.     prcou1=prcou1 et (prog (flot icou1));
  371.     jmoy1=jmoy1 + (SUPTAB . 'RESULTATS');
  372.   fin BCOU1;
  373.   evjcou1 = evol BLEU manu 'CONTOUR' prcou1 'J' prj1 ;
  374.   si(graph); dess evjcou1; fins;
  375. * extraction d un J moyen
  376.   jmoy1 = jmoy1 / ncou1;
  377. *cas rapide (on ne teste poas toutes les options)  
  378. sino;
  379.   si(flfin);    ncou1 = 3;
  380.   sino;         ncou1 = 1;
  381.   fins;   
  382.   SUPTAB.'COUCHE' =ncou1;
  383.   G_THETA SUPTAB;
  384.   jmoy1 = SUPTAB . 'RESULTATS';
  385. fins;
  386. *erreur relative / J   
  387.  errG1 = abs (1. - (jmoy1 / Gth1));
  388.  
  389.  
  390. ***   Calcul des FIC par methode integrale G_THETA   ***
  391. KTAB = TABL;
  392. KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
  393. KTAB .'PSI' = psi0;
  394. KTAB .'PHI' = phi0;
  395. KTAB .'FRONT_FISSURE' = ptip1 ; 
  396. KTAB .'MODELE' = mod1tot;
  397. KTAB .'CARACTERISTIQUES' = mat1tot;
  398. KTAB .'SOLUTION_RESO' = u0;
  399. KTAB .'CHARGEMENTS_MECANIQUES' = pre1;
  400.    ncou1 = 1;
  401. KTAB . 'COUCHE' = ncou1;
  402.  
  403.   G_THETA KTAB;
  404.  
  405.   K1G = KTAB . 'RESULTATS' . 'I';
  406.   K2G = KTAB . 'RESULTATS' . 'II';
  407.   errK1G = ABS (1 - (ABS ( K1G/K1ref)));
  408.   errK2G = ABS (1 - (ABS ( K2G/K2ref)));
  409.   GK1K2 = Estar1 * ((K1G**2) + (K2G**2)) ;
  410.  
  411.  
  412. *******************************************************
  413. *** erreur? ***
  414.  
  415. mess '------------------------------------------------';
  416. mess (chai '                   G            K1           K2   ');
  417. mess (chai 'reference     ' Gth1  '  ' K1ref '  ' K2ref );
  418. mess (chai '   G_THETA    ' jmoy1 '  ' K1G   '  ' K2G   );
  419. mess (chai '   ERR_REL    ' errG1 '  ' errK1G '  ' errK2G );
  420.  
  421. *OPTI DONN 5 echo 1 trac X;
  422. * Si (errG1 < 0.05);
  423. *   Erre 0;
  424. * Sinon;
  425. *   Erre 5;
  426. * Finsi;
  427.  
  428.  
  429. *OPTI DONN 5 echo 1 trac X;
  430. lcrack1 = lcrack0;
  431.  
  432. *******************************************************
  433. * BOUCLE TEMPORELLE POUR FAIRE PROPAGER LA FISSURE
  434. *******************************************************
  435. it1 = 0 ;
  436. si(flfin);
  437.   si(fllong); nt1 = 18;
  438.   sino; nt1 = 3;
  439.   fins;
  440. sino;
  441.   si(fllong); nt1 = 7;
  442.   sino; nt1 = 3;
  443.   fins;
  444. fins;
  445.  
  446. REPE BT1 nt1;
  447.   it1 = it1 + 1;
  449.   mess '***********************' it1 '*************************';
  450.   TITR (chai 'iteration ' it1);
  451.  
  452. *******************************************************
  453. * propagation élémentaire de fissure
  454.  
  455. *** Increment et angle de propagation ***
  456. * increment manuel
  457. *  da1 = 1.5 * dx1;
  458. * increment auto
  459.   pptip1 = poin sa1 'PROCH' ptip1;
  460.   dx1tip = coor 3 pptip1;
  461.   da1 = 1.5 * dx1tip;
  462. * angle
  463.   si((abs K2G) &lt;eg (1.E-6*(abs K1G))); 
  464.     tet1c = 0.;
  465.   sino;
  466.     tet0c = K1G / K2G;
  467.     tet8c = (tet0c**2) + 8.;
  468.     tet1c = tet0c - ( (sign K2G) * (tet8c**0.5) );
  469.     tet1c = 2. * (atg (0.25 * tet1c));
  470.   fins;
  471.   mess 'propagation avec l angle ' tet1c;
  472. * les calculs ont été réalisés sur la pointe ptip1
  473. * on recupere la direction du repere local de la fissure (CHPOINT)
  474.   ETIP1 = KTAB . 'UTILTET1' . 'DIRECTION1';
  475.   ETIP2 = KTAB . 'UTILTET1' . 'DIRECTION2';
  476. * on fait tourner de l'angle critique 
  477.   DA1VEC = ((cos tet1c) * ETIP1) + ((sin tet1c) * ETIP2);
  478.   DA1VEC = da1 * DA1VEC;
  479. * par symetrie on en deduit ptip2
  480.   DA2VEC = DA1VEC DEDU 180. (0. 0.) 'ROTA' 
  481.            (ptip1 MANU 'POI1') (ptip2 MANU 'POI1');  
  482. * propagation du front de fissure
  483.   ptip1 = ptip1 PLUS DA1VEC;
  484.   ptip2 = ptip2 PLUS DA2VEC;
  485.   lcrack1 = ((ptip2 droi lcrack1) droi ptip1) coul 'ROUG';
  486. *stockage
  487.   TXFEM . it1 = tabl;
  488.   TXFEM . it1 . 'POINTE1' = ptip1;
  489.   TXFEM . it1 . 'POINTE2' = ptip2;
  490.   TXFEM . it1 . 'FISSURE' = lcrack1;
  491.  
  492. *** Actualisations ***
  493. * des level set
  494.   psi1 phi1 = PSIPHI sa1 lcrack1 'DEUX' ptip1 ptip2;
  495.   si(graph);
  496. *opti isov lign;
  497.     TRAC phi1 (sa1 et lcrack1) isolv7 ;
  498.     TRAC psi1 (sa1 et lcrack1) isolv7;
  499. *opti isov surf;
  500.   finsi;
  501. * du modele, de la rigidite,... (inutile pour le materiau)
  502.   Che1X . it1 = TRIE mod1tot psi1 phi1;
  503. * des blocages des ddl X-fem non actifs dans 
  504. * les éléments de transition.
  505.   Rel1X . it1 = RELA mod1tot;
  506.   Rig1tot = RIGI mod1tot mat1tot;
  507.   K1tot = Rig1tot et cl1tot et ( Rel1X . it1);
  508.  
  509.  
  510. *** resolution ***
  511.   u1 = RESO K1tot pre1;
  512.  
  513.   u1phy = XFEM 'RECO' u1 mod1tot ;
  514.   sig1   = SIGM mod1tot u1 mat1tot;
  515.   ucrk1u ucrk1d = XFEM 'FISS' lcrack1 u1 mod1tot ;
  516.   si(graph);
  517.     def1 = (DEFO u1phy stot1 20.) ;
  518.     TRAC sig1 mod1tot def1;
  519.     def1 = def1 et (DEFO ucrk1u lcrack1 20.)
  520.                 et (DEFO ucrk1d lcrack1 20.);
  521.     TRAC def1;
  522.   finsi;
  523.  
  524. ***   Calcul des FIC par methode integrale G_THETA   ***
  525.   KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
  526.   KTAB .'PSI' = psi1;
  527.   KTAB .'PHI' = phi1;
  528.   KTAB .'FRONT_FISSURE' = ptip1 ; 
  529.   KTAB .'MODELE' = mod1tot;
  530.   KTAB .'CARACTERISTIQUES' = mat1tot;
  531.   KTAB .'SOLUTION_RESO' = u1;
  532.   KTAB .'CHARGEMENTS_MECANIQUES' = pre1;
  533.   KTAB . 'COUCHE' = ncou1;
  534. *  repe bcou 6;
  535. *   KTAB . 'COUCHE' = (&bcou) - 1;
  536.   G_THETA KTAB;
  537.   K1G = KTAB . 'RESULTATS' . 'I';
  538.   K2G = KTAB . 'RESULTATS' . 'II';
  539. *   mess (KTAB . 'COUCHE') K1G K2G;
  540.   si(graph);
  541.     q7 = KTAB . CHAMP_THET1;
  542.     trac (vect q7) (sa1 et lcrack1);
  543.   finsi;
  544. *  fin bcou;
  545.  
  546. *si on a des nan, on arrete tout c'est qu'on est à coté de la plaque
  547. si (non (((abs K1G) < 1.E16) et ((abs K2G) < 1.E16) ));
  548.     quit BT1;
  549. finsi;
  550.  
  551. FIN BT1;
  552. *******************************************************
  553.  
  554.  
  555. * OPTI DONN 5 echo 1 trac X;
  556. *******************************************************
  557. *** erreur? (apres propa elementaire) ***
  558. *******************************************************
  559. pref1 = 9.10932E-02 1.10838E-02;
  560. errptip1 = norm (pref1 moin ptip1);
  561. Si (errptip1 < (1.E-2 * dx1));
  562.   Erre 0;
  563. Sinon;
  564.   Erre 5;
  565. Finsi;
  566.  
  567.  
  568. *******************************************************
  569. ***   autres tests divers relatifs aux XFEM         ***
  570. *******************************************************
  571.  
  572. *** test de la syntaxe 2 de TRIELE ***
  573.  
  574. mod0a = MODE sa1 mecanique elastique isotrope xq4R 'CONSTITUANT' 'A';
  575. TRIE mod0a (che1X . 0);
  576. mod0tot = mod0a et mod1b;
  577. mat0tot = MATE mod0tot 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  578. Rig0tot = RIGI mod0tot mat0tot;
  579. Rel0 = RELA mod0tot;
  580. K0tot = Rig0tot et cl1tot et Rel0;
  581. u0 = RESO K0tot pre1;
  582. sig0 = SIGM mod0tot u0 mat0tot;
  583. def0 = ELAS mod0tot sig0 mat0tot;
  584. u0phy = XFEM 'RECO' u0 mod0tot ;
  585. si(graph); 
  586.   TRAC sig0 mod0tot (DEFO u0phy stot1 20.);
  587. fins;
  588.  
  589. * OPTI DONN 5 echo 1 trac X;
  590. fin;
  591.  
  592. *** test de la sauvegarde des xfem (a prevoir) ***
  593. * opti sauv 'XFEM01.sauv';
  594. * sauv;
  595. * 
  596. * opti rest 'XFEM01.sauv';
  597. * rest;
  598.  
  599.  
  600.  
  601.  
  602.  

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