Télécharger la plaquette Cast3M


fr en

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.  
  881.  
  882.  

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