Télécharger la plaquette Cast3M


fr en

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

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