Télécharger la plaquette Cast3M


fr en

Télécharger tubedeto2d1.dgibi

Retour à la liste

Numérotation des lignes : 

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

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