Télécharger pq1.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : pq1.dgibi
  2. *---------------------------------------------------------
  3. * Cas test de non régression sur le controle des maillages
  4. *---------------------------------------------------------
  5. * Pressurisation d'une enceinte type Phébus
  6. *
  7. * Le maillage correspond à une enceinte cylindrique
  8. * d'environ 10 m3 avec un mur en contact avec la
  9. * paroi verticale (10 cm d'acier)
  10. * Tout le volume est initialement a 1bar et 40oC
  11. * et la température du mur est mise à 60oC
  12. * au début du calcul. On calcule la pressurisation
  13. * de cette enceinte sur 50 secondes en injectant un débit
  14. * de 50g/s de vapeur à 150oC. Ce test (un peu long)
  15. * verifie le demarrage de la condensation a 20"
  16. * verifie l'evolution moyenne de la temprature gaz
  17. * verifie la pression max a 50"
  18. * verifie la vitesse max a 50" (Convection naturelle)
  19. * verifie la masse d'eau a 50"
  20. * Auteurs : E. Studer, Novembre 1999
  21. * revisite Mars 2001
  22. *--------------------------------------------------------
  23. COMPLET = FAUX ;
  24. *COMPLET = VRAI ;
  25. GRAPH = VRAI ;
  26. GRAPH = FAUX ;
  27. 'OPTI' 'TRAC' 'X' ;
  28. *'OPTI' 'TRAC' 'PSC' ;
  29.  
  30. 'SI' COMPLET ;
  31.  
  32. nbit=100;
  33. DT0 = 1. ;
  34. n1 = 1 ; n2 = 4 ; n3 = 4 ;
  35. n4 = 8 ; nn = 2 ;
  36.  
  37. 'SINON' ;
  38.  
  39. nbit= 5 ;
  40. DT0 = 10. ;
  41. n1 = 1 ; n2 = 2 ; n3 = 4 ;
  42. n4 = 4 ; nn = 1 ;
  43.  
  44. 'FINSI' ;
  45.  
  46. *--------------------------------------------------------
  47. * Definition du maillage de l'enceinte cylindrique
  48. *
  49. 'OPTI' 'DIME' 3 'ELEM' 'CU20' ;
  50.  
  51. ri = 1.052 ; sp = 0.10 ; hc = 4.163 ;
  52.  
  53.  
  54.  
  55. epsi = 1.000e-2 ; ;
  56.  
  57. p0 = 0.000 0.000 0.000 ;
  58. px = -1000.000 0.000 0.000 ;
  59. py = 0.000 -1000.000 0.000 ;
  60. pz = 0.000 0.000 1000.000 ;
  61. cd = 0.000 0.000 -20.000 ;
  62. ph0 = 0.000 0.000 hc ;
  63. phx = ri 0.000 hc ;
  64. phy = 0.000 ri hc ;
  65.  
  66. fg1 = 0.25 ;
  67. fg2 = fg1 * (2.0 ** 0.5) / 2. ;
  68.  
  69. p1 = (ri*fg1) 0.000 0.000 ;
  70. p2 = (ri*fg2) (ri*fg2) 0.000 ;
  71. p3 = 0.000 (ri*fg1) 0.000 ;
  72. p4 = ri 0.000 0.000 ;
  73. p5 = 0.000 ri 0.000 ;
  74. p6 = (ri+sp) 0.000 0.000 ;
  75. p7 = 0.000 (ri+sp) 0.000 ;
  76.  
  77. * Hauteur de l'enceinte
  78. h1 = 4.163 ;
  79. * Vecteur de translation
  80. v1 = 0. 0. h1 ;
  81.  
  82. l1 = 'DROI' p0 p1 n1 ;
  83. l2 = 'DROI' p1 p2 n1 ;
  84. l3 = 'DROI' p2 p3 n1 ;
  85. l4 = 'DROI' p3 p0 n1 ;
  86. l5 = 'CERC' p4 p0 p5 (2*n1) ;
  87. l6 = 'CERC' p6 p0 p7 (2*n1) ;
  88.  
  89. basf0= 'DALL' l1 l2 l3 l4 'PLAN' ;
  90. basf1=('REGL' (l2 'ET' l3) l5 n2) ;
  91. l44= cote 2 basf1;
  92. ax4= (inve l4) et l44 ;
  93. l11= cote 4 basf1;
  94. ax1= l11 et (inve l1) ;
  95.  
  96. basf = basf0 'ET' ('REGL' (l2 'ET' l3) l5 n2) ;
  97. 'ELIM' basf epsi ;
  98. basf = basf 'ET' ('SYME' basf 'DROI' p0 p3) ;
  99. ax11 = ('SYME' ax1 'DROI' p0 p3) 'ET' (inve ax1) ;
  100. 'ELIM' basf epsi ;
  101. basf = basf 'ET' ('SYME' basf 'DROI' p0 p1) ;
  102. ax44 = (inve ax4) 'ET' ('SYME' ax4 'DROI' p0 p4) ;
  103. 'ELIM' basf epsi ;
  104.  
  105. basm = 'REGL' l5 l6 n3 ;
  106. basm = basm 'ET' ('SYME' basm 'DROI' p0 p3) ;
  107. 'ELIM' basm epsi ;
  108. basm = basm 'ET' ('SYME' basm 'DROI' p0 p1) ;
  109. 'ELIM' basm epsi ;
  110.  
  111. * Creation du volume
  112.  
  113. dx = ri / 2. ;
  114. nz1 = ('ENTIER' ( h1 / dx ))*nn ;
  115.  
  116. bas = basf 'VOLU' nz1 'TRAN' v1 ;
  117. mbas = basm 'VOLU' nz1 'TRAN' v1 ;
  118. mbas = mbas 'COUL' 'ROUG' ;
  119. plan1 = ax11 'TRAN' nz1 v1 ;
  120. plan4 = ax44 'TRAN' nz1 v1 ;
  121. 'ELIM' (bas et plan1 et plan4) epsi ;
  122.  
  123. mt = bas ;
  124. wall = mbas ;
  125. elim (mt et wall) epsi ;
  126.  
  127. * Localisation d'une brèche éventuelle au bas de l'enceinte
  128.  
  129. pjg = 'POIN' basf 'PROC' (0.000 0.000 0.000) ;
  130. jg = ('ELEM' basf 'APPUIE' 'LARGEMENT' pjg) 'COUL' 'VERT' ;
  131.  
  132. *--------------------------------------------------------------------
  133. * Fin de la définition du maillage
  134. *--------------------------------------------------------------------
  135.  
  136.  
  137. *--------------------------------------------------------------------
  138. * Début de l'initialisation de la procédure ENCEINTE : table RXT
  139. *--------------------------------------------------------------------
  140.  
  141. rxt = 'TABLE' ;
  142. rxt . 'VERSION' = 'V0' ;
  143.  
  144. *-- Nom du volume fluide
  145. rxt . 'vtf' = mt ;
  146. rxt . 'epsi' = epsi ;
  147.  
  148. *-- Definition des murs de l'enceinte : ici un seul mur
  149. *-- en ACIER dont on traitera la thermique dans l'épaisseur
  150. *-- et que l'on initialise a 40oC
  151.  
  152. *-- On definit d'abord le matériau ACIER avec sa conductivite
  153. *-- thermique LAMBDA (W/m/K) et le produit ro*Cp (J/m3/K)
  154. rxt . 'THERMP' = VRAI ;
  155. rxt . 'vtp' = wall ;
  156. rxt . 'LAMBDA' = 15. ;
  157. rxt . 'ROCP' = 3.9E6 ;
  158. rxt . 'Tp0' = 60. ;
  159. rxt . 'ECHAN' = 10. ;
  160.  
  161. *-- On definit un point interne au maillage pour imposer la valeur de
  162. *-- la pression
  163. rxt . 'pi' = (0.0 0.0 0.5) ;
  164.  
  165. *-- On indique que le calcul comporte de la vapeur d'eau
  166. rxt . 'VAPEUR' = VRAI ;
  167. rxt . 'Yvap0' = 0.0023 ;
  168.  
  169. *-- Conditions initiales dans l'enceinte de test
  170. rxt . 'TF0' = 40.0 ;
  171. rxt . 'PT0' = 1.00000e5 ;
  172.  
  173. *-- Definition du scenario à la breche
  174.  
  175. rxt . 'Breches' = 'TABLE' ;
  176. rxt . 'Breches' . 'A' = 'TABLE' ;
  177. rxt . 'Breches' . 'A' . 'Maillage' = jg ;
  178. rxt . 'Breches' . 'A' . 'diru' = (0. 0. 1.) ;
  179. rxt . 'Breches' . 'A' . 'scenario' = 'TABLE' ;
  180. rxt . 'Breches' . 'A' . 'scenario' . 't' = PROG 0.0 1000. ;
  181. rxt . 'Breches' . 'A' . 'scenario' . 'qair' = PROG 0.0 0.0 ;
  182. rxt . 'Breches' . 'A' . 'scenario' . 'qeau' = PROG 0.05 0.05 ;
  183. rxt . 'Breches' . 'A' . 'scenario' . 'tinj' = PROG 150. 150. ;
  184.  
  185. *-- On impose la viscosite turbulente (m2/s)
  186. rxt . 'MODTURB' = 'NUTURB' ;
  187. rxt . 'NUT' = 1.e-2 ;
  188.  
  189. *-- On impose le pas de temps (s)
  190. rxt . 'DT0' = DT0 ;
  191.  
  192. *-- On active le recalcul automatique du préconditionnement
  193. *-- toutes les 5 itérations
  194. rxt . 'FRPREC' = 5 ;
  195. rxt . 'DETMAT' = VRAI ;
  196. *rxt.'TRTF' = FAUX ;
  197.  
  198.  
  199. Si VRAI;
  200. *-- On lance le calcul sur 20 itérations d'une seconde
  201. rxt . 'GRAPH' = GRAPH ;
  202.  
  203. EXECRXT 2 rxt ;
  204. EXECRXT (nbit - 2) rxt ;
  205.  
  206. list rxt.'TIC'.'Tfm' ;
  207. list rxt.'TIC'.'PT' ;
  208. list rxt.'TIC'.'Qc' ;
  209. list rxt.'TIC'.'LMAXU';
  210.  
  211. ltfm = 'PROG'
  212. 40.000 65.434 73.712 81.158 86.809 90.721 ;
  213. lPT = 'PROG'
  214. 1.00000E+05 1.06009E+05 1.21199E+05 1.29592E+05 1.36461E+05 1.43987E+05 ;
  215. Lqc = 'PROG'
  216. 0.0000 0.0000 0.0000 4.06133E-04 2.65247E-03 4.88852E-03 ;
  217. Lmaxu = 'PROG'
  218. 0.0000 0.81320 2.0813 2.6901 2.2350 2.5277 ;
  219.  
  220.  
  221. tic=rxt.'TIC' ;
  222. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  223. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  224. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  225. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  226.  
  227. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  228.  
  229. ierr = 0 ;
  230. Si (ERtf '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  231. Si (ERPT '>' 1.e-3) ; ierr = ierr '+' 1 ; Finsi ;
  232. Si (ERQc '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  233. Si (ERum '>' 1.e-2) ; ierr = ierr '+' 1 ; Finsi ;
  234.  
  235. ******************************************************************
  236.  
  237. rxt.'vtf' = (chan 'MACRO' mt) ;
  238. rxt.'vtp' = (chan 'MACRO' wall) ;
  239. rxt . 'Breches' . 'A' . 'Maillage' = (chan 'MACRO' jg) ;
  240.  
  241. rxt = enlev rxt 'GEO';
  242. rxt = enlev rxt 'TBT';
  243. rxt = enlev rxt 'TIC';
  244.  
  245. EXECRXT 2 rxt ;
  246. EXECRXT (nbit - 2) rxt ;
  247.  
  248. list rxt.'TIC'.'Tfm' ;
  249. list rxt.'TIC'.'PT' ;
  250. list rxt.'TIC'.'Qc' ;
  251. list rxt.'TIC'.'LMAXU';
  252.  
  253. ltfm = 'PROG'
  254. 40.000 65.434 73.716 81.163 86.907 90.942 ;
  255. lPT = 'PROG'
  256. 1.00000E+05 1.06006E+05 1.21166E+05 1.29578E+05 1.36479E+05 1.44098E+05 ;
  257. Lqc = 'PROG'
  258. 0.0000 0.0000 0.0000 3.11548E-04 2.41832E-03 4.44576E-03 ;
  259. Lmaxu = 'PROG'
  260. 0.0000 0.81320 2.0840 2.6866 2.2313 2.5160 ;
  261.  
  262.  
  263. tic=rxt.'TIC' ;
  264. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  265. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  266. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  267. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  268.  
  269. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  270.  
  271. Si (ERtf '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  272. Si (ERPT '>' 1.e-3) ; ierr = ierr '+' 1 ; Finsi ;
  273. Si (ERQc '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  274. Si (ERum '>' 1.e-2) ; ierr = ierr '+' 1 ; Finsi ;
  275.  
  276.  
  277. ******************************************************************
  278.  
  279. rxt.'vtf' = (chan 'QUAF' mt) ;
  280. rxt.'vtp' = (chan 'QUAF' wall) ;
  281. rxt . 'Breches' . 'A' . 'Maillage' = (chan 'QUAF' jg) ;
  282. rxt = enlev rxt 'GEO';
  283. rxt = enlev rxt 'TBT';
  284. rxt = enlev rxt 'TIC';
  285.  
  286. EXECRXT 2 rxt ;
  287. EXECRXT (nbit - 2) rxt ;
  288.  
  289. list rxt.'TIC'.'Tfm' ;
  290. list rxt.'TIC'.'PT' ;
  291. list rxt.'TIC'.'Qc' ;
  292. list rxt.'TIC'.'LMAXU';
  293.  
  294. ltfm = 'PROG'
  295. 40.000 65.434 73.712 81.158 86.809 90.721 ;
  296. lPT = 'PROG'
  297. 1.00000E+05 1.06009E+05 1.21199E+05 1.29592E+05 1.36461E+05 1.43987E+05;
  298. Lqc = 'PROG'
  299. 0.0000 0.0000 0.0000 4.06133E-04 2.65247E-03 4.88852E-03;
  300. Lmaxu = 'PROG'
  301. 0.0000 0.81320 2.0813 2.6901 2.2350 2.5277 ;
  302.  
  303. tic=rxt.'TIC' ;
  304. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  305. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  306. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  307. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  308.  
  309. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  310.  
  311. Si (ERtf '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  312. Si (ERPT '>' 1.e-3) ; ierr = ierr '+' 1 ; Finsi ;
  313. Si (ERQc '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  314. Si (ERum '>' 1.e-2) ; ierr = ierr '+' 1 ; Finsi ;
  315.  
  316. finsi;
  317.  
  318. ******************************************************************
  319. rxt.'vtf' = (chan 'LINE' mt) ;
  320. rxt.'vtp' = (chan 'LINE' wall) ;
  321. rxt . 'Breches' . 'A' . 'Maillage' = (chan 'LINE' jg) ;
  322. rxt = enlev rxt 'GEO';
  323. rxt = enlev rxt 'TBT';
  324. rxt = enlev rxt 'TIC';
  325.  
  326. EXECRXT 2 rxt ;
  327. EXECRXT (nbit - 2) rxt ;
  328.  
  329. list rxt.'TIC'.'Tfm' ;
  330. list rxt.'TIC'.'PT' ;
  331. list rxt.'TIC'.'Qc' ;
  332. list rxt.'TIC'.'LMAXU';
  333.  
  334. ltfm = 'PROG'
  335. 40.000 67.363 75.733 83.257 88.699 92.375 ;
  336. lPT = 'PROG'
  337. 1.00000E+05 1.06452E+05 1.22785E+05 1.31710E+05 1.39047E+05 1.46818E+05;
  338. Lqc = 'PROG'
  339. 0.0000 0.0000 0.0000 1.04818E-03 3.74210E-03 6.39247E-03 ;
  340. Lmaxu = 'PROG'
  341. 0.0000 0.81320 2.0604 2.5919 2.1446 2.4526 ;
  342.  
  343. tic=rxt.'TIC' ;
  344. ERtf=somm( abs (ltfm - tic.'Tfm') )/ 80. ;
  345. ERPT=somm( abs (lPT - tic.'PT' ) ) /1.e5 ;
  346. ERQc=somm( abs (lqc - tic.'Qc' ) ) ;
  347. ERum=somm( abs (Lmaxu - tic.'LMAXU' ) )/2. ;
  348.  
  349. Mess 'ERtf=' ERtf 'ERPT=' ERPT 'ERQc=' ERQc 'ERum=' ERum ;
  350.  
  351. Si (ERtf '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  352. Si (ERPT '>' 1.e-3) ; ierr = ierr '+' 1 ; Finsi ;
  353. Si (ERQc '>' 1.e-4) ; ierr = ierr '+' 1 ; Finsi ;
  354. Si (ERum '>' 1.e-2) ; ierr = ierr '+' 1 ; Finsi ;
  355.  
  356.  
  357.  
  358. Si GRAPH ;
  359.  
  360. $vtf=rxt.'GEO'.'$vtf' ;
  361. mt =doma $vtf maillage ;
  362. Mpl1=chan 'QUAF' plan1 ;
  363. Mpl4=chan 'QUAF' plan4 ;
  364. ELIM (mt et Mpl1 et Mpl4) epsi ;
  365. $mpl1= mode Mpl1 'NAVIER_STOKES' MACRO ;
  366. $mpl4= mode Mpl4 'NAVIER_STOKES' MACRO ;
  367. plan1= doma $mpl1 maillage ;
  368. plan4= doma $mpl4 maillage ;
  369. plan=plan1 et plan4 ;
  370.  
  371. paroif = rxt.'GEO'.'paroif';
  372. rho=rxt.'TIC'.'RHO';
  373. rvp=rxt.'TIC'.'RVP';
  374. tf=rxt.'TIC'.'TF';
  375. un=rxt.'TIC'.'UN' ;
  376. un1=redu un plan ;
  377. ung= vect un1 1. ux uy uz jaune ;
  378. trace ung plan ;
  379. opti isov suli ;
  380. trace rho plan 'TITRE'' Rho' ;
  381. trace rvp plan 'TITRE'' Rvp' ;
  382. trace tf plan 'TITRE'' Tf ' ;
  383.  
  384. trace rho paroif 'TITRE'' Rho' ;
  385. trace rvp paroif 'TITRE'' Rvp' ;
  386. trace Tf paroif 'TITRE'' Tf ' ;
  387. Fcond = tic.'Fcondw';
  388. trace Fcond paroif 'TITRE' ' Flux de condensation Kg / m**2 ' ;
  389.  
  390. axe = p0 d nz1 (p0 plus v1) ;
  391. axe = chan axe 'QUAF' ;
  392. elim (axe et mt) epsi ;
  393.  
  394. evr= evol 'CHPO' rho axe ;
  395. dess evr 'TITRE' ' Rho axe ';
  396.  
  397. evt = evol 'CHPO' Tf axe ;
  398. dess evt 'TITRE' ' Température axe ';
  399.  
  400. evvp= evol 'CHPO' rvp axe ;
  401. dess evvp 'TITRE' ' rvp axe ';
  402.  
  403.  
  404.  
  405. 'FINSI' ;
  406.  
  407.  
  408. 'SI' (IERR '>' 0) ;
  409. 'MESS' 'Il y a des problemes !!!' ;
  410. 'ERRE' 5 ;
  411. 'SINO';
  412. 'MESS' 'Tout s est bien passe!' ;
  413. 'FINS' ;
  414.  
  415.  
  416.  
  417. 'FIN' ;
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  

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