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

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