Télécharger xfem3d_01.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : xfem3d_01.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *******************************************************
  5. *opti echo 0;
  6.  
  7. MESS '' ;
  8. MESS ' CALCUL ELASTIQUE DE PROPAGATION STATIQUE DE FISSURE ';
  9. MESS ' en 3D XFEM ' ;
  10. MESS '' ;
  11.  
  12.  
  13. *******************************************************
  14. *** Options
  15.  
  16. *** options de calcul
  17. opti dime 3 elem CUB8 mode trid;
  18. *** options de trace
  19. graph = faux ;
  20. *graph = vrai ;
  21. opti 'TRAC' PSC ;
  22. *opti 'TRAC' X ;
  23.  
  24.  
  25. *******************************************************
  26. *** Maillage : cube sain 1 x 1 x 1
  27.  
  28. * finesse de la discretisation
  29. * n1 = 36;
  30. * n1 = 24;
  31. * n1 = 12;
  32. n1 = 8;
  33.  
  34. * geoemtrie
  35. L1 = 1. ;
  36. dx1 = (L1) / (flot n1) ;
  37. isov1 = dx1 * (prog -6. PAS 0.5 6.);
  38.  
  39. * face avant (a)
  40. DENS (dx1) ;
  41. pa1 = (-0.5*L1) 0. (-0.5*L1) ;
  42. pa2 = ( 0.5*L1) 0. (-0.5*L1) ;
  43. pa3 = ( 0.5*L1) 0. ( 0.5*L1) ;
  44. pa4 = (-0.5*L1) 0. ( 0.5*L1) ;
  45. dra1 = pa1 droi n1 pa2;
  46.  
  47. * repartition 0.25 // 0.5 // 0.25
  48. vz1 = 0. 0. (0.25*L1);
  49. sa1 = dra1 TRAN (n1/4) vz1;
  50. dra10 = sa1 cote 3;
  51. sa2 = dra10 TRAN (n1/2) (2.*vz1);
  52. dra20 = sa2 cote 3;
  53. sa3 = dra20 TRAN (n1/4) vz1;
  54. dra3 = sa3 cote 3;
  55.  
  56. * volume = cube
  57. ex = 1. 0. 0.;
  58. ey = 0. 1. 0.;
  59. ez = 0. 0. 1.;
  60. ny1 = n1;
  61.  
  62. vol1 = sa2 VOLU 'TRAN' ny1 (L1*ey);
  63. vol2 = sa1 VOLU 'TRAN' ny1 (L1*ey);
  64. vol3 = sa3 VOLU 'TRAN' ny1 (L1*ey);
  65. dx1elim = 1.E-5*dx1;
  66. ELIM vol1 vol2 dx1elim;
  67. ELIM vol1 vol3 dx1elim;
  68. vol23 = vol2 et vol3;
  69. vol0 = vol1 et vol23;
  70.  
  71. * qq recup
  72. pa0 = poin vol1 'PROCH' (0. 0. 0.);
  73. pa0x = poin vol1 'PROCH' (0. L1 0.);
  74. pa0z = (poin vol1 'PROCH' ((-0.5*L1) L1 0.))
  75. et (poin vol1 'PROCH' ((0.5*L1) L1 0.));
  76. sb1 = (dra3 TRAN ny1 (L1*ey)) coul 'BLEU';
  77. ELIM sb1 vol3 (1.E-4 * dx1);
  78. sc1 = (dra1 TRAN ny1 (L1*ey)) coul 'OCEA';
  79. ELIM sc1 vol2 (1.E-4 * dx1);
  80.  
  81. con0 = aret vol0;
  82. env0 = ENVE vol0;
  83.  
  84. eye1 = 100. -50. 30. ;
  85. eye2 = 0. -0. 100. ;
  86. si(graph);
  87. trac eye1 (vol0 et sb1 et sc1) 'CACH';
  88. fins;
  89. mess (nbno vol23) (nbno vol1) (nbno vol0);
  90. mess (nbel vol23) (nbel vol1) (nbel vol0);
  91.  
  92.  
  93. **********************************************************
  94. * tables et indice de propagation
  95. it = 0;
  96. pcrack1 = tabl ;
  97. crack1 = tabl;
  98. che1X = tabl;
  99. Rel1X = tabl;
  100.  
  101.  
  102. ***********************************************************
  103. *** Maillage de la fissure
  104.  
  105. *** fissure semi elliptique aplatie
  106.  
  107. a0 = 0.50;
  108. rapab0 = 0.2;
  109. b0 = a0 * rapab0;
  110. bbord = 0.20;
  111. * finesse du maillage de la fissure
  112. da0 = dx1 / 4.;
  113. na0 = enti (a0 / da0);
  114. nb0 = enti (b0 / da0);
  115. &ELEM = VALE 'ELEM';
  116. opti elem TRI3;
  117. * da000 = rapab0 * da0;
  118. da000 = da0;
  119. * points de base
  120. DENS (da000);
  121. pf1 = a0 0. 0. ;
  122. pf2 = a0 bbord 0. ;
  123. pf4 = (-1.*a0) bbord 0. ;
  124. pf5 = (-1.*a0) 0. 0. ;
  125. DENS (da0);
  126. pf0 = 0. 0. 0.;
  127. pf3 = 0. bbord 0.;
  128. * lignes de base
  129. lf1 = (pf5 droi pf0) droi pf1;
  130. lf2 = (pf4 droi pf3) droi pf2;
  131. * profil elliptique
  132. un2 = manu lf2 'CHPO' 'SCAL' 1.;
  133. x2 = coor lf2 1;
  134. vy2 = b0 * ((un2 - ((x2/a0)**2))**0.5);
  135. vy2 = vy2 nomc 'UY';
  136. lfron1 = (lf2 plus vy2) coul 'ORAN';
  137. * surface reglee
  138. sfron1 = (regl lfron1 lf1) coul 'ROUG';
  139.  
  140. * tres important
  141.  
  142. si(graph);
  143. trac eye2 sfron1;
  144. trac eye1 (con0 et sfron1 ) 'CACH';
  145. trac eye1 (con0 et lfron1) 'CACH';
  146. fins;
  147.  
  148. pcrack1.it = lfron1;
  149. crack1.it = sfron1 ;
  150.  
  151. psi1 phi1 = PSIP vol1 crack1.it 'DEUX' pcrack1.it;
  152.  
  153. si(graph);
  154. trac phi1 vol0 isov1 'COUP' pf0 (0. 0. 1.) pf3;
  155. trac psi1 vol0 isov1 'COUP' pf0 pf1 pf3;
  156. fins;
  157.  
  158.  
  159. *opti donn 5 ;
  160. **********************************************************
  161. *** MODELE & MATERIAU ***
  162. **********************************************************
  163.  
  164. * coef elastique
  165. Ey0 = 2.E11;
  166. nu0 = 0.3 ;
  167. rho0= 7800. ;
  168. mu0 = (0.5 * Ey0) / (1. + nu0) ;
  169. Estar0 = Ey0 / (1. - (nu0**2));
  170.  
  171. * Loi de Paris (on convertit C : mm -> m)
  172. mParis = 3.;
  173. CParis = 1.E-9 * 1.E-3;
  174.  
  175.  
  176. **********************************************************
  177. *** Modele & materiau ***
  178. *
  179. * zone de propagation (X-FEM)
  180. mod1 = MODE vol1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'XC8R';
  181. mat1 = MATE mod1 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
  182. che1X.it = TRIE mod1 psi1 phi1 ;
  183.  
  184. * constructionsion des blocages des ddl X-fem non actifs dans
  185. * les éléments de transition.
  186. Rel1x . it = RELA mod1;
  187. *
  188. * elements standards
  189. mod23 = MODE vol23 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'CUB8';
  190. mat23 = MATE mod23 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
  191. *
  192. mod1tot = mod1 et mod23;
  193. mat1tot = mat1 et mat23 ;
  194.  
  195.  
  196.  
  197. **********************************************************
  198. *** CL et DEPLACEMENTS IMPOSES ***
  199. **********************************************************
  200. syy0 = -100.E6;
  201.  
  202.  
  203. *** MODE I ***
  204. * CL mode I
  205. cl1 = bloq pa0 'UX' 'UY' 'UZ';
  206. cl2 = bloq pa0x 'UZ';
  207. cl3 = bloq pa0z 'UX' 'UZ';
  208. cl1tot = (cl1 ET cl2 ET cl3) ;
  209. * chargement mode I
  210. for1 = PRES 'MASS' mod23 sb1 syy0;
  211. for2 = PRES 'MASS' mod23 sc1 syy0;
  212. for1tot = for1 et for2;
  213.  
  214. si(graph);
  215. trac (vect for1tot 'FORC' 'BLEU')
  216. (con0 et pcrack1.it) 'CACH';
  217. fins;
  218.  
  219.  
  220. *** matrice elastique et assemblage
  221. Rig1tot = RIGI mod1tot mat1tot;
  222. K1tot = Rig1tot et cl1tot et (Rel1x . it);
  223.  
  224.  
  225. *******************************************************
  226. * INITIALISATIONS POUR LA BOUCLE TEMPORELLE
  227. *******************************************************
  228. logda = FAUX;
  229. u1 = tabl;
  230. u1 . 0 = RESO K1tot (0.*for1tot);
  231.  
  232. TKchpo1 = tabl; TK1ev = tabl; TK2ev = tabl; TK3ev = tabl;
  233. *******************************************************
  234. * BOUCLE TEMPORELLE POUR FAIRE PROPAGER LA FISSURE
  235. *******************************************************
  236. *nt1=n1;
  237. nt1=2;
  238.  
  239. REPE BT1 nt1;
  240.  
  241. it = it + 1;
  242. mess '***********************' it '*************************';
  243. TITR (chai 'iteration ' it);
  244.  
  245. *******************************************************
  246. * PROPAGATION ÉLÉMENTAIRE DE FISSURE
  247. *******************************************************
  248. si(logda);
  249.  
  250. *** Angle de propagation ***
  251. un1 = manu 'CHPO' pcrack1 . (it - 1) 'SCAL' 1.;
  252. * angle
  253. sgnK2 = (masq K2G 'EGSUPE' 0.) - (masq K2G 'INFERIEUR' 0.);
  254. * modif pour etre tranquille
  255. masK2 = masq (abs K2G) 'EGINFE' (1.E-4*(abs K1G));
  256. K2G = (masK2) + ((un1 - masK2) * K2G);
  257. tet0c = K1G / K2G;
  258. tet8c = (tet0c**2) + (8. * un1);
  259. tet1c = tet0c - ( sgnK2 * (tet8c**0.5) );
  260. tet1c = 2. * (atg (0.25 * tet1c));
  261. tet1c = (0. * masK2 ) + ((un1 - masK2) * tet1c);
  262. * on recupere la direction du repere local de la fissure (CHPOINT)
  263. ETIP1 = KTAB . 'UTILTET1' . 'DIRECTION1';
  264. ETIP2 = KTAB . 'UTILTET1' . 'DIRECTION2';
  265. * on fait tourner de l'angle critique
  266. DA1VEC = ((cos tet1c) * ETIP1) + ((sin tet1c) * ETIP2);
  267. *trac (vect DA1VEC) (pcrack1 . it et con0);
  268.  
  269. *** vitesse de propagation ***
  270. K1eff = ( (K1G**2) + (K2G**2) )**0.5;
  271. K1eff = ( (K1G**2) + (K2G**2) + ((Estar0/mu0)*(K3G**2)) )**0.5;
  272. * Passage en MPa.m^0.5
  273. K1eff = K1eff * 1.E-6;
  274. * quel est le dn correspondant a da1min ?
  275. K1min= mini K1eff 'ABS';
  276. dadn = CParis * (K1min**mParis);
  277. dn = da1min / dadn;
  278. dadn = (CParis * dn) * (K1eff**mParis) ;
  279. DA1VEC = dadn * DA1VEC;
  280.  
  281. *** on cree un maillage par translation ****
  282. trac (vect DA1VEC 1) (pcrack1 . (it - 1) et con0);
  283. pcrk1 = pcrack1 . (it - 1) PLUS DA1VEC;
  284.  
  285. *** quelques chpo utiles ***
  286. * on calcul la distance signee a la paroi
  287. dis0 = PSIP pcrack1 . (it - 1) env0;
  288. dis1 = PSIP pcrk1 env0;
  289. * on cree un chpoint permettant de savoir si un point est a l interieur
  290. unvol1 = manu 'CHML' vol0 'SCAL' 1.;
  291. uncrk1 = PROI pcrk1 unvol1;
  292.  
  293. *** Etape 1 : traitement des segments "a cheval" sur le bord ***
  294. * boucle sur les segments de pcrk1
  295. dismall1 = 1.E-6 * da0;
  296. ncrk1 = nbel pcrk1;
  297. icrk1 = 0;
  298. repe Bcrk1 ncrk1;
  299. icrk1 = icrk1 + 1;
  300. * recup des segments et des points
  301. pcrk1i = ELEM pcrk1 icrk1;
  302. pt1a = poin pcrk1i 'INITIAL'; pt1b = poin pcrk1i 'FINAL';
  303. * test pour savoir si a cheval
  304. un1a = extr uncrk1 pt1a 'SCAL'; un1b = extr uncrk1 pt1b 'SCAL';
  305. si((abs(un1a - un1b)) >eg 0.5);
  306. phi1a = EXTR dis1 pt1a 'PHI'; phi1b = EXTR dis1 pt1b 'PHI';
  307. si(un1a < 0.5);
  308. * il faut ramener le point pt1a
  309. v1a = phi1a / (phi1a - phi1b);
  310. v1a = v1a * (pt1b moin pt1a);
  311. DEPL pt1a 'PLUS' v1a;
  312. sino;
  313. * il faut ramener le point pt1b
  314. v1b = phi1b / (phi1b - phi1a);
  315. v1b = v1b * (pt1a moin pt1b);
  316. DEPL pt1b 'PLUS' v1b;
  317. fins;
  318. fins;
  319. fin Bcrk1;
  320.  
  321. *** Etape 2 : traitement des points encore a l'exterieur ***
  322. * mise a jour de uncrk1
  323. uncrk1 = PROI pcrk1 unvol1;
  324. icrk1 = 0;
  325. repe Bcrk2 ncrk1;
  326. icrk1 = icrk1 + 1;
  327. * recup des segments et des points
  328. pcrk1i = ELEM pcrk1 icrk1;
  329. pt1a = poin pcrk1i 'INITIAL';
  330. * pt1a a l exterieur ?
  331. un1a = EXTR uncrk1 pt1a 'SCAL';
  332. si(un1a < 0.5);
  333. pcrk0i = ELEM pcrack1 . (it - 1) icrk1;
  334. pt0a = poin pcrk0i 'INITIAL';
  335. * distance signée associée
  336. phi0a = EXTR dis0 pt0a 'PHI';
  337. phi1a = EXTR dis1 pt1a 'PHI';
  338. * il faut ramener le point pt1a
  339. v1a = phi1a / (phi0a - phi1a);
  340. v1a = v1a * (pt1a moin pt0a);
  341. DEPL pt1a 'PLUS' v1a;
  342. fins;
  343. fin Bcrk2;
  344. * sans oublier le dernier point du dernier segment
  345. pt1a = poin pcrk1i 'FINAL';
  346. * pt1a a l exterieur ?
  347. un1a = EXTR uncrk1 pt1a 'SCAL';
  348. si(un1a < 0.5);
  349. * recup des segments et des points du front precedent
  350. pcrk0i = ELEM pcrack1 . (it - 1) ncrk1;
  351. pt0a = poin pcrk0i 'FINAL';
  352. * distance signée associée
  353. phi0a = EXTR dis0 pt0a 'PHI';
  354. phi1a = EXTR dis1 pt1a 'PHI';
  355. * il faut ramener le point pt1a
  356. v1a = phi1a / (phi0a - phi1a);
  357. v1a = v1a * (pt1a moin pt0a);
  358. DEPL pt1a 'PLUS' v1a;
  359. fins;
  360.  
  361. *** si tous les points du front sont sortis
  362. * c est qu on a tout cassé !!! => on quitte
  363. meshun1 = extr uncrk1 'MAILLAGE';
  364. si(ega (TYPE meshun1) 'ENTIER'); quit BT1; fins;
  365.  
  366. *** Etape 3 : points sur le bord qui sont entrés a linterieur ***
  367. * mise a jour de dis1
  368. dis1 = PSIP pcrk1 env0;
  369. * 1er point du 1er segment
  370. pcrk0i = ELEM pcrack1 . (it - 1) 1;
  371. pt0a = poin pcrk0i 'INITIAL';
  372. phi0a = EXTR dis0 pt0a 'PHI';
  373. si((abs phi0a) &lt;eg dismall1);
  374. pcrk1i = ELEM pcrk1 1;
  375. pt1a = poin pcrk1i 'INITIAL';
  376. phi1a = EXTR dis1 pt1a 'PHI';
  377. si((abs phi1a) > dismall1);
  378. pt0b = poin pcrk0i 'FINAL';
  379. phi0b= EXTR dis0 pt0a 'PHI';
  380. pt1b = poin pcrk1i 'FINAL';
  381. phi1b= EXTR dis1 pt1b 'PHI';
  382. v1a = phi1a / (phi1b - phi1a) ;
  383. v1a = v1a * (pt1a moin pt1b);
  384. DEPL pt1a 'PLUS' v1a;
  385. fins;
  386. fins;
  387. * dernier point du dernier segment
  388. pcrk0i = ELEM pcrack1 . (it - 1) ncrk1;
  389. pt0b = poin pcrk0i 'FINAL';
  390. phi0b= EXTR dis0 pt0a 'PHI';
  391. si((abs phi0b) &lt;eg dismall1);
  392. pcrk1i = ELEM pcrk1 ncrk1;
  393. pt1b = poin pcrk1i 'FINAL';
  394. phi1b= EXTR dis1 pt1b 'PHI';
  395. si((abs phi1b) > dismall1);
  396. pt0a = poin pcrk0i 'INITIAL';
  397. phi0a = EXTR dis0 pt0a 'PHI';
  398. pt1a = poin pcrk1i 'INITIAL';
  399. phi1a = EXTR dis1 pt1a 'PHI';
  400. v1b = phi1b / (phi1a - phi1b) ;
  401. v1b = v1b * (pt1b moin pt1a);
  402. DEPL pt1b 'PLUS' v1b;
  403. fins;
  404. fins;
  405.  
  406. *** on recupere le vrai nouveau front ***
  407. pcrack1 . it = ELEM pcrk1 'APPUYE' 'STRICT' meshun1;
  408. si((flot (it / 2)) ega ((flot it) / 2.));
  409. pcrack1 . it = pcrack1 . it coul 'ROUG';
  410. sinon;
  411. pcrack1 . it = pcrack1 . it coul 'ORAN';
  412. fins;
  413.  
  414. *** on maille la surface fissuree ***
  415. &ELEM = VALE 'ELEM';
  416. opti 'ELEM' TRI3;
  417. nregl1 = 1 + (enti (da1min / da0));
  418. sfiss1 = regl pcrk1 nregl1 pcrack1 . (it - 1) ;
  419. *** on retraite la surface fissuree cree au debut ***
  420. ELIM sfiss1 dismall1;
  421. sfiss1 = (REGE sfiss1) elem 'TRI3';
  422. crack1 . it = sfiss1 et crack1 . (it - 1) ;
  423.  
  424. *** fin de l avancee du maillage du front **********************
  425.  
  426.  
  427. *** Actualisations ***
  428. * ...des level set
  429. psi1 phi1 = PSIP vol1 crack1.it 'DEUX' pcrack1.it;
  430. si(graph);
  431. trac phi1 vol0 isov1 'COUP' pf0 (0. 0. 1.) pf3;
  432. trac psi1 vol0 isov1 'COUP' pf0 pf1 pf3;
  433. fins;
  434.  
  435. * ...du modele, de la rigidite,... (inutile pour le materiau)
  436. Che1X . it = TRIE mod1tot psi1 phi1;
  437.  
  438. * ... des blocages des ddl X-fem non actifs dans
  439. * les éléments de transition.
  440. Rel1x . it = RELA mod1tot;
  441. Rig1tot = RIGI mod1tot mat1tot;
  442. K1tot = Rig1tot et cl1tot et (Rel1x . it) ;
  443.  
  444.  
  445. *******************************************************
  446. * PAS DE PROPA
  447. *******************************************************
  448. sino;
  449.  
  450. crack1 . it = crack1 . (it - 1);
  451. pcrack1 . it = pcrack1 . (it - 1);
  452.  
  453. fins;
  454.  
  455. trac (con0 et crack1 . it) eye2;
  456.  
  457. ****************************
  458. *** resolution ***
  459. u1 . it = RESO K1tot for1tot;
  460.  
  461. u1phy = XFEM 'RECO' u1 . it mod1tot ;
  462. sig1 = SIGM mod1tot u1 . it mat1tot line;
  463. ucrk1u ucrk1d = XFEM 'FISS' crack1 . it u1 . it mod1tot ;
  464. si(graph);
  465. * amp1 = 25.;
  466. amp1 = 10.;
  467. def1 = (DEFO u1phy vol0 amp1) ;
  468. TRAC sig1 mod1tot def1 eye1;
  469. def2 = (DEFO u1phy con0 amp1)
  470. et (DEFO ucrk1u crack1 . it amp1)
  471. et (DEFO ucrk1d crack1 . it amp1);
  472. TRAC def2 eye1;
  473. def3 = def1
  474. et (DEFO ucrk1u crack1 . it amp1)
  475. et (DEFO ucrk1d crack1 . it amp1);
  476. TRAC def3 eye1 'COUP' pf0 (0. 0. 1.) pf3 'CACH';
  477. finsi;
  478.  
  479. *** Calcul des FIC par methode integrale G_THETA ***
  480. KTAB = tabl ;
  481. KTAB . 'OBJECTIF' = MOT 'DECOUPLAGE';
  482. KTAB .'PSI' = psi1;
  483. KTAB .'PHI' = phi1;
  484. KTAB .'FRONT_FISSURE' = pcrack1 . it ;
  485. KTAB .'MODELE' = mod1tot;
  486. KTAB .'CARACTERISTIQUES' = mat1tot;
  487. KTAB .'SOLUTION_RESO' = u1 . it;
  488. KTAB .'CHARGEMENTS_MECANIQUES' = for1tot;
  489. KTAB . 'COUCHE' = 1;
  490. * KTAB . 'COUCHE' = 2;
  491.  
  492. G_THETA KTAB;
  493.  
  494. Kchpo1 = KTAB . 'CHPO_RESULTATS' ;
  495. TKchpo1 . it = Kchpo1;
  496.  
  497. K1ev = EVOL BLEU 'CHPO' Kchpo1 'K1' pcrack1 . it;
  498. K2ev = EVOL DEFA 'CHPO' Kchpo1 'K2' pcrack1 . it;
  499. K3ev = EVOL ROUG 'CHPO' Kchpo1 'K3' pcrack1 . it;
  500. TK1ev . it = K1ev;
  501. TK2ev . it = K2ev;
  502. TK3ev . it = K3ev;
  503.  
  504. S1pr = EXTR K1ev 'ABSC';
  505. K1pr = EXTR K1ev 'ORDO';
  506. K2pr = EXTR K2ev 'ORDO';
  507. K3pr = EXTR K3ev 'ORDO';
  508.  
  509. * K1eqpr = ((K1pr**2) + (K2pr**2))**0.5;
  510. K1eqpr = ( (K1pr**2) + (K2pr**2) + ((Estar0/mu0)*(K3pr**2)) )**0.5;
  511.  
  512. K1eqev = EVOL 'VERT' 'MANU' 's' S1pr 'K1eq' K1eqpr;
  513. si(it ega 1);
  514. tdess1 = tabl;
  515. tdess1 . 1 = 'MOT' 'MARQ PLUS TIRR ';
  516. tdess1 . 2 = 'MOT' 'MARQ CROI TIRR ';
  517. tdess1 . 3 = 'MOT' 'TIRM ';
  518. fins;
  519. dess (K1ev et K2ev et K3ev et K1eqev) tdess1 'TITR' 'SIF';
  520.  
  521. dadnev = CParis * (K1eqev**mParis);
  522. dess dadnev 'TITR' 'DADN';
  523.  
  524. K1G = EXCO Kchpo1 'K1';
  525. K2G = EXCO Kchpo1 'K2';
  526. K3G = EXCO Kchpo1 'K3';
  527.  
  528. GTAB = tabl ;
  529. GTAB .'OBJECTIF' = MOT 'J ';
  530. GTAB .'PSI' = psi1;
  531. GTAB .'PHI' = phi1;
  532. GTAB .'FRONT_FISSURE' = pcrack1 . it ;
  533. GTAB .'MODELE' = mod1tot;
  534. GTAB .'CARACTERISTIQUES' = mat1tot;
  535. GTAB .'SOLUTION_RESO' = u1 . it;
  536. GTAB .'CHARGEMENTS_MECANIQUES' = for1tot;
  537. GTAB .'COUCHE' = 1;
  538.  
  539. G_THETA GTAB;
  540. Jchpo1 = GTAB . CHPO_RESULTATS;
  541. Jev1 = EVOL ROSE 'CHPO' Jchpo1 'J' pcrack1 . it;
  542. dess Jev1;
  543.  
  544.  
  545. *** Increment de propagation ***
  546. logda = VRAI;
  547. * increment manuel
  548. da1min = 1.4 * dx1;
  549.  
  550.  
  551. * OPTI DONN 5 echo 1 trac X;
  552. FIN BT1;
  553. *******************************************************
  554.  
  555.  
  556. *******************************************************
  557. *** erreur? (avant propa elementaire) ***
  558. *******************************************************
  559. K1ev1 = TK1ev . 1;
  560. K2ev1 = TK2ev . 1;
  561. K3ev1 = TK3ev . 1;
  562. K1ev1m = maxi (extr K1ev1 'ORDO');
  563. K2ev1m = maxi (extr K2ev1 'ORDO');
  564. K3ev1m = maxi (extr K3ev1 'ORDO');
  565. * mess K1ev1m K2ev1m K3ev1m;
  566. K1ref1m = 1.59717E+08;
  567. K2ref1m = 0.;
  568. K3ref1m = 0.;
  569.  
  570. errK1 = abs ((K1ev1m - K1ref1m) / K1ref1m);
  571. errK2 = abs ((K2ev1m - K2ref1m) / K1ref1m) ;
  572. errK3 = abs ((K2ev1m - K2ref1m) / K1ref1m);
  573. errK1 = errK1 + errK2 + errK3;
  574.  
  575. *******************************************************
  576. *** erreur? (apres propa elementaire) ***
  577. *******************************************************
  578. Jmaxref1 = 8.14855E+05;
  579. Jminref1 = 4.75952E+05;
  580. Jmin1 = mini Jchpo1;
  581. Jmax1 = maxi Jchpo1;
  582. errJ1 = (abs ((Jmin1 - Jminref1) /Jminref1))
  583. + (abs ((Jmax1 - Jmaxref1) /Jmaxref1));
  584.  
  585.  
  586. mess 'erreur relative sur K1+K2+K3 et sur J';
  587. mess ' ' errK1 errJ1;
  588. *******************************************************
  589. *2% d ecart acceptee
  590. *******************************************************
  591. Si ((errJ1 < 2.E-2) ou (errK1 < 2.E-2));
  592. Erre 0;
  593. Sinon;
  594. Erre 5;
  595. Finsi;
  596.  
  597.  
  598.  
  599. FIN ;
  600.  
  601.  
  602.  
  603.  
  604.  
  605.  
  606.  
  607.  
  608.  
  609.  
  610.  
  611.  
  612.  
  613.  
  614.  

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