Télécharger tubedeto3d1.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : tubedeto3d1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ****************************************************************
  5. ***** PROPAGATION D UNE DETONATION DANS UN TUBE *****
  6. ***** MODELE CREBCOM *****
  7. ***** A. BECCANTINI, SFME/LTMF, 15.01.02 *****
  8. ****************************************************************
  9. *
  10. * The file is divided into 5 parts
  11. *
  12. * 1) mesh
  13. * 2) initial conditions and gas properties
  14. * 3) parameters used in the computation
  15. * 4) the main part (where the computation is realised)
  16. * 5) the post-treatment
  17. *
  18.  
  19. 'OPTION' 'ECHO' 0 'DIME' 3
  20. 'TRAC' 'X' ;
  21.  
  22. ****************************************************************
  23. ***** 1D mesh *****
  24. ****************************************************************
  25. *
  26. *
  27. * A6 ---- A5 ------------------- A4 ----------------------- A8
  28. * | | | |
  29. * | AR | ZONE1 | ZONE2 |
  30. * | | | |
  31. * A1 ---- A2 ------------------- A3 ----------------------- A7
  32. * | | | |
  33. * | L1 | L2 | L3 |
  34. * |<---->|<--------------------->|<----------------------->|
  35. *
  36. * AR = activation region
  37. * ZONE1
  38. * ZONE2
  39. *
  40.  
  41. RAF = 1 ;
  42.  
  43. NL2SL1 = 49 ;
  44. NL3SL1 = 100 ;
  45. L1 = 0.4 ;
  46.  
  47. L2 = NL2SL1 '*' L1 ;
  48. L3 = NL3SL1 '*' L1 ;
  49.  
  50. * If NL2SL1 = 49 and NL3SL1 = 50, then L2 = 19.6, L3 = 20.0 ;
  51.  
  52. N1 = RAF ;
  53. N2 = NL2SL1 '*' N1 ;
  54. N3 = NL3SL1 '*' N1 ;
  55.  
  56. * For sake of simplicity, we will only consider one mesh in y-direction
  57.  
  58. DX = L1 '/' N1 ;
  59. NY = 4 ;
  60. NZ = NY ;
  61. H1 = NY '*' DX ;
  62.  
  63. A1 = 0.0 0.0 0.0 ;
  64. A2 = L1 0.0 0.0 ;
  65. A3 = (L1 '+' L2) 0.0 0.0 ;
  66. A4 = (L1 '+' L2) H1 0.0 ;
  67. A5 = L1 H1 0.0 ;
  68. A6 = 0.0 H1 0.0 ;
  69. A7 = (L1 '+' L2 '+' L3) 0.0 0.0 ;
  70. A8 = (L1 '+' L2 '+' L3) H1 0.0 ;
  71.  
  72. 'OPTION' 'ELEM' 'SEG2' ;
  73. A1A2 = A1 'DROIT' N1 A2 ;
  74. A2A3 = A2 'DROIT' N2 A3 ;
  75. A3A4 = A3 'DROIT' NY A4 ;
  76. A4A5 = A4 'DROIT' N2 A5 ;
  77. A5A2 = A5 'DROIT' NY A2 ;
  78. A2A5 = 'INVERSE' A5A2 ;
  79. A5A6 = A5 'DROIT' N1 A6 ;
  80. A6A1 = A6 'DROIT' NY A1 ;
  81. A3A7 = A3 'DROIT' N3 A7 ;
  82. A7A8 = A7 'DROIT' NY A8 ;
  83. A8A4 = A8 'DROIT' N3 A4 ;
  84. A4A3 = 'INVERSE' A3A4 ;
  85.  
  86. *
  87. **** 2D mesh
  88. *
  89.  
  90. * DOM1 is the activation region
  91. * DOM2 is zone 1
  92. * DOM3 is zone 2
  93.  
  94. 'OPTION' 'ELEM' 'QUA4' ;
  95. DOM1P = 'TRANSLATION' A2A5 N1 ((-1.0 '*' L1) 0.0 0.0) ;
  96. DOM2P = 'TRANSLATION' A2A5 N2 (L2 0.0 0.0) ;
  97. DOM3P = 'TRANSLATION' A3A4 N3 (L3 0.0 0.0) ;
  98.  
  99. * DOMTOT = total region
  100.  
  101. DOMTOTP = DOM1P 'ET' DOM2P 'ET' DOM3P ;
  102. 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOTP ;
  103. 'ELIMINATION' DOMTOTP (1.0D-3 '/' RAF) (A1A2 'ET' A2A3 'ET'
  104. A3A4 'ET' A4A5 'ET' A5A2 'ET' A5A6 'ET' A6A1 'ET' A3A7
  105. 'ET' A7A8 'ET' A8A4 'ET' A4A3) ;
  106.  
  107. 'OPTION' 'ELEM' 'CUB8' ;
  108. DOM1 = 'VOLUME' DOM1P 'TRANSLATION' NZ (0.0 0.0 H1) ;
  109. DOM2 = 'VOLUME' DOM2P 'TRANSLATION' NZ (0.0 0.0 H1) ;
  110. DOM3 = 'VOLUME' DOM3P 'TRANSLATION' NZ (0.0 0.0 H1) ;
  111.  
  112. DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
  113. 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;
  114.  
  115. *
  116. **** The tables 'DOMAINE'
  117. *
  118.  
  119. 'SI' VRAI ;
  120. * EULER model
  121. $DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
  122. $DOM1 = 'MODELISER' DOM1 'EULER' ;
  123. $DOM2 = 'MODELISER' DOM2 'EULER' ;
  124. $DOM3 = 'MODELISER' DOM3 'EULER' ;
  125.  
  126. TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
  127. TDOM1 = 'DOMA' $DOM1 'VF' ;
  128. TDOM2 = 'DOMA' $DOM2 'VF' ;
  129. TDOM3 = 'DOMA' $DOM3 'VF' ;
  130.  
  131. QDOMTOT = TDOMTOT . 'QUAF' ;
  132. QDOM1 = TDOM1 . 'QUAF' ;
  133. QDOM2 = TDOM2 . 'QUAF' ;
  134. QDOM3 = TDOM3 . 'QUAF' ;
  135.  
  136. 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM1 ;
  137. 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM2 ;
  138. 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM3 ;
  139. 'SINON' ;
  140. * NAVIER_STOKES model
  141. DOM1 = 'CHANGER' 'QUAF' DOM1 ;
  142. DOM2 = 'CHANGER' 'QUAF' DOM2 ;
  143. DOM3 = 'CHANGER' 'QUAF' DOM3 ;
  144.  
  145. DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
  146. 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;
  147.  
  148. $DOMTOT = 'MODELISER' DOMTOT 'NAVIER_STOKES' 'LINE' ;
  149. $DOM1 = 'MODELISER' DOM1 'NAVIER_STOKES' 'LINE' ;
  150. $DOM2 = 'MODELISER' DOM2 'NAVIER_STOKES' 'LINE' ;
  151. $DOM3 = 'MODELISER' DOM3 'NAVIER_STOKES' 'LINE' ;
  152.  
  153. DOMTOT = 'DOMA' $DOMTOT 'MAILLAGE' ;
  154. DOM1 = 'DOMA' $DOM1 'MAILLAGE' ;
  155. DOM2 = 'DOMA' $DOM2 'MAILLAGE' ;
  156. DOM3 = 'DOMA' $DOM3 'MAILLAGE' ;
  157.  
  158. TDOMTOT = 'DOMA' $DOMTOT 'TABLE' ;
  159. TDOMTOT . 'PRECONDI' = 1 ;
  160. 'FINSI' ;
  161.  
  162.  
  163. *
  164. **** MOD1 = model (created just to display the CHAMELEMs)
  165. *
  166.  
  167. MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
  168.  
  169. *
  170. **** Line for graphics
  171. *
  172.  
  173. 'OPTION' 'ELEM' 'SEG2' ;
  174. P1 = (DX '/' 2) (DX '/' 2) (DX '/' 2) ;
  175. P2 = ((L1 '+' L2 '+' L3) '-' (DX '/' 2)) (DX '/' 2) (DX '/' 2) ;
  176. LIGEVO = (P1 'DROIT' (N1 '+' N2 '+' N3 '-' 1) P2) 'COULEUR' 'ROUG' ;
  177. 'ELIMINATION' LIGEVO ('DOMA' $DOMTOT 'CENTRE') (0.001 '/' RAF) ;
  178.  
  179. * 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
  180. * 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;
  181.  
  182. ****************************************************************
  183. ****************************************************************
  184. ***** Initial conditions and gas properties. *****
  185. ****************************************************************
  186. ****************************************************************
  187.  
  188. *
  189. **** 0 from a numerical point of view
  190. *
  191.  
  192. zero = 1.0d-9 ;
  193.  
  194. *
  195. **** (Non-Uniform) initial conditions.(DOM1 et DOM2)
  196. * We compute the conservative variables.
  197. *
  198. *
  199. * Contact discontinuity between DOM2 and DOM3 (same pressure,
  200. * different temperature)
  201. *
  202.  
  203. PINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
  204. ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 99700 'NATU' 'DISCRET') 'ET'
  205. ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 99700
  206. 'NATU' 'DISCRET') ;
  207. TINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
  208. ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 285. 'NATU' 'DISCRET') 'ET'
  209. ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 400
  210. 'NATU' 'DISCRET') ;
  211. WINI = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 3 'UX' 0.0 'UY' 0.0
  212. 'UZ' 0.0 ;
  213.  
  214. * We suppose that the combustion is complete everywhere, i.e. the
  215. * progress variable is 1 after the detonation transition.
  216.  
  217. CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0 ;
  218.  
  219. *
  220. **** From initial conditions on molar fractions, we duduce
  221. * the data for the virtual mixture of burned/unburned
  222. * gases
  223.  
  224. *
  225. **** Molar fractions (before combustion)
  226. *
  227.  
  228. XH21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
  229. ('DOMA' $DOM2 'CENTRE')) 1 'H2' 0.3 'NATU' 'DISCRET') 'ET'
  230. ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'H2' 0.2
  231. 'NATU' 'DISCRET') ;
  232. XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' zero
  233. 'NATURE' 'DISCRET' ;
  234. XO21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
  235. ('DOMA' $DOM2 'CENTRE')) 1 'O2' 0.147 'NATU' 'DISCRET') 'ET'
  236. ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'O2' 0.2
  237. 'NATU' 'DISCRET') ;
  238. XN = XH21 'ET' XH2O1 'ET' XO21 ;
  239. CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0D0) ;
  240. XN21 = CHUN '-' ('PSCAL' XN CHUN ('MOTS' 'H2' 'O2' 'H2O')
  241. ('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
  242. XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
  243. XN = XN 'ET' XN21 ;
  244.  
  245. *
  246. ***** Gas properties
  247. *
  248.  
  249. PGAZ = 'TABLE' ;
  250.  
  251. * Polynomial degree of specific heats
  252.  
  253. PGAZ . 'NORD' = 4 ;
  254.  
  255. * Species explicitly treated in the Euler Equations
  256.  
  257. PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
  258.  
  259. * Species non explicitly treated
  260.  
  261. PGAZ . 'ESPNEULE' = 'N2 ';
  262.  
  263. * Single gases properties
  264.  
  265. PGAZ . 'H2 ' = 'TABLE' ;
  266. PGAZ . 'H2O ' = 'TABLE' ;
  267. PGAZ . 'N2 ' = 'TABLE' ;
  268. PGAZ . 'O2 ' = 'TABLE' ;
  269.  
  270. * R (J/Kg/K)
  271.  
  272. mH2 = 2. '*' 1.00797E-3 ;
  273. mo2 = 2. '*' 15.9994E-3 ;
  274. mH2O = mh2 '+' (0.5 '*' mo2) ;
  275. mN2 = 2 '*' 14.0067E-3 ;
  276. RGAS = 8.31441 ;
  277.  
  278. PGAZ . 'H2 ' . 'R' = RGAS '/' mh2 ;
  279. PGAZ . 'H2O ' . 'R' = RGAS '/' mh2o ;
  280. PGAZ . 'N2 ' . 'R' = RGAS '/' mn2 ;
  281. PGAZ . 'O2 ' . 'R' = RGAS '/' mo2 ;
  282.  
  283. * Polynomials regressions coefficients
  284.  
  285. PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
  286. -2.37281455E-07 1.84701105E-11 ;
  287. PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
  288. -1.82753232E-08 2.44485692E-12 ;
  289. PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
  290. 8.78233606E-09 -3.05514485E-13 ;
  291. PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
  292. 2.33636971E-08 -1.53304905E-12;
  293.  
  294. * Formation enthalpies at 0K (J/Kg))
  295. *
  296. * h_i(0K) = h_i(T0) '-' \int_0^{T0} cp_i(x) dx
  297. * = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
  298. *
  299.  
  300. PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
  301. PGAZ . 'H2O ' . 'H0K' = -1.395D7 ;
  302. PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
  303. PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
  304.  
  305. * Passive scalars names
  306.  
  307. PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ';
  308.  
  309. *
  310.  
  311. * 'OPTION' 'SAUVER' 'gas_properties.sauv' ;
  312. * 'SAUVER' PGAZ ;
  313.  
  314. *
  315. **** We compute the mass fractions
  316. *
  317.  
  318. MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' Mh2
  319. 'NATURE' 'DISCRET') 'ET'
  320. ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' Mo2
  321. 'NATURE' 'DISCRET') 'ET'
  322. ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' Mh2o
  323. 'NATURE' 'DISCRET') 'ET'
  324. ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'N2' Mn2
  325. 'NATURE' 'DISCRET') ;
  326.  
  327. MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  328. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  329.  
  330. YH21 = (XH21 '*' MH2) '/' MTOT ;
  331. YO21 = (XO21 '*' MO2) '/' MTOT ;
  332. YN21 = (XN21 '*' MN2) '/' MTOT ;
  333. YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;
  334.  
  335. *
  336. **** Burned gas
  337. * Here x refers to the mole producted per 1 mole of unburned gas
  338. * N.B. Number of moles changes after the combustion!
  339. *
  340.  
  341. *
  342. * H2MAG = 1 in the region where H2 is the majority species.
  343. * Conversely H2MAG = 0
  344. *
  345.  
  346. H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
  347. ('EXCO' 'O2' XN) ;
  348. O2MAG = CHUN '-' H2MAG ;
  349.  
  350. XO22a = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' zero ;
  351. XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
  352.  
  353. XO22a = XO22a ;
  354. XH22a = XH22a ;
  355.  
  356. XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' zero ;
  357. XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (xh21 '-' xh22b))) ;
  358.  
  359. XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
  360. XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;
  361.  
  362. XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
  363. XN22 = XN21 ;
  364.  
  365. *
  366. **** mtot = weight of a mole of unburned gas
  367. * mtot2 = weight of a mole of unburned gas after the combustion
  368. * Of course mtot2 = mtot
  369. *
  370.  
  371. MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
  372. ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  373. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  374.  
  375. mtotcel = 'MAXIMUM' MTOT ;
  376.  
  377. 'SI' ( ('MAXIMUM' (mtot2 '-' mtot) 'ABS') '>' (zero '*' mtotcel)) ;
  378. 'MESSAGE' 'Probleme dan le calcul du gaz brule' ;
  379. 'ERREUR' 21 ;
  380. 'FINSI' ;
  381.  
  382. *
  383. **** xtot2 = mole of burnet gas per mole of unburned gas
  384. *
  385.  
  386. xtot2 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
  387. CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  388. ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
  389.  
  390. *
  391. **** x = real molar fractions
  392. *
  393.  
  394. xh22 = xh22 '/' xtot2 ;
  395. xo22 = xo22 '/' xtot2 ;
  396. xh2o2 = xh2o2 '/' xtot2 ;
  397. xn22 = xn22 '/' xtot2 ;
  398.  
  399.  
  400. xtot3 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
  401. CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  402. ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
  403.  
  404. 'SI' ( ('MAXIMUM' (xtot3 '-' 1.0D0) 'ABS' ) '>'
  405. zero) ;
  406. 'MESSAGE' 'Probleme dan le calcul du gaz brule' ;
  407. 'ERREUR' 21 ;
  408. 'FINSI' ;
  409.  
  410. *
  411. **** mtot2 = weight of 1 mole of burned gas
  412. *
  413.  
  414. MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
  415. ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  416. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  417.  
  418. *
  419. **** Mass fractions of burned gas after the complete combustion
  420. *
  421.  
  422. YH22 = (XH22 '*' MH2) '/' MTOT2 ;
  423. YO22 = (XO22 '*' MO2) '/' MTOT2 ;
  424. YN22 = (XN22 '*' MN2) '/' MTOT2 ;
  425. YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
  426.  
  427. *
  428. **** The conservative variables
  429. *
  430.  
  431. RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
  432. (('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
  433. (('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
  434. (('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
  435. ;
  436.  
  437. RN = PINI '/' (RINI '*' TINI) ;
  438.  
  439. * Internal energy
  440.  
  441. eini = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 0.0 ;
  442. Tfois = Tini ;
  443.  
  444. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  445.  
  446. acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  447. &BL1)) '+'
  448. (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  449. &BL1)) '+'
  450. (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  451. &BL1)) '+'
  452. (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  453. &BL1)) ;
  454.  
  455. eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
  456.  
  457. Tfois = Tfois * tini ;
  458.  
  459. 'FIN' BL1 ;
  460.  
  461. RETN = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
  462. ('MOUS' 'UX' 'UY' 'UZ') ('MOTS' 'UX' 'UY' 'UZ')))) ;
  463. GN = RN '*' WINI ;
  464. YN = YH21 '+' YO21 '+' YH2O1 ;
  465. RYN = RN '*' YN ;
  466.  
  467. *
  468. **** K0/DX = chemical source term
  469. *
  470. * i.e. dcsi/dt = K0 '/' DX . {criterium function}
  471. * where
  472. * {criterium function} = 1.0 or 0.0
  473. *
  474. * d(rho etot)/dt = -rho (dcsi/dt) (h_f '-' h_i)
  475. *
  476. * h_f reference enthalpy (0K) of the burned gas
  477. * h_i reference enthalpy (0K) of the unburned gas
  478. *
  479. **** Mass fraction of H2 before and after the combustion
  480. *
  481. * CSIMAX = 1.0 -> complete combustion
  482. *
  483.  
  484. YININ = yh21 ;
  485. YFINN = YININ '+' ((yh22 '-' yh21) '*' csimax) ;
  486. K0N = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 500 ;
  487.  
  488. *
  489. **** Passive scalars
  490. *
  491.  
  492. RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
  493. ('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
  494. ('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
  495.  
  496. 'SI' VRAI ;
  497. * Test
  498. VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
  499. PCEL = 'MAXIMUM' PINI ;
  500. TCEL = 'MAXIMUM' TINI ;
  501. 'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
  502. 'MESSAGE' 'Probleme calcul energie totale' ;
  503. 'ERREUR' 21 ;
  504. 'FINSI' ;
  505. 'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
  506. 'MESSAGE' 'Probleme calcul energie totale' ;
  507. 'ERREUR' 21 ;
  508. 'FINSI' ;
  509. 'FINSI' ;
  510.  
  511. ****************************************************************
  512. ****************************************************************
  513. ***** Parameters for the computations ******
  514. ****************************************************************
  515. ****************************************************************
  516.  
  517. * Upwind scheme
  518. * METO = 'VLH' ;
  519. METO = 'SS' ;
  520.  
  521. * Iterations
  522. * Final time
  523. * Safety Factor fot the time step
  524. * Second order reconstruction?
  525.  
  526. NITER = 10000 ;
  527. TFINAL = 5.0D-3 ;
  528. SAFFAC = 0.7 ;
  529. LOGSO = VRAI ;
  530.  
  531. *
  532. * EPS_CC linked to CREBCOM criterion
  533. *
  534.  
  535. EPS_CC = 0.8 ;
  536.  
  537. * 'OPTION' 'SAUVER' 'parameters.sauv' ;
  538. * 'SAUVER' METO NITER TFINAL SAFFAC EPS_CC LOGSO ;
  539.  
  540. ****************************************************************
  541. ****************************************************************
  542. ***** The computation ******
  543. ****************************************************************
  544. ****************************************************************
  545.  
  546. * Mesh dimension (in 3D!)
  547.  
  548. DELTAXN = ('DOMA' $DOMTOT 'VOLUME') '**' (1. '/' 3.) ;
  549. YN = RYN '/' RN ;
  550.  
  551. *
  552. **** We impose we have combustion in the ignition region (DOM1)
  553. * We can use the operator 'FLAM', with epsilon=0 and
  554. * DELTAT >> DELTATC
  555. * where DELTATC is the characteristic time step linked to
  556. * the source term.
  557. * Since DELTAXN is constant, DELTATC does not change during the
  558. * computation.
  559. *
  560.  
  561. SN = RSN '/' RN ;
  562. K0N = 'EXCO' 'K0' SN ;
  563. YININ = 'EXCO' 'H2IN' SN 'H2' ;
  564. YFINN = 'EXCO' 'H2FI' SN 'H2' ;
  565.  
  566. DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
  567. LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
  568. LCOEF = 'PROG' 1.0 0.5 -1.0 ;
  569.  
  570. DELTARE DELTARY = 'FLAM' 'CREBCOM2' $DOM1 PGAZ LMOT1 LCOEF
  571. ('REDU' RN ('DOMA' $DOM1 'CENTRE'))
  572. ('REDU' YN ('DOMA' $DOM1 'CENTRE'))
  573. ('REDU' YININ ('DOMA' $DOM1 'CENTRE'))
  574. ('REDU' YFINN ('DOMA' $DOM1 'CENTRE'))
  575. ('REDU' K0N ('DOMA' $DOM1 'CENTRE'))
  576. ('REDU' DELTAXN ('DOMA' $DOM1 'CENTRE')) 0.0 (1D4 '*' DELTATC) 0.3 ;
  577.  
  578. RYN = RYN '+' DELTARY ;
  579. YN = RYN '/' RN ;
  580. RETN = RETN '+' DELTARE ;
  581.  
  582.  
  583. RN0 = 'COPIER' RN ;
  584. GN0 = 'COPIER' GN ;
  585. RETN0 = 'COPIER' RETN ;
  586. RYN0 = 'COPIER' RYN ;
  587. RSN0 = 'COPIER' RSN ;
  588.  
  589. *
  590. **** Parameter for the time loop
  591. *
  592.  
  593. * Names of the residuum components
  594.  
  595. LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RUZ' 'RETN' 'H2' 'O2' 'H2O' 'H2IN'
  596. 'H2FI' 'K0' ;
  597.  
  598. *
  599. **** Geometrical coefficient to compute gradients
  600. *
  601.  
  602. GRADR CACCA COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  603. ('MOTS' 'SCAL') RN ;
  604.  
  605. GRADV CACCA COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  606. ('MOTS' 'UX' 'UY' 'UZ') GN ;
  607.  
  608. TPS = 0.0 ;
  609.  
  610. *
  611. **** The chemistry time step
  612. *
  613. *
  614.  
  615. DT_CHEM = SAFFAC '*' DELTATC ;
  616.  
  617. *
  618. **** Temporal loop
  619. *
  620.  
  621. 'MESSAGE' ;
  622. 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  623. 'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
  624. 'MESSAGE' ;
  625.  
  626. VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
  627.  
  628. TNM1 = 'COPIER' TN ;
  629.  
  630. 'TEMPS' 'ZERO' ;
  631. 'REPETER' BL1 NITER ;
  632.  
  633. *
  634. **** Primitive variables
  635. *
  636.  
  637. VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN
  638. TNM1 ;
  639.  
  640. TNM1 = 'COPIER' TN ;
  641.  
  642. 'SI' LOGSO ;
  643.  
  644. GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  645. ('MOTS' 'SCAL') RN 'GRADGEO' COEFSCAL ;
  646.  
  647. GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  648. ('MOTS' 'SCAL') PN 'GRADGEO' COEFSCAL ;
  649.  
  650. GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  651. ('MOTS' 'UX' 'UY' 'UZ') VN 'GRADGEO' COEFVECT ;
  652.  
  653. GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  654. ('MOTS' 'H2' 'O2' 'H2O') YN 'GRADGEO' COEFSCAL ;
  655.  
  656. GRADS ALS = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  657. ('MOTS' 'H2IN' 'H2FI' 'K0') SN 'GRADGEO' COEFSCAL ;
  658.  
  659. ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 2 1
  660. $DOMTOT PGAZ
  661. RN GRADR ALR
  662. VN GRADV ALV
  663. PN GRADP ALP
  664. YN GRADY ALY
  665. SN GRADS ALS ;
  666.  
  667. 'SINON' ;
  668.  
  669. ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 1 1
  670. $DOMTOT PGAZ
  671. RN
  672. VN
  673. PN
  674. YN
  675. SN ;
  676.  
  677. 'FINSI' ;
  678.  
  679. RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
  680. $DOMTOT PGAZ LISTINCO ROF VITF PF YF SF ;
  681.  
  682. DT_CON = SAFFAC '*' DELTAT ;
  683. *
  684. **** The time step linked to tps
  685. *
  686.  
  687. DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  688.  
  689. *
  690. **** Total time step
  691. *
  692.  
  693. DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM) ;
  694.  
  695. *
  696. **** Increment of the variables (convection)
  697. *
  698.  
  699. RESIDU = DTMIN '*' RESIDU ;
  700.  
  701. DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  702. DGN = 'EXCO' ('MOTS' 'RUX' 'RUY' 'RUZ') RESIDU ('MOTS' 'UX' 'UY'
  703. 'UZ') ;
  704. DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  705. DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
  706. ('MOTS' 'H2' 'O2' 'H2O') ;
  707. DRSN = 'EXCO' ('MOTS' 'H2IN' 'H2FI' 'K0') RESIDU
  708. ('MOTS' 'H2IN' 'H2FI' 'K0') ;
  709.  
  710.  
  711. TPS = TPS '+' DTMIN ;
  712. RN = RN '+' DRN ;
  713. GN = GN '+' DGN ;
  714. RETN = RETN '+' DRETN ;
  715. RYN = RYN '+' DRYN ;
  716. RSN = RSN '+' DRSN ;
  717. YN = RYN '/' RN ;
  718. SN = RSN '/' RN ;
  719.  
  720. *
  721. **** Increment of the variables (chemical source)
  722. *
  723.  
  724. K0N = 'EXCO' 'K0' SN 'SCAL' ;
  725. YININ = 'EXCO' 'H2IN' SN 'H2' ;
  726. YFINN = 'EXCO' 'H2FI' SN 'H2' ;
  727.  
  728. DRETN DRYN = 'FLAM' 'CREBCOM2' $DOMTOT PGAZ LMOT1 LCOEF RN YN
  729. YININ YFINN K0N DELTAXN EPS_CC DTMIN 0.3 ;
  730.  
  731. RYN = RYN '+' DRYN ;
  732. RETN = RETN '+' DRETN ;
  733.  
  734.  
  735. 'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  736. 'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
  737. 'FINSI' ;
  738.  
  739. 'SI' (TPS > TFINAL) ;
  740. 'QUITTER' BL1 ;
  741. 'FINSI' ;
  742.  
  743. 'FIN' BL1 ;
  744. 'TEMPS' ;
  745.  
  746. * 'OPTION' 'SAUVER' ('CHAINE' 'result.sauv_' RAF 'tps_' TPS ) ;
  747. * 'SAUVER' ;
  748.  
  749. ****************************************************************
  750. ****************************************************************
  751. ***** The post-treatment ******
  752. ****************************************************************
  753. ****************************************************************
  754.  
  755. * 'OPTI' 'REST' 'result.sauv_1tps_5' ;
  756. * 'REST' ;
  757.  
  758. GRAPH = VRAI ;
  759. GRAPH = FAUX ;
  760.  
  761. *
  762. **** The mesh
  763. *
  764.  
  765. 'SI' GRAPH ;
  766. 'TRACER' DOMTOT 'TITR' ('CHAINE' 'Maillage: '
  767. ('NBEL' DOMTOT) ' elts');
  768. 'FINSI' ;
  769.  
  770. *
  771. **** Initial conditions
  772. *
  773.  
  774. MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
  775.  
  776. 'SI' GRAPH ;
  777.  
  778. VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP'
  779. PGAZ RN0 GN0 RETN0 RYN0 RSN0 ;
  780.  
  781. CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN0 ;
  782. CHM_VN = 'KCHA' $DOMTOT 'CHAM' VN ;
  783. CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
  784. CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
  785. CHM_YN = 'KCHA' $DOMTOT 'CHAM' YN ;
  786. CHM_SN = 'KCHA' $DOMTOT 'CHAM' SN ;
  787.  
  788. 'TRAC' CHM_RN MOD1
  789. 'TITR' ('CHAINE' 'rho at t=' 0.0) ;
  790. 'TRAC' CHM_VN MOD1
  791. 'TITR' ('CHAINE' 'v at t= ' 0.0) ;
  792. 'TRAC' CHM_PN MOD1
  793. 'TITR' ('CHAINE' 'p at t= ' 0.0) ;
  794. 'TRAC' CHM_TN MOD1
  795. 'TITR' ('CHAINE' 'T at t= ' 0.0) ;
  796. 'TRAC' CHM_YN MOD1
  797. 'TITR' ('CHAINE' 'y at t= ' 0.0) ;
  798. 'TRAC' CHM_SN MOD1
  799. 'TITR' ('CHAINE' 's at t= ' 0.0) ;
  800.  
  801. 'FINSI' ;
  802.  
  803. *
  804. **** Results
  805. *
  806.  
  807. VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
  808.  
  809. 'SI' GRAPH ;
  810.  
  811. CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
  812. CHM_VN = 'KCHA' $DOMTOT 'CHAM' VN ;
  813. CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
  814. CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
  815. CHM_YN = 'KCHA' $DOMTOT 'CHAM' YN ;
  816. CHM_SN = 'KCHA' $DOMTOT 'CHAM' SN ;
  817.  
  818. 'TRAC' CHM_RN MOD1
  819. 'TITR' ('CHAINE' 'rho at t=' TPS) ;
  820. 'TRAC' CHM_VN MOD1
  821. 'TITR' ('CHAINE' 'v at t= ' TPS) ;
  822. 'TRAC' CHM_PN MOD1
  823. 'TITR' ('CHAINE' 'p at t= ' TPS) ;
  824. 'TRAC' CHM_TN MOD1
  825. 'TITR' ('CHAINE' 'T at t= ' TPS) ;
  826. 'TRAC' CHM_YN MOD1
  827. 'TITR' ('CHAINE' 'y at t= ' TPS) ;
  828. 'TRAC' CHM_SN MOD1
  829. 'TITR' ('CHAINE' 's at t= ' TPS) ;
  830.  
  831. 'FINSI' ;
  832.  
  833. *
  834. *** Evolution Objects
  835. *
  836.  
  837. SS = 'COORDONNEE' 1 LIGEVO;
  838.  
  839. evxx = 'EVOL' 'CHPO' ss LIGEVO ;
  840. lxx = 'EXTRAIRE' evxx 'ORDO' ;
  841.  
  842. x0 = 'MINIMUM' lxx;
  843. x1 = 'MAXIMUM' lxx ;
  844.  
  845. titolo = 'CHAINE' ' at t = ' TPS ' Methode = ' METO
  846. ' SF = ' SAFFAC ;
  847.  
  848. VNN = 'EXCO' 'UX' VN ;
  849. VNT = 'EXCO' 'UY' VN ;
  850.  
  851. * rho
  852.  
  853. evro = 'EVOL' 'CHPO' RN LIGEVO ;
  854. lro = 'EXTRAIRE' evro 'ORDO' ;
  855. evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro ;
  856. tro = 'CHAINE' 'ro ' titolo ;
  857.  
  858. * u
  859.  
  860. evu = 'EVOL' 'CHPO' VNN LIGEVO ;
  861. lu = 'EXTRAIRE' evu 'ORDO' ;
  862. evu = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
  863. tu = 'CHAINE' 'u ' titolo ;
  864.  
  865. * yini
  866.  
  867. YININ = 'EXCO' 'H2IN' SN ;
  868. evyini = 'EVOL' 'CHPO' YININ LIGEVO ;
  869. lyini = 'EXTRAIRE' evyini 'ORDO' ;
  870. evyini = 'EVOL' 'MANU' 'x' lxx 'yini' lyini ;
  871. tyini = 'CHAINE' 'yini ' titolo ;
  872.  
  873. * yfin
  874.  
  875. YFINN = 'EXCO' 'H2FI' SN ;
  876. evyfin = 'EVOL' 'CHPO' YFINN LIGEVO ;
  877. lyfin = 'EXTRAIRE' evyfin 'ORDO' ;
  878. evyfin = 'EVOL' 'MANU' 'x' lxx 'yfin' lyfin ;
  879. tyfin = 'CHAINE' 'yfin ' titolo ;
  880.  
  881. * K0
  882.  
  883. K0N = 'EXCO' 'K0' SN ;
  884. evk0 = 'EVOL' 'CHPO' K0N LIGEVO ;
  885. lk0 = 'EXTRAIRE' evk0 'ORDO' ;
  886. evk0 = 'EVOL' 'MANU' 'x' lxx 'k0' lk0 ;
  887. tk0 = 'CHAINE' 'k0 ' titolo ;
  888.  
  889. * csi
  890.  
  891. CSIN = (('EXCO' 'H2' YN) '-' YININ) '/'
  892. (YFINN '-' YININ) ;
  893. evcsi = 'EVOL' 'CHPO' CSIN LIGEVO ;
  894. lcsi = 'EXTRAIRE' evcsi 'ORDO' ;
  895. evcsi = 'EVOL' 'MANU' 'x' lxx 'csi' lcsi ;
  896. tcsi = 'CHAINE' 'csi ' titolo ;
  897.  
  898. * v
  899.  
  900. evv = 'EVOL' 'CHPO' VNT LIGEVO ;
  901. lv = 'EXTRAIRE' evv 'ORDO' ;
  902. evv = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
  903. tv = 'CHAINE' 'v ' titolo ;
  904.  
  905.  
  906. * p
  907.  
  908. evp = 'EVOL' 'CHPO' PN LIGEVO ;
  909. lp = 'EXTRAIRE' evp 'ORDO' ;
  910. evp = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
  911. tp = 'CHAINE' 'p ' titolo ;
  912.  
  913.  
  914. * T
  915.  
  916. evT = 'EVOL' 'CHPO' TN LIGEVO ;
  917. lT = 'EXTRAIRE' evT 'ORDO' ;
  918. evT = 'EVOL' 'MANU' 'x' lxx 'T' lT ;
  919. tT = 'CHAINE' 'T ' titolo ;
  920.  
  921.  
  922. 'SI' GRAPH ;
  923. 'DESSIN' evro 'TITRE' tro 'XBOR' x0 x1 'MIMA' ;
  924. 'DESSIN' evp 'TITRE' tp 'XBOR' x0 x1 'MIMA' ;
  925. 'DESSIN' evu 'TITRE' tu 'XBOR' x0 x1 'MIMA' ;
  926. 'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1 'MIMA' ;
  927. 'DESSIN' evT 'TITRE' tT 'XBOR' x0 x1 'MIMA' ;
  928. 'DESSIN' evcsi 'TITRE' tcsi 'XBOR' x0 x1 'MIMA' ;
  929. 'DESSIN' evk0 'TITRE' tk0 'XBOR' x0 x1 'MIMA' ;
  930. 'DESSIN' evyini 'TITRE' tyini 'XBOR' x0 x1 'MIMA' ;
  931. 'DESSIN' evyfin 'TITRE' tyfin 'XBOR' x0 x1 'MIMA' ;
  932. 'FINSI' ;
  933.  
  934.  
  935. * Test
  936.  
  937. PVN = 3.2E6 ;
  938. PCJ = 1.1E6 ;
  939. RHOVN = 3.7 ;
  940. RHOCJ = 1.6 ;
  941.  
  942. PMAX = 'MAXIMUM' PN ;
  943. 'SI' ((PMAX '>' PVN) 'OU' (PMAX '<' PCJ)) ;
  944. 'MESSAGE' 'Pression = ???' ;
  945. 'ERREUR' 5 ;
  946. 'FINSI' ;
  947.  
  948. RHOMAX = 'MAXIMUM' RN ;
  949. 'SI' ((RHOMAX '>' RHOVN) 'OU' (RHOMAX '<' RHOCJ)) ;
  950. 'MESSAGE' 'Densite = ???' ;
  951. 'ERREUR' 5 ;
  952. 'FINSI' ;
  953.  
  954. 'FIN' ;
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961.  
  962.  
  963.  
  964.  
  965.  
  966.  
  967.  
  968.  

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