Télécharger rupt21.dgibi

Retour à la liste

Numérotation des lignes :

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

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