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

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