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

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