Télécharger comb.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : comb.dgibi
  2. *************************************************************************
  3. *
  4. * NOM : comb.dgibi
  5. *
  6. * DESCRIPTION : We compute the AICC (Adiabatic Isochoric Complete
  7. * Combustion), the AIBCC, (Adiabatic IsoBaric Complete
  8. * Combustion), the CJ (Chapman - Jouguet) and VN (von
  9. * Neumann) states of a ZND detonation for several mixtures
  10. * at NIST-normal conditions involving H2 and air.
  11. * We check that some physical properties are satisfied.
  12. * Comparison of the results obtained using the DETO
  13. * operator and the ones given by the CREBCOM combustion
  14. * model (operators PRIM and FLAM).
  15. *
  16. * LANGAGE : GIBIANE-CAST3M
  17. * AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
  18. * beccantini@cea.fr
  19. * DATE : 28/09/2006
  20. *************************************************************************
  21. *
  22.  
  23. 'OPTION' 'ECHO' 1 'DIME' 2
  24. 'TRAC' 'X' 'ELEM' 'QUA4' ;
  25.  
  26. ****************************************************************
  27. ***** 1D mesh *****
  28. ****************************************************************
  29. *
  30. *
  31. * A4 ---- A3
  32. * | |
  33. * | COM |
  34. * | |
  35. * A1 ---- A2
  36. *
  37.  
  38. A1 = 0.0 0.0 ;
  39. A2 = 1.0 0.0 ;
  40.  
  41. A1A2 = A1 'DROIT' 1 A2 ;
  42. DOM1 = A1A2 'TRANSLATION' 1 (0.0 1.0) ;
  43.  
  44. *
  45. **** The table 'DOMAINE'
  46. *
  47.  
  48. $DOM1 = 'MODELISER' DOM1 'EULER' ;
  49. TDOM1 = 'DOMA' $DOM1 'VF' ;
  50. QDOM1 = TDOM1 . 'QUAF' ;
  51.  
  52.  
  53. *
  54. ***** Initial conditions and gas properties.
  55. *
  56.  
  57. *
  58. **** 0 from a numerical point of view
  59. *
  60.  
  61. zero = 1.0d-9 ;
  62.  
  63. *
  64.  
  65. PINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.013D5
  66. 'NATU' 'DISCRET' ;
  67. TINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 298.15
  68. 'NATU' 'DISCRET' ;
  69. WINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 2 'UX' 0.0 'UY' 0.0
  70. 'NATU' 'DISCRET' ;
  71.  
  72. * We suppose that the combustion is complete everywhere
  73.  
  74. CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0 ;
  75.  
  76. *
  77. ***** Gas properties
  78. *
  79.  
  80. PGAZ = 'TABLE' ;
  81.  
  82. * Polynomial degree of specific heats
  83.  
  84. PGAZ . 'NORD' = 4 ;
  85.  
  86. * Species explicitly treated in the Euler Equations
  87.  
  88. PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
  89.  
  90. * Species non explicitly treated
  91.  
  92. PGAZ . 'ESPNEULE' = 'N2 ';
  93.  
  94. * Single gases properties
  95.  
  96. PGAZ . 'H2 ' = 'TABLE' ;
  97. PGAZ . 'H2O ' = 'TABLE' ;
  98. PGAZ . 'N2 ' = 'TABLE' ;
  99. PGAZ . 'O2 ' = 'TABLE' ;
  100.  
  101. * R (J/Kg/K)
  102.  
  103. mH2 = 2. '*' 1.00797E-3 ;
  104. mo2 = 2. '*' 15.9994E-3 ;
  105. mH2O = mh2 '+' (0.5 '*' mo2) ;
  106. mN2 = 2 '*' 14.0067E-3 ;
  107. RGAS = 8.31441 ;
  108.  
  109. PGAZ . 'H2 ' . 'R' = RGAS '/' mh2 ;
  110. PGAZ . 'H2O ' . 'R' = RGAS '/' mh2o ;
  111. PGAZ . 'N2 ' . 'R' = RGAS '/' mn2 ;
  112. PGAZ . 'O2 ' . 'R' = RGAS '/' mo2 ;
  113.  
  114. * Polynomials regressions coefficients
  115.  
  116. PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
  117. -2.37281455E-07 1.84701105E-11 ;
  118. PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
  119. -1.82753232E-08 2.44485692E-12 ;
  120. PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
  121. 8.78233606E-09 -3.05514485E-13 ;
  122. PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
  123. 2.33636971E-08 -1.53304905E-12;
  124.  
  125. * Formation enthalpies at 0K (J/Kg)
  126. *
  127. * h_i(0K) = h_i(T0) '-' \int_0^{T0} cp_i(x) dx
  128. * = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
  129. *
  130.  
  131. PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
  132. PGAZ . 'H2O ' . 'H0K' = -1.395D7 ;
  133. PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
  134. PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
  135.  
  136. * Passive scalars names
  137.  
  138. PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ';
  139.  
  140. *
  141. **** Initial conditions on molar fractions of hydrogen.
  142. *
  143.  
  144. LISTXH2 = PROG 0.08 0.134 0.188 0.242 0.29581 0.350 ;
  145. LISXH2F = PROG ;
  146. LISXO2F = PROG ;
  147. LISXWAF = PROG ;
  148. LISXN2F = PROG ;
  149. LISYH2F = PROG ;
  150. LISYO2F = PROG ;
  151. LISYWAF = PROG ;
  152. LISYN2F = PROG ;
  153. * Inside
  154. LRINSI = PROG ;
  155. LTINSI = PROG ;
  156. LPINSI = PROG ;
  157. LCINSI = PROG ;
  158. LRGINSI = PROG ;
  159. LGINSI = PROG ;
  160. LGAINSI = PROG ;
  161. * VN conditions
  162. LRVN = PROG ;
  163. LTVN = PROG ;
  164. LPVN = PROG ;
  165. LCVN = PROG ;
  166. LRGVN = PROG ;
  167. LGVN = PROG ;
  168. LGAVN = PROG ;
  169. * AICC
  170. LRAICC = PROG ;
  171. LTAICC = PROG ;
  172. LPAICC = PROG ;
  173. LCAICC = PROG ;
  174. LRGAICC = PROG ;
  175. LGAICC = PROG ;
  176. LGAAICC = PROG ;
  177. * AIBCC
  178. LRAIBC = PROG ;
  179. LTAIBC = PROG ;
  180. LPAIBC = PROG ;
  181. LCAIBC = PROG ;
  182. LRGAIBC = PROG ;
  183. LGAIBC = PROG ;
  184. LGAAIBC = PROG ;
  185. * CJ conditions
  186. LRCJ = PROG ;
  187. LTCJ = PROG ;
  188. LPCJ = PROG ;
  189. LCCJ = PROG ;
  190. LRGCJ = PROG ;
  191. LGCJ = PROG ;
  192. LGACJ = PROG ;
  193. * Released energy per mass unit
  194. LQ = PROG ;
  195. *
  196.  
  197.  
  198. REPE BLCALC (DIME LISTXH2) ;
  199.  
  200. ************************************************************************
  201. ************************************************************************
  202. ***** LOOP to determine the initial conditions, AICC, AIBCC, CJ *****
  203. ***** anv VN states *****
  204. ************************************************************************
  205. ************************************************************************
  206.  
  207. *
  208. **** Molar fractions (before combustion)
  209. *
  210.  
  211. cell = (EXTR LISTXH2 &BLCALC) ;
  212. cell1 = (1. - cell) / (1. + (0.79/0.21)) ;
  213. XH21 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' cell
  214. 'NATU' 'DISCRET' ;
  215. XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' zero
  216. 'NATURE' 'DISCRET' ;
  217. XO21 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1
  218. 'O2' cell1 'NATU' 'DISCRET' ;
  219. XN = XH21 'ET' XH2O1 'ET' XO21 ;
  220. CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0D0) ;
  221. XN21 = CHUN '-' ('PSCAL' XN CHUN ('MOTS' 'H2' 'O2' 'H2O')
  222. ('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
  223. XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
  224. XN = XN 'ET' XN21 ;
  225.  
  226. *
  227. **** We compute the mass fractions
  228. *
  229.  
  230. MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' Mh2
  231. 'NATURE' 'DISCRET') 'ET'
  232. ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' Mo2
  233. 'NATURE' 'DISCRET') 'ET'
  234. ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' Mh2o
  235. 'NATURE' 'DISCRET') 'ET'
  236. ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'N2' Mn2
  237. 'NATURE' 'DISCRET') ;
  238.  
  239. MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  240. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  241.  
  242. YH21 = (XH21 '*' MH2) '/' MTOT ;
  243. YO21 = (XO21 '*' MO2) '/' MTOT ;
  244. YN21 = (XN21 '*' MN2) '/' MTOT ;
  245. YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;
  246.  
  247. *
  248. **** Burnt gas
  249. * Here X refers to the mole producted per 1 mole of unburned gas
  250. * N.B. The number of moles changes after the combustion!
  251. *
  252.  
  253. * H2MAG = 1 in the region where H2 is the majority species.
  254. * Conversely H2MAG = 0
  255.  
  256. H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
  257. ('EXCO' 'O2' XN) ;
  258. O2MAG = CHUN '-' H2MAG ;
  259.  
  260. XO22a = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' zero ;
  261. XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
  262.  
  263. XO22a = XO22a ;
  264. XH22a = XH22a ;
  265.  
  266. XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' zero ;
  267. XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (xh21 '-' xh22b))) ;
  268.  
  269. XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
  270. XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;
  271.  
  272. XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
  273. XN22 = XN21 ;
  274. *
  275. **** mtot = weight of a mole of unburned gas
  276. * mtot2 = weight of a mole of unburned gas after the combustion
  277. * Of course mtot2 = mtot
  278. *
  279.  
  280. MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
  281. ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  282. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  283.  
  284. mtotcel = 'MAXIMUM' MTOT ;
  285.  
  286.  
  287. 'SI' ( ('MAXIMUM' (mtot2 '-' mtot) 'ABS') '>' (zero '*' mtotcel)) ;
  288. 'MESSAGE' 'Problem in the computation of the burnt gas (1).' ;
  289. 'ERREUR' 21 ;
  290. 'FINSI' ;
  291.  
  292. *
  293. **** xtot2 = mole of burnt gas per mole of unburned gas
  294. *
  295.  
  296. xtot2 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
  297. CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  298. ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
  299.  
  300. *
  301. **** X = real molar fractions
  302. *
  303.  
  304. xh22 = xh22 '/' xtot2 ;
  305. xo22 = xo22 '/' xtot2 ;
  306. xh2o2 = xh2o2 '/' xtot2 ;
  307. xn22 = xn22 '/' xtot2 ;
  308.  
  309. LISXH2F = LISXH2F 'ET' (PROG (MAXI XH22)) ;
  310. LISXO2F = LISXO2F 'ET' (PROG (MAXI XO22)) ;
  311. LISXWAF = LISXWAF 'ET' (PROG (MAXI XH2O2)) ;
  312. LISXN2F = LISXN2F 'ET' (PROG (MAXI XN22)) ;
  313.  
  314.  
  315. xtot3 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
  316. CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  317. ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
  318.  
  319. 'SI' ( ('MAXIMUM' (xtot3 '-' 1.0D0) 'ABS' ) '>'
  320. zero) ;
  321. 'MESSAGE' 'Problem in the computation of the burnt gas (2).' ;
  322. 'ERREUR' 21 ;
  323. 'FINSI' ;
  324.  
  325. *
  326. **** mtot2 = weight of 1 mole of burned gas
  327. *
  328.  
  329. MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
  330. ('MOTS' 'H2' 'O2' 'H2O' 'N2')
  331. ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
  332.  
  333. *
  334. **** Mass fractions of burned gas after the complete combustion
  335. *
  336.  
  337. YH22 = (XH22 '*' MH2) '/' MTOT2 ;
  338. YO22 = (XO22 '*' MO2) '/' MTOT2 ;
  339. YN22 = (XN22 '*' MN2) '/' MTOT2 ;
  340. YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
  341.  
  342. LISYH2F = LISYH2F 'ET' (PROG (MAXI YH22)) ;
  343. LISYO2F = LISYO2F 'ET' (PROG (MAXI YO22)) ;
  344. LISYWAF = LISYWAF 'ET' (PROG (MAXI YH2O2)) ;
  345. LISYN2F = LISYN2F 'ET' (PROG (MAXI YN22)) ;
  346.  
  347.  
  348. *
  349. **** The conservative variables
  350. *
  351.  
  352. * Gas constant
  353.  
  354. RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
  355. (('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
  356. (('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
  357. (('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
  358. ;
  359.  
  360. * Density
  361.  
  362. RNINI = PINI '/' (RINI '*' TINI) ;
  363. RN = 'COPIER' RNINI ;
  364.  
  365. * Internal energy
  366.  
  367. eini = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
  368. Tfois = Tini ;
  369.  
  370. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  371.  
  372. acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  373. &BL1)) '+'
  374. (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  375. &BL1)) '+'
  376. (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  377. &BL1)) '+'
  378. (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  379. &BL1)) ;
  380.  
  381. eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
  382.  
  383. Tfois = Tfois * tini ;
  384.  
  385. 'FIN' BL1 ;
  386.  
  387. RETNINI = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
  388. ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
  389. RETN = RETNINI ;
  390. GN = RN '*' WINI ;
  391. YN = YH21 '+' YO21 '+' YH2O1 ;
  392. RYN = RN '*' YN ;
  393. *
  394. **** K0/DX = chemical source term
  395. *
  396. * i.e. dcsi/dt = K0 '/' DX . {criterium function}
  397. * where
  398. * {criterium function} = 1.0 or 0.0
  399. *
  400. * d(rho etot)/dt = -rho (dcsi/dt) (h_f '-' h_i)
  401. *
  402. * h_f reference enthalpy (0K) of the burned gas
  403. * h_i reference enthalpy (0K) of the unburned gas
  404. *
  405. **** Mass fraction of H2 before and after the combustion
  406. *
  407. * CSIMAX = 1.0 -> complete combustion
  408. *
  409.  
  410. YININ = yh21 ;
  411. YFINN = YININ '+' ((yh22 '-' yh21) '*' csimax) ;
  412. K0N = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 500 ;
  413.  
  414. *
  415. **** Passive scalars
  416. *
  417.  
  418. RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
  419. ('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
  420. ('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
  421.  
  422. * Test of prim
  423. VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
  424. PCEL = 'MAXIMUM' PINI ;
  425. TCEL = 'MAXIMUM' TINI ;
  426. 'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
  427. 'MESSAGE' 'Problem in the computation of the energy (1).'
  428. 'ERREUR' 21 ;
  429. 'FINSI' ;
  430. 'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
  431. 'MESSAGE' 'Problem in the computation of the energy (2).'
  432. 'ERREUR' 21 ;
  433. 'FINSI' ;
  434.  
  435. PINSI1 = MAXI P ;
  436. TINSI1 = MAXI T ;
  437. RINSI1 = MAXI RN ;
  438. GAMINSI = MAXI GAMN ;
  439. GAMAINSI = 1 + (MAXI (P / RETN)) ;
  440. RGAS = MAXI (P / (RN * T)) ;
  441.  
  442. LRINSI = LRINSI 'ET' (PROG RINSI1) ;
  443. LTINSI = LTINSI 'ET' (PROG TINSI1) ;
  444. LPINSI = LPINSI 'ET' (PROG PINSI1) ;
  445. LRGINSI = LRGINSI 'ET' (PROG RGAS) ;
  446. LGINSI = LGINSI 'ET' (PROG GAMINSI) ;
  447. LGAINSI = LGAINSI 'ET' (PROG GAMAINSI) ;
  448. LCINSI = LCINSI 'ET' (PROG ((GAMINSI * PINSI1 / RINSI1) ** 0.5)) ;
  449.  
  450. *
  451. * VN state
  452. *
  453. * VN state via DETO operator (not really complete combustion...)
  454. *
  455. CHPO1 = 'MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 6
  456. 'H2' (MAXI XH21)
  457. 'O2' (MAXI XO21)
  458. 'N2' (1. - ((MAXI XH21) + (MAXI XO21)))
  459. 'H2O' 0.0
  460. 'P' (MAXI PINI) 'T' (MAXI TINI) ;
  461.  
  462. CHPCJ CHPVN CHPAICC = 'DETO' CHPO1 ;
  463.  
  464. PVN = (EXCO CHPVN 'PZND') ;
  465. RVN = (EXCO CHPVN 'RZND') ;
  466. TVN = (EXCO CHPVN 'TZND') ;
  467. DVN = (EXCO CHPCJ 'VCJ') ;
  468. * VVN via RH conditions
  469. VVN = DVN '*' (RNINI / RVN) ;
  470. VVN = DVN - VVN ;
  471. VVN = (NOMC VVN 'UX') +
  472. ('MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 1 'UY' 0.0) ;
  473.  
  474. * Internal energy
  475.  
  476. EVN = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
  477. Tfois = TVN ;
  478.  
  479. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  480.  
  481. acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  482. &BL1)) '+'
  483. (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  484. &BL1)) '+'
  485. (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  486. &BL1)) '+'
  487. (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  488. &BL1)) ;
  489.  
  490. eVN = eVN '+' (acel '*' (Tfois '/' &BL1)) ;
  491.  
  492. Tfois = Tfois * tVN ;
  493.  
  494. 'FIN' BL1 ;
  495.  
  496. RETNVN = RVN '*' (eVN '+' (0.5 '*' ('PSCAL' VVN VVN
  497. ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
  498.  
  499. V1VN P1VN T1VN YNVN SNVN GAMVN = 'PRIM' 'PERFTEMP' PGAZ RVN (RVN * VVN)
  500. RETNVN (RVN * YN) RSN ;
  501.  
  502. ERRO = MAXI ((TVN - T1VN) / T1VN) ABS ;
  503.  
  504. 'SI' (ERRO > 1.0d-3) ;
  505. MESS 'Problem in VN' ;
  506. ERREUR 21 ;
  507. FINS ;
  508. *
  509. * We check that the RH on the energy is satisfied
  510. *
  511. VVN = DVN '*' (RNINI / RVN) ;
  512. HVN = (EVN + (PVN / RVN)) + (0.5 * VVN * VVN) ;
  513. HINI = (EINI + (PINI / RNINI)) + (0.5 * DVN * DVN) ;
  514. ERRO = MAXI ((HVN - HINI) / HVN) ABS ;
  515. 'SI' (ERRO > 1.0d-2) ;
  516. MESS 'Problem in VN (2)' ;
  517. ERREUR 21 ;
  518. FINS ;
  519.  
  520. PVN1 = MAXI PVN ;
  521. TVN1 = MAXI TVN ;
  522. RVN1 = MAXI RVN ;
  523. GAMVN = MAXI GAMVN ;
  524. GAMAVN = 1 + (MAXI (PVN / (RVN * EVN))) ;
  525. RGAS = MAXI (PVN / (RVN * TVN)) ;
  526.  
  527. LRVN = LRVN 'ET' (PROG RVN1) ;
  528. LTVN = LTVN 'ET' (PROG TVN1) ;
  529. LPVN = LPVN 'ET' (PROG PVN1) ;
  530. LRGVN = LRGVN 'ET' (PROG RGAS) ;
  531. LGVN = LGVN 'ET' (PROG GAMVN) ;
  532. LGAVN = LGAVN 'ET' (PROG GAMAVN) ;
  533. LCVN = LCVN 'ET' (PROG ((GAMVN * PVN1 / RVN1) ** 0.5)) ;
  534.  
  535. *
  536. ***** Combustion via the operator FLAM
  537. *
  538.  
  539. DELTAXN = ('DOMA' $DOM1 'VOLUME') '**' (1. '/' 2.) ;
  540. YN = RYN '/' RN ;
  541.  
  542. *
  543. **** We want the complete combustion in DOM1.
  544. * We can use the operator 'FLAM', with epsilon=0 and
  545. * DELTAT >> DELTATC
  546. * where DELTATC is the characteristic time step linked to
  547. * the source term.
  548. *
  549.  
  550. SN = RSN '/' RN ;
  551. K0N = 'EXCO' 'K0' SN ;
  552. YININ = 'EXCO' 'H2IN' SN 'H2' ;
  553. YFINN = 'EXCO' 'H2FI' SN 'H2' ;
  554.  
  555. DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
  556. LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
  557. LCOEF = 'PROG' 1.0 0.5 -1.0 ;
  558.  
  559. DELTARE DELTARY = 'FLAM' 'CREBCOM2' $DOM1 PGAZ LMOT1 LCOEF
  560. RN YN YININ YFINN K0N DELTAXN
  561. 0.0 (1D4 '*' DELTATC) 0.3 ;
  562.  
  563. Q = ((MAXI (DELTARE / RN))) ;
  564. Q1 = (('NOMC' 'SCAL' (yh21 '-' yh22)) * (PGAZ .'H2 ' . 'H0K')) +
  565. (('NOMC' 'SCAL' (yo21 '-' yo22)) * (PGAZ .'O2 ' . 'H0K')) +
  566. (('NOMC' 'SCAL' (yn21 '-' yn22)) * (PGAZ .'N2 ' . 'H0K')) +
  567. (('NOMC' 'SCAL' (yh2o1 '-' yh2o2)) * (PGAZ .'H2O ' . 'H0K')) ;
  568.  
  569. ERRO = MAXI ((Q - Q1) / Q1) ABS ;
  570. 'SI' (ERRO > 1.0d-2) ;
  571. MESS 'Problem in Q' ;
  572. ERREUR 21 ;
  573. FINS ;
  574.  
  575. LQ = LQ 'ET' (PROG Q) ;
  576.  
  577. RYN = RYN '+' DELTARY ;
  578. YN = RYN '/' RN ;
  579. RETN = RETN '+' DELTARE ;
  580. *
  581. *
  582.  
  583. *
  584. VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
  585. PAICC1 = MAXI PN ;
  586. TAICC1 = MAXI TN ;
  587. RAICC1 = MAXI RN ;
  588. GAMAICC = MAXI GAMMAN ;
  589. GAMAAICC = 1 + (MAXI (PN / RETN)) ;
  590. RGAS = MAXI (PN / (RN * TN)) ;
  591. *
  592. * AICC, CJ via DETO operator (not really complete combustion...)
  593. *
  594. CHPO1 = 'MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 6
  595. 'H2' (MAXI XH21)
  596. 'O2' (MAXI XO21)
  597. 'N2' (1. - ((MAXI XH21) + (MAXI XO21)))
  598. 'H2O' 0.0
  599. 'P' (MAXI PINI) 'T' (MAXI TINI) ;
  600.  
  601. CHPCJ CHPZND CHPAICC = 'DETO' CHPO1 ;
  602. PAICC2 = MAXI (EXCO CHPAICC 'PAIC') ;
  603. RAICC2 = MAXI (EXCO CHPAICC 'RAIC') ;
  604. TAICC2 = MAXI (EXCO CHPAICC 'TAIC') ;
  605. *
  606. LERRO = ABS (PROG ((PAICC1 - PAICC2) / PAICC1)
  607. ((RAICC1 - RAICC2) / RAICC1)
  608. ((TAICC1 - TAICC2) / TAICC1)) ;
  609. ERRO = MAXI LERRO ;
  610. 'SI' (ERRO > 1.0d-2) ;
  611. MESS 'Problem in AICC' ;
  612. ERREUR 21 ;
  613. FINS ;
  614. LRAICC = LRAICC 'ET' (PROG RAICC1) ;
  615. LTAICC = LTAICC 'ET' (PROG TAICC1) ;
  616. LPAICC = LPAICC 'ET' (PROG PAICC1) ;
  617. LRGAICC = LRGAICC 'ET' (PROG RGAS) ;
  618. LGAICC = LGAICC 'ET' (PROG GAMAICC) ;
  619. LGAAICC = LGAAICC 'ET' (PROG GAMAAICC) ;
  620. LCAICC = LCAICC 'ET' (PROG ((GAMAICC * PAICC1 / RAICC1) ** 0.5)) ;
  621. *
  622. * CJ
  623. *
  624. PCJ = (EXCO CHPCJ 'PCJ') ;
  625. RCJ = (EXCO CHPCJ 'RCJ') ;
  626. TCJ = (EXCO CHPCJ 'TCJ') ;
  627. DCJ = (EXCO CHPCJ 'VCJ') ;
  628.  
  629. Tfois = TCJ ;
  630. eCJ = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
  631.  
  632.  
  633. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  634.  
  635. acel = (('NOMC' 'SCAL' yh22) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  636. &BL1)) '+'
  637. (('NOMC' 'SCAL' yo22) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  638. &BL1)) '+'
  639. (('NOMC' 'SCAL' yh2o2) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  640. &BL1)) '+'
  641. (('NOMC' 'SCAL' yn22) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  642. &BL1)) ;
  643.  
  644. eCJ = eCJ '+' (acel '*' (Tfois '/' &BL1)) ;
  645.  
  646. Tfois = Tfois * TCJ ;
  647.  
  648. 'FIN' BL1 ;
  649.  
  650. RETNCJ = RCJ '*' eCJ ;
  651.  
  652. RWCJ = RCJ * VN * 0 ;
  653.  
  654. V1CJ P1CJ T1CJ YNCJ SNCJ GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RCJ RWCJ
  655. RETNCJ (RCJ * YN) (RCJ * SN) ;
  656. CCJ = (GAMMAN * (P1CJ / RCJ)) '**' 0.5 ;
  657. VCJ = DCJ * (RNINI / RCJ) ;
  658. erro = MAXI ((VCJ / CCJ) - 1.0) ABS ;
  659.  
  660. * We check weather the speed in the shock reference is equal to the speed of
  661. * sound.
  662. 'SI' (erro > 2.0d-2) ;
  663. MESS 'Problem in CJ' ;
  664. ERREUR 21 ;
  665. FINS ;
  666. GAMACJ = 1 + (MAXI (PCJ / RETNCJ)) ;
  667.  
  668. LRCJ = LRCJ 'ET' (PROG (maxi RCJ)) ;
  669. LTCJ = LTCJ 'ET' (PROG (maxi TCJ)) ;
  670. LPCJ = LPCJ 'ET' (PROG (maxi PCJ)) ;
  671. LRGCJ = LRGCJ 'ET' (PROG RGAS) ;
  672. LGCJ = LGCJ 'ET' (PROG (MAXI GAMMAN)) ;
  673. LGACJ = LGACJ 'ET' (PROG GAMACJ) ;
  674. LCCJ = LCCJ 'ET' (PROG (maxi CCJ)) ;
  675.  
  676. *
  677. * AIBCC
  678. *
  679. * Gas constant
  680. *
  681. YO2 = 'MAXIMUM' ('EXCO' 'O2' YN) ;
  682. YH2 = 'MAXIMUM' ('EXCO' 'H2' YN) ;
  683. YH2O = 'MAXIMUM' ('EXCO' 'H2O' YN) ;
  684. YN2 = (-1 '*' ((YO2 + YH2 + YH2O) '-' 1)) ;
  685. R = ((PGAZ . 'H2O ' . 'R') '*' YH2O) '+'
  686. ((PGAZ . 'O2 ' . 'R') '*' YO2) '+'
  687. ((PGAZ . 'H2 ' . 'R') '*' YH2) '+'
  688. ((PGAZ . 'N2 ' . 'R') '*' YN2) ;
  689. T0 = 500 ;
  690. T1 = 3500 ;
  691. hfin = 'MAXIMUM' ((RETN '/' RN) '+' (PINI '/' RN)) ;
  692. * hfin = EINI '+' q '+' (PINI '/' RN) ;
  693. *
  694. * Regula falsi
  695. *
  696. 'REPETER' BLITER 100 ;
  697. T = (T0 '+' T1) '/' 2. ;
  698. eini = 0.0 ;
  699. Tfois = T ;
  700. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  701. acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  702. &BL1)) '+'
  703. (YN2 '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  704. &BL1)) '+'
  705. (YO2 '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  706. &BL1)) '+'
  707. (YH2 '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  708. &BL1)) ;
  709. eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
  710. Tfois = Tfois * T ;
  711. 'FIN' BL1 ;
  712. hini = eini '+' (R * T) ;
  713. 'MESS' ('CHAINE' 'h, T ' hini ' ' T) ;
  714. 'SI' (< ('ABS' (hini '-' hfin)) 1.0D-4) ;
  715. 'MESSAGE' OK ;
  716. 'QUITTER' BLITER ;
  717. 'SINON' ;
  718. 'SI' (hini > hfin) ;
  719. T1 = T ;
  720. 'SINON' ;
  721. T0 = T ;
  722. 'FINSI' ;
  723. 'FINSI' ;
  724. 'FIN' BLITER ;
  725. *
  726. TAIBC1=T;
  727. PAIBC1=(MAXI PINI) ;
  728. RNAIBCC = PINI '/' (R '*' T) ;
  729. RAIBC1=(MAXI RNAIBCC) ;
  730. eini = 0.0 ;
  731. Tfois = T ;
  732. 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
  733. acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
  734. &BL1)) '+'
  735. (YN2 '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
  736. &BL1)) '+'
  737. (YO2 '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
  738. &BL1)) '+'
  739. (YH2 '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
  740. &BL1)) ;
  741. eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
  742. Tfois = Tfois * T ;
  743. 'FIN' BL1 ;
  744. RETAIBCC = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' eini) '*'
  745. RNAIBCC ;
  746. RYNAIBCC = RNAIBCC '*' YN ;
  747. VN2 P2 T2 Y2 S2 GAMN = 'PRIM' 'PERFTEMP' PGAZ RNAIBCC GN RETAIBCC
  748. RYNAIBCC RSN ;
  749. TAIBC2= MAXI T2 ;
  750. PAIBC2= MAXI P2 ;
  751. LERRO = ABS (PROG ((PAIBC1 - PAIBC2) / PAIBC1)
  752. ((TAIBC1 - TAIBC2) / TAIBC1)) ;
  753. ERRO = MAXI LERRO ;
  754. 'SI' (ERRO > 1.0d-2) ;
  755. MESS 'Problem in AIBCC' ;
  756. ERREUR 21 ;
  757. FINS ;
  758. LRAIBC = LRAIBC 'ET' (PROG RAIBC1) ;
  759. LTAIBC = LTAIBC 'ET' (PROG TAIBC1) ;
  760. LPAIBC = LPAIBC 'ET' (PROG PAIBC1) ;
  761. GAMAIBC = MAXI GAMN ;
  762. GAMAAIBC = 1 + (MAXI (P2 / RETAIBCC)) ;
  763. RGAS = MAXI (P2 / (RNAIBCC * T2)) ;
  764. LRGAIBC = LRGAIBC 'ET' (PROG RGAS) ;
  765. LGAIBC = LGAIBC 'ET' (PROG GAMAIBC) ;
  766. LGAAIBC = LGAAIBC 'ET' (PROG GAMAAIBC) ;
  767. LCAIBC = LCAIBC 'ET' (PROG ((GAMAIBC * PAIBC1 / RAIBC1) ** 0.5)) ;
  768. *
  769. *************************************************************************
  770. *************************************************************************
  771. *********************** END OF THE LOOP *********************************
  772. *************************************************************************
  773. *************************************************************************
  774. *
  775. FIN BLCALC ;
  776.  
  777. *
  778. **** Some LaTex tables.
  779. *
  780. SI FAUX ;
  781. OPTI ECHO 0 ;
  782. REPE BLCALC (DIME LISTXH2) ;
  783. *
  784. R = EXTR LRINSI &BLCALC ;
  785. P = EXTR LPINSI &BLCALC ;
  786. T = EXTR LTINSI &BLCALC ;
  787. C = EXTR LCINSI &BLCALC ;
  788. GAM = EXTR LGINSI &BLCALC ;
  789. GAMA = EXTR LGAINSI &BLCALC ;
  790. RGAS = EXTR LRGINSI &BLCALC ;
  791. AA = CHAI &blcalc ' & Inside & ' 'FORMAT' '(E9.3)' R ' & '
  792. 'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
  793. 'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
  794. 'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
  795. MESS AA ;
  796. *
  797. R = EXTR LRAICC &BLCALC ;
  798. P = EXTR LPAICC &BLCALC ;
  799. T = EXTR LTAICC &BLCALC ;
  800. C = EXTR LCAICC &BLCALC ;
  801. GAM = EXTR LGAICC &BLCALC ;
  802. GAMA = EXTR LGAAICC &BLCALC ;
  803. RGAS = EXTR LRGAICC &BLCALC ;
  804. AA = CHAI &blcalc ' & AICC & ' 'FORMAT' '(E9.3)' R ' & '
  805. 'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
  806. 'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
  807. 'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
  808. MESS AA ;
  809. *
  810. R = EXTR LRAIBC &BLCALC ;
  811. P = EXTR LPAIBC &BLCALC ;
  812. T = EXTR LTAIBC &BLCALC ;
  813. C = EXTR LCAIBC &BLCALC ;
  814. GAM = EXTR LGAIBC &BLCALC ;
  815. GAMA = EXTR LGAAIBC &BLCALC ;
  816. RGAS = EXTR LRGAIBC &BLCALC ;
  817. AA = CHAI &blcalc ' & AIBCC & ' 'FORMAT' '(E9.3)' R ' & '
  818. 'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
  819. 'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
  820. 'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
  821. MESS AA ;
  822. *
  823. R = EXTR LRCJ &BLCALC ;
  824. P = EXTR LPCJ &BLCALC ;
  825. T = EXTR LTCJ &BLCALC ;
  826. C = EXTR LCCJ &BLCALC ;
  827. GAM = EXTR LGCJ &BLCALC ;
  828. GAMA = EXTR LGACJ &BLCALC ;
  829. RGAS = EXTR LRGCJ &BLCALC ;
  830. AA = CHAI &blcalc ' & CJ & ' 'FORMAT' '(E9.3)' R ' & '
  831. 'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
  832. 'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
  833. 'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
  834. MESS AA ;
  835. *
  836. R = EXTR LRVN &BLCALC ;
  837. P = EXTR LPVN &BLCALC ;
  838. T = EXTR LTVN &BLCALC ;
  839. C = EXTR LCVN &BLCALC ;
  840. GAM = EXTR LGVN &BLCALC ;
  841. GAMA = EXTR LGAVN &BLCALC ;
  842. RGAS = EXTR LRGVN &BLCALC ;
  843. AA = CHAI &blcalc ' & VN & ' 'FORMAT' '(E9.3)' R ' & '
  844. 'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
  845. 'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
  846. 'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
  847. MESS AA ;
  848. FIN BLCALC ;
  849. OPTI ECHO 1 ;
  850. OPTI ECHO 0 ;
  851. REPE BLCALC (DIME LISTXH2) ;
  852. *
  853. XH2 = EXTR LISXH2F &BLCALC ;
  854. XO2 = EXTR LISXO2F &BLCALC ;
  855. XWA = EXTR LISXWAF &BLCALC ;
  856. XN2 = EXTR LISXN2F &BLCALC ;
  857. YH2 = EXTR LISYH2F &BLCALC ;
  858. YO2 = EXTR LISYO2F &BLCALC ;
  859. YWA = EXTR LISYWAF &BLCALC ;
  860. YN2 = EXTR LISYN2F &BLCALC ;
  861. Q = EXTR LQ &BLCALC ;
  862. * AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)' XH2 ' & '
  863. * 'FORMAT' '(E9.3)' XO2 ' & ' 'FORMAT' '(E9.3)' XWA ' & '
  864. * 'FORMAT' '(E9.3)' XN2 ' & ' 'FORMAT' '(E9.3)' YH2 ' & '
  865. * 'FORMAT' '(E9.3)' YO2 ' & ' 'FORMAT' '(E9.3)' YWA ' & '
  866. * 'FORMAT' '(E9.3)' YN2 ' \\';
  867. AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)' XH2 ' & '
  868. 'FORMAT' '(E9.3)' XO2 ' & ' 'FORMAT' '(E9.3)' XWA ' & '
  869. 'FORMAT' '(E9.3)' YH2 ' & '
  870. 'FORMAT' '(E9.3)' YO2 ' & ' 'FORMAT' '(E9.3)' YWA ' & '
  871. 'FORMAT' '(E9.3)' Q ' \\';
  872. MESS AA ;
  873. MESS '\hline' ;
  874. FIN BLCALC ;
  875. OPTI ECHO 1 ;
  876. FINSI ;
  877. FIN ;
  878.  
  879.  
  880.  

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