Télécharger rupt20.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : rupt20.dgibi
  2. opti epsi lineaire ;
  3. ************************************************************************
  4. ************************************************************************
  5.  
  6. *******************************************************
  7. * Test rupt20.dgibi: Jeux de données *
  8. * --------------------------------- *
  9. *******************************************************
  10. * CAS TEST DU 15/12/15 PROVENANCE : TEST
  11.  
  12. *Cas test de validation pour le calcul de J sous plusieurs chargement
  13. *avec les procedures g_theta.procedur et g_calcul.procedur
  14. *
  15. *- chargement en traction
  16. *- chargement avec pression sur levres
  17. *- chargement thermique
  18.  
  19. *Calcul en dimension 3 avec des elements CUB8 sur un maillage complet
  20. *non symetrique
  21.  
  22. opti dime 3 elem cub8 echo 1 ;
  23. ************************
  24. *Données paramètriques :
  25. ************************
  26. * a : profondeur de la fissure *
  27. * t : epaisseur du tube *
  28. * ri, re : rayon interne/externe *
  29. * h : hauteur du tube *
  30.  
  31. h = 1. ;
  32. t = 60.e-3 ;
  33. a = t / 5. ;
  34. ri = t * 5.;
  35. re = ri + t;
  36.  
  37. *NB D'ELEMENT ET LE DEGREE DE REVOLUTION
  38. n1 = 1;
  39. deg1 = 0.5 ;
  40.  
  41. *POINTS POUR L'AXE DE REVOLUTION
  42. p0 = (0. 0. 0.);
  43. py = (0. 1. 0.);
  44.  
  45. *COORDONNEE DE LA POINTE DE LA FISSURE
  46. pf = (a + ri) 0. 0.;
  47.  
  48. *NOMBRE D'ELEMENTS AUTOUR DE LA POINTE DE LA FISSURE (1 et 2 COUT)
  49. nfiss = 10;
  50.  
  51. *TAILLE D'UN ELEMENT DE LA 1ERE ET 2EME COUTURE*
  52. tel = 200e-6 ; tel2 = 400e-6 ;
  53. *Facteur d'agrandissement de la taille du derafinement
  54. ttel2 = 4.*tel2 ;
  55.  
  56. *LONGUEUR DE LA 1ERE ET 2EME COUTURE*
  57. lc1 = nfiss * tel ; lc2 = tel2 * nfiss;
  58.  
  59. *NIVEAU DE CHARGEMENT
  60. p0T = -400. ; p0P = 400. ; dt0 = 300.;
  61.  
  62. *=============================================================
  63. **************************************************************
  64. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  65. * DEBUT DU MAILLAGE
  66. *=============================================================
  67. **************************************************************
  68. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  69.  
  70.  
  71. **************************************************************
  72. ********************** 1ERE COUTURE ************************
  73. ************ (Autour de la pointe de la fissure) ***********
  74. **************************************************************
  75.  
  76. p1cbd = pf plus (lc1 0. 0.) ;
  77. p1chd = pf plus (lc1 lc1 0.) ;
  78. pf1 = pf plus (0. lc1 0.) ;
  79. p1chg = pf1 moin (lc1 0. 0.) ;
  80. p1cbg = pf moin (lc1 0. 0.) ;
  81.  
  82. d1ch = droi (2*nfiss) p1chg p1chd;
  83. d1cg = droi (nfiss) p1cbg p1chg ;
  84. d1cd = droi (nfiss) p1cbd p1chd ;
  85. d1cbg = droi (nfiss) p1cbg pf ;
  86. d1cbd = droi (nfiss) pf p1cbd ;
  87.  
  88. cout1 = regl nfiss d1ch (d1cbg et d1cbd) ;
  89. *cout1 = coul jaun cout1 ;
  90.  
  91. **************************************************************
  92. ********************** 2EME COUTURE ************************
  93. ************ (Autour de la pointe de la fissure) ***********
  94. ************** Partie au-dessus de la fissure **************
  95. **************************************************************
  96.  
  97. p2cbd = pf plus (lc2 0. 0.) ;
  98. p2chd = pf plus (lc2 lc2 0.) ;
  99. pf2 = pf plus (0. lc2 0.) ;
  100. p2chg = pf2 moin (lc2 0. 0.) ;
  101. p2cbg = pf moin (lc2 0. 0.) ;
  102.  
  103. d2ch = droi (2*nfiss) p2chg p2chd;
  104. d2cg = droi (nfiss) p2cbg p2chg ;
  105. d2cd = droi (nfiss) p2cbd p2chd ;
  106.  
  107. cout2 = regl nfiss (d1cd et d1ch et d1cg) (d2cd et d2ch et d2cg) ;
  108. *cout2 = coul vert cout2 ;
  109.  
  110. cout1et2 = cout1 et cout2;
  111.  
  112. **************************************************************
  113. ***************** DERAFINEMENT DES COUTURES ****************
  114. ************** Partie au-dessus de la fissure **************
  115. **************************************************************
  116.  
  117. *------------------( DERAF A 4 ELEMENT )-------------------
  118.  
  119. pid1 = p2chg moin (0. ttel2 0.) ;
  120. pid2 = pid1 plus (0. tel2 0.) ;
  121. pid3 = pid2 plus (0. tel2 0.) ;
  122. pid4 = pid3 plus (0. tel2 0.) ;
  123. pid5 = pid4 plus (0. tel2 0.) ;
  124. pid6 = pid2 moin (tel2 0. 0.) ;
  125. pid7 = pid3 moin (tel2 0. 0.) ;
  126. pid8 = pid4 moin (tel2 0. 0.) ;
  127. pid9 = pid1 moin (ttel2 0. 0.) ;
  128. pid10 = pid3 moin (ttel2 0. 0.) ;
  129. pid11 = pid5 moin (ttel2 0. 0.) ;
  130.  
  131. did1 = droi 1 pid1 pid2 ;
  132. did2 = droi 1 pid2 pid3 ;
  133. did3 = droi 1 pid3 pid4 ;
  134. did4 = droi 1 pid4 pid5 ;
  135. did5 = droi 1 pid9 pid6 ;
  136. did6 = droi 1 pid6 pid7 ;
  137. did7 = droi 1 pid7 pid8 ;
  138. did8 = droi 1 pid11 pid8 ;
  139. did9 = droi 1 pid10 pid7 ;
  140.  
  141. si1 = (regl 1 did1 did5) et (regl 1 did2 did6) et
  142. (regl 1 did3 did7) et (regl 1 did4 (inve did8)) et
  143. (regl 1 did8 did9) et (regl 1 did9 did5) ;
  144. elim si1 1.e-5 ;
  145.  
  146. *trac si1;
  147.  
  148. *------------------( DERAF A 3 ELEMENT )-------------------
  149.  
  150. pad1 = pf moin (lc2 0. 0.) ;
  151. pad2 = pad1 plus (0. tel2 0.) ;
  152. pad3 = pad2 plus (0. tel2 0.) ;
  153. pad4 = pad3 plus (0. tel2 0.) ;
  154. pad5 = pad2 moin (tel2 0. 0.) ;
  155. pad6 = pad3 moin (tel2 0. 0.) ;
  156. pad7 = pad1 moin (ttel2 0. 0.) ;
  157. pad8 = pad4 moin (ttel2 0. 0.) ;
  158.  
  159. dad1 = droi 1 pad1 pad2 ;
  160. dad2 = droi 1 pad2 pad3 ;
  161. dad3 = droi 1 pad3 pad4 ;
  162. dad4 = droi 1 pad7 pad5 ;
  163. dad5 = droi 1 pad5 pad6 ;
  164. dad6 = droi 1 pad6 pad8 ;
  165.  
  166. sa1 = (regl 1 dad1 dad4) et (regl 1 dad2 dad5) et
  167. (regl 1 dad3 dad6) et (regl 1 dad4 (inve dad6));
  168. saa1 = sa1 ;
  169. repe i0 1 ;
  170. ssa1 = sa1 plus ( 0. (3.*&i0*tel2) 0.) ;
  171. fin i0 ;
  172. sa1 = sa1 et ssa1 ;
  173. elim sa1 1.e-5 ;
  174.  
  175. *trac sa1;
  176.  
  177. *---------------------- PARTIE GAUCHE -----------------------
  178. sig = sa1 et si1 ; elim sig 1.e-5 ;
  179.  
  180. *---------------------- PARTIE DROITE -----------------------
  181. sid = sig syme droi ((coor 1 pf) 0. 0.) ((coor 1 pf) lc2 0.) ;
  182. elim sid 1.e-5 ;
  183.  
  184. *---------------------- PARTIE HAUTE -----------------------
  185.  
  186. *lignes diagonales pour la symetrie
  187. p_diagod = p2chd plus (lc1 lc1 0.);
  188. p_diago = p2chg moin (lc1 0. 0.);
  189. p_diagog = p_diago plus (0. lc1 0.);
  190.  
  191. d_diagog = droi 1 p1chg p_diagog;
  192. d_diagod = droi 1 p1chd p_diagod;
  193.  
  194. sihg = sig syme droi p1chg p_diagog ;
  195. elim sihg 1.e-5 ;
  196. sihd = sid syme droi p1chd p_diagod ;
  197. elim sihd 1.e-5 ;
  198.  
  199. sih = sihd et sihg ; elim sih 1.e-5 ;
  200.  
  201. *trac sih;
  202.  
  203. *---------------------- PARTIE COIN -----------------------
  204.  
  205. dg = droi 1 pid11 p2chg;
  206. dcg = dg tran 1 (0. ttel2 0.);
  207. dcd = dcg syme droi ((coor 1 pf) 0. 0.) ((coor 1 pf) lc2 0.) ;
  208. sic = dcd et dcg; elim sic 1.e-5;
  209. cout3 = sig et sid et sih et sic ; elim cout3 1.e-5 ;
  210.  
  211. couttot = cout1 et cout2 et cout3;
  212. elim d1cbd couttot 1.e-5 ;
  213.  
  214. **************************************************************
  215. ********************* RESTE DU MAILLAGE ********************
  216. ************** Partie au-dessus de la fissure **************
  217. **************************************************************
  218.  
  219. *Partie de gauche
  220. *----------------
  221. pt1 = mini (coor 1 couttot) 0. 0.;
  222. pt2 = (mini (coor 1 couttot)) (maxi (coor 2 couttot)) 0.;
  223. ptpartg = cout3 poin droit pt1 pt2 1.e-5 ;
  224. d_partg = (cont couttot) elem appuye strictement ptpartg ;
  225. pri = ri 0. 0.;
  226. pg = d_partg tran
  227. (((coor 1 pri)-(mini(coor 1 couttot))) 0. 0.) dini 1.6e-3 dfin 3.2e-3;
  228.  
  229. *Partie de droite
  230. *----------------
  231. pt3 = maxi (coor 1 couttot) 0. 0.;
  232. pt4 = (maxi (coor 1 couttot)) (maxi (coor 2 couttot)) 0.;
  233. ptpartd = couttot poin droit pt3 pt4 1.e-5 ;
  234. dpartd = (cont couttot) elem appuye strictement ptpartd ;
  235. pre = re 0. 0.;
  236. pd = dpartd tran
  237. (((coor 1 pre)-(maxi(coor 1 couttot))) 0. 0.) dini 1.6e-3 dfin 3.2e-3;
  238. bascout = pg et pd et couttot ; elim bascout 1.e-5 ;
  239.  
  240. *Partie du haut
  241. *--------------
  242. p5 = (mini (coor 1 bascout)) (maxi (coor 2 bascout)) 0.;
  243. p6 = (maxi (coor 1 bascout)) (maxi (coor 2 bascout)) 0.;
  244. pt_parth = bascout poin droit p5 p6 1.e-5 ;
  245. dpartd = (cont bascout) elem appuye strictement pt_parth ;
  246. ph = dpartd tran (0. ((h/2.) - (maxi(coor 2 bascout))) 0.)
  247. dini 1.6e-3 dfin t;
  248.  
  249. *Lèvre et Ligament sur partie haute
  250. *----------------------------------
  251. struhau = ph et bascout ;
  252. pri = struhau poin proc (ri 0. 0.) ;
  253. pre = struhau poin proc (re 0. 0.) ;
  254. prh = struhau poin proc (re (0.5 * h) 0.) ;
  255. ll1 = struhau poin droit pri pre 1e-7 ;
  256. ll2 = (cont struhau) elem appuye strictement ll1 ;
  257. llvsup = ll2 elem comp pri pf ;
  258. llig = ll2 elem comp pre pf ;
  259.  
  260. *Structure haute
  261. *---------------
  262. vhau = struhau volu n1 rota deg1 p0 py ;
  263. levsup = llvsup rota n1 deg1 p0 py ;
  264. slig = llig rota n1 deg1 p0 py ;
  265. elim (vhau et slig et levsup)1.e-7 ;
  266.  
  267. *Structure basse
  268. *---------------
  269. pri2 = pri plus (0. 0. 1.) ;
  270. vbas = vhau syme plan pri pri2 pre ;
  271. levinf = levsup plus (0. 0. 0.) ;
  272. levinf = INVE levinf ;
  273. elim (vbas et levinf et slig) 1.e-7 ;
  274.  
  275. ******************************************************************
  276.  
  277. *bord haut
  278. zm = maxi (coor 2 vhau) ;
  279. sab5 = POIN vhau PLAN (0. zm 0.) (0. zm 1.)(1. zm 0.) 1e-6;
  280. sursup = ELEM (ENVE vhau) APPU STRIC sab5;
  281. sursup = INVE sursup;
  282.  
  283. *Surface bas
  284. zm = mini (coor 2 vbas) ;
  285. sab6 = POIN vbas PLAN (0. zm 0.) (0. zm 1.)(1. zm 0.) 1e-6;
  286. surbas = ELEM (ENVE vbas) APPU STRIC sab6;
  287.  
  288. *Trois points sur la surface f2
  289. zm = mini (vhau coor 3) ;
  290. PC = vhau poin proc (ri 0. zm) ;
  291. PD = vhau poin proc (re 0. zm) ;
  292. PE = vhau poin proc (ri (h / 2.) zm) ;
  293.  
  294. *TRAC (VBAS et VHAU);
  295.  
  296. *=============================================================
  297. **************************************************************
  298. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  299. * FIN DU MAILLAGE
  300. *=============================================================
  301. **************************************************************
  302. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  303.  
  304.  
  305. *=============================================================
  306. **************************************************************
  307. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  308. * PARTIE CALCULS
  309. *=============================================================
  310. **************************************************************
  311. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  312.  
  313. * PROPRIETE MATERIAUX A 300°C
  314. E0 = 185e3; nu0 = 0.3; alfa0 = 13.08e-6;
  315.  
  316. v0 = vbas et vhau;
  317. FRONT = (cont levsup) elem appuye strictement slig ;
  318.  
  319. mo0 = mode v0 mecanique elastique isotrope ;
  320. ma0 = mate mo0 YOUN E0 nu nu0 alph alfa0 ;
  321. rg0 = rigi mo0 ma0 ;
  322.  
  323. *CONDITIONS AUX LIMITES
  324. *Blocages
  325. cl1 = bloq uy surbas ;
  326. cl2 = rela ense uy sursup ;
  327. cl3 = v0 symt depl pri pre prh 1.e-6 ;
  328. cl4 = v0 symt depl PC PD PE 1.e-6 ;
  329. cl0 = cl1 et cl2 et cl3 et cl4 ;
  330.  
  331. *Traction uniaxiale (via un modele de pression)
  332. moph = MODE sursup 'CHARGEMENT' 'PRESSION' 'CONS' 'HAUT' ;
  333.  
  334. *Pression sur les levres (via un modele de pression)
  335. mopl = MODE (levsup ET levinf) 'CHARGEMENT' 'PRESSION' 'CONS' 'LEVRES' ;
  336.  
  337. *Elevation de temperature
  338. chr = ((v0 coor 1) - ri) / (re - ri) ;
  339. cht0 = nomc 'T' (dt0 * chr);
  340.  
  341.  
  342. *****************************************
  343. ********* SOLUTIONS ANALYTIQUES *********
  344. *****************************************
  345. *Fonctions d'influence
  346. i0 = 1.211 ;
  347. i1 = 0.718 ;
  348.  
  349. *Contraintes imposées pour le gradient de temperature
  350. sig0 = ((E0*alfa0*dt0)/(1-nu0)) * (ri/(3*t)) *
  351. ((2*(re**2))/(ri*(re+ri)) - 1);
  352. sig1 = -1. * ((E0*alfa0*dt0)/(1-nu0));
  353.  
  354. * J analytiques
  355. JT = (1-(nu0**2)) * ((i0*(-1.*p0T)*((pi*a)**(1./2.)))**2) / E0;
  356. JP = (1-(nu0**2)) * ((i0*(-1.*p0P)*((pi*a)**(1./2.)))**2) / E0;
  357. JTH = (1-(nu0**2)) *
  358. ((((i0*sig0)+ (i1*sig1*(a/t)))*((pi*a)**(1./2.))) **2)/E0;
  359.  
  360.  
  361. ****************************************************
  362. * CALCUL ELASTIQUE AVEC RESO - CALCUL DE J ELASTIQUE
  363. ****************************************************
  364.  
  365. * Construction des second membres
  366. maph = pres moph 'PRES' p0T ;
  367. f0T = BSIG moph maph ;
  368. mapl = pres mopl 'PRES' p0P ;
  369. f0P = BSIG mopl mapl ;
  370.  
  371. sgth0 = THET mo0 ma0 cht0 ;
  372. f0TH = BSIG mo0 sgth0 ;
  373.  
  374. * RESOLUTION ELASTIQUE DES 3 PROBLEMES
  375. utestT utestP utestTH = RESO (rg0 ET cl0) f0T f0P f0TH;
  376.  
  377. *PROCEDURE G_THETA
  378. *cas 1 : traction seule
  379. tabJel = table ;
  380. tabJel . 'MODELE' = mo0 ET moph ;
  381. tabJel . 'CARACTERISTIQUES' = ma0 ;
  382. tabJel . 'PRESSION' = maph ;
  383. tabJel . 'BLOCAGES_MECANIQUES' = cl0 ;
  384. tabJel . 'SOLUTION_RESO' = utestT ;
  385. tabJel . 'OBJECTIF' = MOT 'J' ;
  386. tabJel . 'LEVRE_SUPERIEURE' = levsup ;
  387. tabJel . 'LEVRE_INFERIEURE' = levinf ;
  388. tabJel . 'FRONT_FISSURE' = FRONT ;
  389. tabJel . 'COUCHE' = 5 ;
  390. g_theta tabJel ;
  391. JelT1 = tabJel.resultats.global ;
  392.  
  393. *cas 2 : pression sur les levres
  394. tabJel . 'MODELE' = mo0 ET mopl ;
  395. tabJel . 'CARACTERISTIQUES' = ma0 ;
  396. tabJel . 'PRESSION' = mapl ;
  397. tabJel . 'SOLUTION_RESO' = utestP ;
  398. g_theta tabJel ;
  399. JelP1 = tabJel.resultats.global ;
  400.  
  401. *cas 3 : gradient de temperature
  402. tabJel . 'MODELE' = mo0 ;
  403. tabJel . 'CARACTERISTIQUES' = ma0 ;
  404. tabJel . 'TEMPERATURES' = cht0 ;
  405. tabJel . 'SOLUTION_RESO' = utestTH ;
  406. g_theta tabJel ;
  407. JelTH1 = tabJel.resultats.global ;
  408.  
  409. *Erreurs sur J entre la solution analytique et le MEF
  410. errT1 = ((JelT1-JT)/JT)*100.;
  411. errP1 = ((JelP1-JP)/JP)*100.;
  412. errTH1 = ((JelTH1-JTH)/JTH)*100.;
  413.  
  414.  
  415.  
  416.  
  417.  
  418. *******************************************************
  419. * CALCUL ELASTIQUE AVEC PASAPAS - CALCUL DE J ELASTIQUE
  420. *******************************************************
  421.  
  422. * Chargements de pression (obligatoires si modele de pression)
  423. evph = EVOL 'MANU' 'TEMP' (PROG 0. 1. 2. 3.)
  424. 'PRES' (PROG 0. 1. 0. 0.) ;
  425. chaph = CHAR 'PRES' maph evph ;
  426.  
  427. evpl = EVOL 'MANU' 'TEMP' (PROG 0. 1. 2. 3.)
  428. 'PRES' (PROG 0. 0. 1. 0.) ;
  429. chapl = CHAR 'PRES' mapl evpl ;
  430.  
  431. * Chargement thermique
  432. chath = CHAR 'T' cht0 (EVOL 'MANU' (PROG 0. 1. 2. 3.)
  433. (PROG 0. 0. 0. 1.)) ;
  434.  
  435. *RESOLUTION AVEC PASAPAS DES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS)
  436. *AU PAS 1 : Traction seule
  437. *AU PAS 2 : Pression sur les levres
  438. *AU PAS 3 : Gradient de temperature
  439. tabT = TABL ;
  440. tabT . 'MODELE' = mo0 ET moph ET mopl ;
  441. tabT . 'CARACTERISTIQUES' = ma0 ;
  442. tabT . 'BLOCAGES_MECANIQUES' = cl0 ;
  443. tabT . 'CHARGEMENT' = chaph ET chapl ET chath ;
  444. tabT . 'TEMPS_CALCULES' = PROG 1. 2. 3. ;
  445. PASAPAS tabT ;
  446.  
  447. *PROCEDURE G_THETA POUR LES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS)
  448. *ATTENTION, IL FAUT RETIRER LE CHARGERMENT MECA DE PRESSION SUR LES
  449. *LEVRES ET UTILISER LE CHARGEMENT PLEV
  450. *PROCEDURE G_THETA POUR LES 3 PROBLEMES (UN A CHAQUE PAS DE TEMPS)
  451. tabJel = TABL ;
  452. tabJel . 'SOLUTION_PASAPAS' = tabT ;
  453. tabJel . 'OBJECTIF' = MOT 'J' ;
  454. tabJel . 'LEVRE_SUPERIEURE' = levsup ;
  455. tabJel . 'LEVRE_INFERIEURE' = levinf ;
  456. tabJel . 'FRONT_FISSURE' = FRONT ;
  457. tabJel . 'COUCHE' = 5 ;
  458. g_theta tabJel ;
  459. JelT2 = tabJel.resultats. 1 . global ;
  460. JelP2 = tabJel.resultats. 2 . global ;
  461. JelTH2 = tabJel.resultats. 3 . global ;
  462. *Erreurs sur J : solution analytique VS calcul PASAPAS + G_THETA
  463. errT2 = ((JelT2-JT)/JT)*100.;
  464. errP2 = ((JelP2-JP)/JP)*100.;
  465. errTH2 = ((JelTH2-JTH)/JTH)*100.;
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472. ****************************************
  473. * AFFICHAGE DES RESULTATS ET DES ERREURS
  474. ****************************************
  475. SAUT 5 'LIGNE' ;
  476. mess 'Solution Theorique : ' JT JP JTH ;
  477. mess 'Solution MEF (RESO) : ' JelT1 JelP1 JelTH1 ;
  478. mess 'Erreur en % : ' errT1 errP1 errTH1 ;
  479. mess 'Solution MEF (PASAPAS) : ' JelT2 JelP2 JelTH2 ;
  480. mess 'Erreur en % : ' errT2 errP2 errTH2 ;
  481.  
  482. * Test sur les erreurs
  483. errT = MAXI 'ABS' (PROG errT1 errT2) ;
  484. si ((abs errT) > 0.5) ;
  485. erre 'Erreur sur le calcul de JelT' ;
  486. fins ;
  487. errP = MAXI 'ABS' (PROG errP1 errP2) ;
  488. si ((abs errP) > 0.5) ;
  489. erre 'Erreur sur le calcul de JelP' ;
  490. fins ;
  491. errTH = MAXI 'ABS' (PROG errTH1 errTH2) ;
  492. si ((abs errTH) > 0.2) ;
  493. erre 'Erreur sur le calcul de JelTH' ;
  494. fins ;
  495.  
  496. FIN ;
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  

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