Télécharger xfem3d_02.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : xfem3d_02.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *******************************************************
  5. **-------------------------------------------------------------------
  6. MESS '' ;
  7. MESS ' CALCUL ELASTIQUE DYNAMIQUE ';
  8. MESS ' en 3D XFEM ' ;
  9. MESS ' Unites: N MPa mm g ms' ;
  10. MESS '' ;
  11. **--------------------------------------------------------------------
  12. *******************************************************
  13. *** Options de calcul
  14. opti dime 3 elem CUB8 mode trid;
  15. * ajout de option epsilon lineaire pour la precision des test!
  16. OPTION epsilon lineaire;
  17. opti 'TRAC' PSC ;
  18. * opti 'TRAC' X ;
  19.  
  20. *******************************************************
  21. *** options de trace
  22.  
  23. isov1 = prog -0.7 pas 0.1 0.7;
  24.  
  25. *******************************************************
  26. *** maillage
  27. pa0 = 0. 0. 0. ; pb0 = 3. 0. 0. ;
  28. liab = pa0 droi 12 pb0;
  29. su1a = liab trans 16 (0. 0. -4.);
  30. liam = COTE 3 su1a;
  31. su1b = liam trans 8 (0. 0. -2.);
  32. liah = COTE 3 su1b;
  33. vol1a = su1a volu tran 12 (0. 6. 0.);
  34. vol1b = su1b volu tran 12 (0. 6. 0.);
  35. vol1 = vol1a et vol1b;
  36. elim 0.001 vol1;
  37. *
  38. *trac cach (vol1 et (liab coul vert) et (liah coul rouge));
  39. *mess (nbno vol1) (nbel vol1) ;
  40. *
  41. pc1 = (enve vol1a) poin proch (3. 6. 0.);
  42. pc2 = (enve vol1b) poin proch (3. 6. -6.);
  43.  
  44. **********************************************************
  45. *** fissure 1 = debouchante (=1 seule pointe) rectiligne
  46. opti elem tri3;
  47. a0 = 3.;
  48. pf1 = 0. 0. -3.;
  49. pf2 = 3. 0. -3.;
  50. pf3 = 3. a0 -3;
  51. pf4 = 0. a0 -3;
  52. ll12 = droi 8 pf1 pf2;
  53. fro23 = droi 10 pf2 pf3;
  54. pcrack2 = droi 8 pf3 pf4;
  55. ll41 = droi 10 pf4 pf1;
  56.  
  57. crack1 = (dall ll12 fro23 pcrack2 ll41 'PLAN') coul roug;
  58.  
  59.  
  60. psi1 phi1 = PSIPHI vol1a crack1 'DEUX' pcrack2;
  61.  
  62. *isov1 = prog -0.001 pas 0.0001 0.001;
  63. * opti isov surf ;
  64. * trac phi1 (vol1a et vol1b) isov1 ((vol1a et vol1b ) et crack1);
  65. * trac psi1 (vol1a et vol1b) isov1 ((vol1a et vol1b ) et crack1);
  66.  
  67. **********************************************************
  68. *** MODELE & MATERIAU ***
  69. **********************************************************
  70. Ey0 = 200000.;
  71. nu0 = 0.3 ;
  72. rho0= 7.8e-3 ;
  73. **********************************************************
  74. *** Modele & materiau ***
  75. *
  76. * elements standards
  77. mod1b = 'MODELISER' vol1b 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'CUB8';
  78. mat1b = 'MATERIAU' mod1b'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
  79. *
  80. **********************************************************
  81. * zone de propagation (X-FEM)
  82. mod1a = 'MODELISER' vol1a 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'XC8R';
  83. * che1X = TRIE Mod1a psi1 phi1 'SAUT' ;
  84. che1X = TRIE Mod1a psi1 phi1 ;
  85. *
  86. modtot = mod1a et mod1b;
  87. mat1a = 'MATERIAU' mod1a 'YOUN' Ey0 'NU' nu0 'RHO' rho0 ;
  88.  
  89. **********************************************************
  90.  
  91.  
  92.  
  93. **********************************************************
  94. *** CL et DEPLACEMENTS IMPOSES ***
  95. **********************************************************
  96.  
  97. cl1 = bloq (pc1 et pc2) 'UZ';
  98. cl2 = bloq liab 'UX' 'UY';
  99. cl3 = bloq liah 'UX' 'UY';
  100. cltot = (cl1 'ET' cl2 'ET' cl3) ;
  101. *
  102. *... chargement
  103. for1 = FORC 'FZ' liah (-3.E3);
  104. for2 = FORC 'FZ' liab (3.E3);
  105. *
  106. xx1 = prog 0. 7.5e-4 10.;
  107. yy1 = prog 0. 1. 1.;
  108. evochar = evol manu 'temps' xx1 'charge' yy1;
  109. char1 = char meca (for1 et for2) evochar;
  110. *
  111. *trac (vect (for1 et for2) 'FX' 'FY' 'FZ' 1.e-2) ((enve vol1)
  112. *et (liah coul vert) et (liab coul rose));
  113. *******************
  114.  
  115. **********************************************************
  116. *** MATRICES ***
  117. **********************************************************
  118.  
  119. *** de rigidite ***
  120. * Xfem
  121. *Rig1a = RIGI Mod1a mat1a;
  122. * standard
  123. *Rig1b = RIGI Mod1b mat1b;
  124.  
  125. *** de masse ***
  126. * Xfem
  127. mas1a = mass mod1a mat1a;
  128. * standard
  129. mas1b = mass mod1b mat1b ;
  130. *
  131. *=======================================================*
  132. * *
  133. * RESOLUTION *
  134. * *
  135. *=======================================================*
  136.  
  137. dt1 = 2.5e-4;
  138. xx0 = prog 0. pas dt1 'NPAS' 4;
  139. *
  140. tab1 = table;
  141. tab1.dynamique = vrai;
  142. tab1.modele = modtot;
  143. tab1.caracteristiques = mat1a et mat1b;
  144. tab1.blocages_mecaniques = cltot;
  145. tab1.chargement = char1;
  146. tab1.temps_calcules = xx0;
  147. pasapas tab1;
  148. *
  149. *post traitement
  150. ind0 = index (tab1.deplacements);
  151. nnn1 = dime ind0;
  152. lwela1 = prog;
  153. lwcin1 = prog;
  154. ltemp1 = prog;
  155. ev1 = table;
  156. def1 = table;
  157. *
  158. ii = 0 ;
  159. repeter bouc1 nnn1;
  160. ii = ii + 1;
  161. pass = ind0.ii;
  162. sig1 = tab1.contraintes.pass;
  163. ut1 = tab1.deplacements.pass;
  164. vt1 = tab1.vitesses.pass;
  165. tt1 = tab1.temps.pass;
  166. etot1 = EPSI ut1 modtot;
  167. *
  168. ev1.ii = evol chpo ut1 'UZ' liah;
  169. def1.ii = defo ut1 vol1 30.;
  170.  
  171. ltemp1 = ltemp1 et (prog tt1);
  172. *...elastique
  173. Wela0 = 'INTG' (0.5 *(ENER sig1 etot1 modtot)) modtot ;
  174. lwela1 = lwela1 et (prog wela0);
  175. *
  176. *...cinetique
  177. Wcin0 = 0.5 * (XTMX vt1 (mas1a 'ET' mas1b)) ;
  178. lwcin1 = lwcin1 et (prog wcin0);
  179. fin bouc1;
  180. *
  181. *
  182. **********************************************************
  183. errmax = 0.;
  184. lerr1= prog 0. -1.20735e-07 2.44871e-06 2.08117e-06
  185. -9.9959e-07 3.38064e-04 9.78166e-04;
  186. lresu = prog 1.e-8 3.49142e-03 1.09155e-02 1.8968e-02 2.44906e-02
  187. 73.8862 27.5316;
  188.  
  189. *comparaison XFEM et STANDARD 3D
  190. ii = 0;
  191. repeter bouc3 nnn1;
  192. ii = ii + 1 ;
  193. ord2 = extr lresu ii;
  194. val0 = extr lerr1 ii;
  195. ord1 = extr ev1.ii 'ORDO';
  196. val1 =( (abs(extr ord1 3)) - ord2)/ ord2;
  197. message 'ord1 ord2 val1' (extr ord1 3) ord2 val1;
  198. si (ii > 1);
  199. message ' valeurs a comparer' val1 val0;
  200. errmax = maxi(prog errmax val1);
  201. finsi;
  202. fin bouc3;
  203. *
  204. evn1 = evol manu 'temps' ltemp1 'ener' lwela1;
  205. evn2 = evol manu 'temps' ltemp1 'ener' lwcin1;
  206.  
  207. ii = ii + 1 ;
  208. ord2 = extr lresu ii;
  209. val0 = extr lerr1 ii;
  210. ord1 = extr evn1 'ORDO';
  211. val1 = abs(( (abs(extr ord1 5)) - ord2)/ ord2);
  212. message 'ord1 ord2 val1' (extr ord1 5) ord2 val1;
  213. message ' valeurs a comparer' val1 val0;
  214. errmax = maxi(prog errmax val1);
  215. *
  216. ii = ii + 1 ;
  217. ord2 = extr lresu ii;
  218. val0 = extr lerr1 ii;
  219. ord1 = extr evn2 'ORDO';
  220. val1 = abs(( (abs(extr ord1 5)) - ord2)/ ord2);
  221. message 'ord1 ord2 val1' (extr ord1 5) ord2 val1;
  222. message ' valeurs a comparer' val1 val0;
  223. errmax = maxi(prog errmax val1);
  224. *
  225. si(errmax >EG 0.001);
  226. ERRE 5;
  227. finsi;
  228. *
  229.  
  230. *opti donn 5;
  231. fin;
  232. **-------------------------------------------------------------------
  233. **-------------------------------------------------------------------
  234. MESS '' ;
  235. MESS ' CALCUL ELASTIQUE DYNAMIQUE ';
  236. MESS ' en 3D standard ' ;
  237. MESS ' Unites: N MPa mm g ms' ;
  238. MESS '' ;
  239. **--------------------------------------------------------------------
  240. *
  241. opti echo 1 ;
  242.  
  243. *******************************************************
  244. *** Options de calcul
  245. opti dime 3 elem CUB8 mode trid;
  246. * opti 'TRAC' PSC ;
  247. opti 'TRAC' X ;
  248. *******************************************************
  249. *** options de tracé
  250. *isov1 = prog -0.7 pas 0.1 0.7;
  251. *******************************************************
  252. *** maillage 1
  253. pa1 = 0. 0. 0. ; pb1 = 3. 0. 0. ;
  254. liab2 = pa1 droi 12 pb1;
  255. *
  256. sy1a = liab2 trans 12 (0. 0. -3.);
  257. bor1 = COTE 3 sy1a;
  258. vol2a = sy1a volu tran 12 (0. 6. 0.);
  259. *
  260. pa2 = 0. 0. -3. ; pb2 = 3. 0. -3. ;
  261. bor2 = pa2 droi 12 pb2;
  262. sy1b = bor2 trans 12 (0. 0. -3.);
  263. vol2b = sy1b volu tran 12 (0. 6. 0.);
  264. *
  265. liah2 = inve ( COTE 3 sy1b);
  266. *
  267. vol2 = vol2a et vol2b;
  268. *
  269. *longueur fissure 3 mm
  270. xfro1 = 2.8;
  271. *
  272. px = poin vol2a plan (1. 0. -3.) (0. 1. -3.) (0. 0. -3.) 0.001;
  273. pla1 = elem (enve vol2a) appuye strictement px;
  274. elim 0.001 (pla1 et bor1);
  275. cha1 = 'COOR' 2 pla1;
  276. geo1 = cha1 'POIN' 'EGSUPE' xfro1;
  277. pla2 = elem pla1 appuye strictement geo1;
  278. *
  279. px = poin vol2b plan (1. 0. -3.) (0. 1. -3.) (0. 0. -3.) 0.001;
  280. plb1 = elem (enve vol2b) appuye strictement px;
  281. elim 0.01 (plb1 et bor2);
  282. chb1 = 'COOR' 2 plb1;
  283. geo2 = chb1 'POIN' 'EGSUPE' xfro1;
  284. plb2 = elem plb1 appuye strictement geo2;
  285. *
  286. **********************************************************
  287. *** fissure 1 = debouchante (=1 seule pointe) rectiligne
  288. pf1 = (cont pla2) poin proch (0. 3. -3.);
  289. pf2 = (cont pla2) poin proch (3. 3. -3.);
  290. *
  291. pcrack2 = (cont pla2) elem compris pf1 pf2;
  292. crack1 = (diff pla1 pla2);
  293. *
  294. **********************************************************
  295. *** Modele
  296. *
  297. Mod2a = vol2a MODE mecanique elastique isotrope ;
  298. Mod2b = vol2b MODE mecanique elastique isotrope;
  299. *
  300. **********************************************************
  301. Ey1 = 200000. ;
  302. nu1 = 0.3 ;
  303. rho1= 7.8e-3 ;
  304. mat2a = MATE Mod2a 'YOUN' Ey1 'NU' nu1 'RHO' rho1 ;
  305. mat2b = MATE Mod2b 'YOUN' Ey1 'NU' nu1 'RHO' rho1 ;
  306. **********************************************************
  307. *CL
  308. pc3 = (enve vol2) poin proch (3. 6. 0.);
  309. pc4 = (enve vol2) poin proch (3. 6. -6.);
  310. *
  311. cl11 = bloq (pc3 et pc4) 'UZ' ;
  312. cl12 = bloq liab2 'UX' 'UY';
  313. cl13 = bloq liah2 'UX' 'UY';
  314. clto2 = cl11 et cl12 et cl13;
  315. *
  316. nb1 = nbno pla2;
  317. pla3 = chan 'POI1' pla2;
  318. ii = 1;
  319. pp0 = pla3 poin ii;
  320. plb3 = plb2 poin proch pp0;
  321. repeter bouc1 (nb1 - 1);
  322. ii = ii + 1 ;
  323. pp0 = pla3 poin ii;
  324. plb3 = plb3 et (plb2 poin proch pp0);
  325. fin bouc1;
  326. *
  327. cond1 = 'RELA' 'UX' pla3 - 'UX' plb3;
  328. cond2 = 'RELA' 'UY' pla3 - 'UY' plb3;
  329. cond3 = 'RELA' 'UZ' pla3 - 'UZ' plb3;
  330. **********************************************************
  331. *CHARGEMENT
  332. fo11 = FORC 'FZ' liah2 (-3.E3);
  333. fo12 = FORC 'FZ' liab2 (3.E3);
  334. *
  335. *trac (vect (fo11 et fo12) 'FX' 'FY' 'FZ' 1.e-2) ((enve vol2)
  336. *et (liah2 coul vert) et (liab2 coul rose));
  337. *opti donn 5;
  338. *
  339. xx2 = prog 0. 7.5e-4 10.;
  340. yy2 = prog 0. 1. 1.;
  341. evocha2 = evol manu 'temps' xx2 'charge' yy2;
  342. char2 = char meca (fo11 et fo12) evocha2;
  343. *
  344. **********************************************************
  345. *** MATRICES ***
  346. **********************************************************
  347. *** de masse ***
  348. mas2a = mass mod2a mat2a;
  349. mas2b = mass mod2b mat2b ;
  350. **********************************************************
  351.  
  352. *=======================================================*
  353. * *
  354. * RESOLUTION *
  355. * *
  356. *=======================================================*
  357. dt1 = 2.5e-4;
  358. xx0 = prog 0. pas dt1 'NPAS' 4;
  359. *
  360. tab2 = table;
  361. tab2.dynamique = vrai;
  362. tab2.modele = mod2a et mod2b;
  363. tab2.caracteristiques = mat2a et mat2b;
  364. tab2.blocages_mecaniques = clto2 et cond1 et cond2 et cond3;
  365. tab2.chargement = char2;
  366. tab2.temps_calcules = xx0;
  367. pasapas tab2;
  368.  
  369. *post traitement
  370. ind0 = index (tab2.deplacements);
  371. nnn2 = dime ind0;
  372. *
  373. lwela2 = prog;
  374. lwcin2 = prog;
  375. ltemp2 = prog;
  376. ev2 = table;
  377. def2 = table;
  378. *
  379. ii = 0;
  380. repeter bouc2 nnn2;
  381. ii = ii + 1;
  382. pass = ind0.ii;
  383. sig2 = tab2.contraintes.pass;
  384. ut2 = tab2.deplacements.pass;
  385. vt2 = tab2.vitesses.pass;
  386. tt2 = tab2.temps.pass;
  387. etot2 = EPSI ut2 (mod2a et mod2b);
  388. *
  389. ev2.ii = evol rose chpo ut2 'UZ' liah2;
  390. def2.ii = defo ut2 vol2 30. rose;
  391.  
  392. ltemp2 = ltemp2 et (prog tt2);
  393. *...elastique
  394. Wela2 = 'INTG' (0.5 *(ENER sig2 etot2 (mod2a et mod2b)))
  395. (mod2a et mod2b) ;
  396. lwela2 = lwela2 et (prog wela2);
  397. *
  398. *...cinetique
  399. Wcin2 = 0.5 * (XTMX vt2 (mas2a 'ET' mas2b)) ;
  400. lwcin2 = lwcin2 et (prog wcin2);
  401. *
  402. fin bouc2;
  403. **********************************************************
  404. errmax = 0.;
  405. lerr1= prog 0. 4.28804e-08 3.46393e-7 2.36123e-6 1.3508e-5 6.27995e-5;
  406.  
  407. *comparaison
  408. ii = 0;
  409. repeter bouc3 nnn2;
  410. ii = ii + 1 ;
  411.  
  412. *comparaison XFEM et STANDARD 3D
  413. val0 = extr lerr1 ii;
  414. *
  415. *trac cach (def1.ii et def2.ii);
  416. *
  417. *dess (ev1.ii et ev2.ii);
  418. *
  419. ord1 = extr ev1.ii 'ORDO';
  420. ord2 = extr ev2.ii 'ORDO';
  421. si (abs(extr ord2 3) < 1.e-8);
  422. val1 = val0;
  423. message 'ord1 ord2 val1' (extr ord1 3) (extr ord2 3) val1;
  424. sinon;
  425. val1 = ((extr ord1 3) -(extr ord2 3))/(extr ord2 3);
  426. message 'ord1 ord2 val1' (extr ord1 3) (extr ord2 3) val1;
  427. message ' valeurs a comparer' val1 val0;
  428. finsi;
  429. errmax = maxi(prog errmax val1);
  430. *
  431. fin bouc3;
  432. *
  433. evn1 = evol manu 'temps' ltemp1 'ener' lwela1;
  434. evn2 = evol manu 'temps' ltemp1 'ener' lwcin1;
  435. evn3 = evol rose manu 'temps' ltemp2 'ener' lwela2;
  436. evn4 = evol rose manu 'temps' ltemp2 'ener' lwcin2;
  437. *dess (evn1 et evn3);
  438. *dess (evn2 et evn4);
  439.  
  440. val0 = 3.37828e-4;
  441. ord1 = extr evn1 'ORDO';
  442. ord2 = extr evn3 'ORDO';
  443. val1 = abs (((extr ord1 5) -(extr ord2 5))/(extr ord2 5));
  444. message 'ord1 ord2 val1' (extr ord1 5) (extr ord2 5) val1;
  445. message ' valeurs a comparer' val1 val0;
  446. errmax = maxi(prog errmax val1);
  447. *
  448. val0 = 9.78566e-4;
  449. ord1 = extr evn2 'ORDO';
  450. ord2 = extr evn4 'ORDO';
  451. val1 = ((extr ord1 5) -(extr ord2 5))/(extr ord2 5);
  452. message 'ord1 ord2 val1' (extr ord1 5) (extr ord2 5) val1;
  453. message ' valeurs a comparer' val1 val0;
  454. errmax = maxi(prog errmax val1);
  455.  
  456. si(errmax >EG 0.001);
  457. ERRE 5;
  458. finsi;
  459. *
  460. opti donn 5;
  461.  
  462. fin;
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  

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