Télécharger tubedeto2d1.dgibi

Retour à la liste

* fichier :  tubedeto2d1.dgibi
************************************************************************
************************************************************************
****************************************************************
*****      PROPAGATION D UNE DETONATION DANS UN TUBE       *****
*****      MODELE CREBCOM                                  *****
*****      A. BECCANTINI, SFME/LTMF, 15.01.02              *****
****************************************************************
*
* The file is divided into 5 parts
*
* 1) mesh
* 2) initial conditions and gas properties
* 3) parameters used in the computation
* 4) the main part (where the computation is realised)
* 5) the post-treatment
*
 
   'OPTION' 'ECHO' 0 'DIME' 2
         'TRAC' 'X' ;
 
****************************************************************
*****      1D mesh                                         *****
****************************************************************
*
*
*  A6 ---- A5 ------------------- A4 ----------------------- A8
*   |      |                       |                         |
*   |  AR  |        ZONE1          |         ZONE2           |
*   |      |                       |                         |
*  A1 ---- A2 ------------------- A3 ----------------------- A7
*   |      |                       |                         |
*   |  L1  |          L2           |           L3            |
*   |<---->|<--------------------->|<----------------------->|
*
*  AR = activation region
*  ZONE1
*  ZONE2
*
 
 RAF = 1  ;
 
 NL2SL1 = 49 ;
 NL3SL1 = 100 ;
 L1 = 0.4 ;
 
 L2 = NL2SL1 '*' L1 ;
 L3 = NL3SL1 '*' L1 ;
 
* If  NL2SL1 = 49 and NL3SL1 = 50, then L2 = 19.6, L3 = 20.0 ;
 
 N1 = RAF ;
 N2 = NL2SL1 '*' N1 ;
 N3 = NL3SL1 '*' N1 ;
 
* For sake of simplicity, we will only consider one mesh in y-direction
 
 DX = L1 '/' N1 ;
 NY = 4 ;
 NZ = NY ;
 H1 = NY '*' DX ;
 
 A1 = 0.0 0.0 ;
 A2 = L1  0.0 ;
 A3 = (L1 '+' L2) 0.0 ;
 A4 = (L1 '+' L2) H1 ;
 A5 = L1  H1 ;
 A6 = 0.0 H1 ;
 A7 = (L1 '+' L2 '+' L3) 0.0 ;
 A8 = (L1 '+' L2 '+' L3) H1 ;
 
 'OPTION' 'ELEM' 'SEG2' ;
 A1A2 = A1 'DROIT' N1 A2 ;
 A2A3 = A2 'DROIT' N2 A3 ;
 A3A4 = A3 'DROIT' NY A4 ;
 A4A5 = A4 'DROIT' N2 A5 ;
 A5A2 = A5 'DROIT' NY A2 ;
 A2A5 = 'INVERSE' A5A2 ;
 A5A6 = A5 'DROIT' N1 A6 ;
 A6A1 = A6 'DROIT' NY A1 ;
 A3A7 = A3 'DROIT' N3 A7 ;
 A7A8 = A7 'DROIT' NY A8 ;
 A8A4 = A8 'DROIT' N3 A4 ;
 A4A3 = 'INVERSE' A3A4 ;
 
*
**** 2D mesh
*
 
* DOM1 is the activation region
* DOM2 is zone 1
* DOM3 is zone 2
 
 'OPTION' 'ELEM' 'QUA4' ;
 DOM1 = 'TRANSLATION'  A2A5 N1 ((-1.0 '*' L1) 0.0) ;
 DOM2 = 'TRANSLATION'  A2A5 N2 (L2 0.0) ;
 DOM3 = 'TRANSLATION'  A3A4 N3 (L3 0.0) ;
 
* DOMTOT = total region
 
 DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;
 
*
**** The tables 'DOMAINE'
*
 
 $DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
 $DOM1   = 'MODELISER' DOM1 'EULER' ;
 $DOM2   = 'MODELISER' DOM2 'EULER' ;
 $DOM3   = 'MODELISER' DOM3 'EULER' ;
 
 TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
 TDOM1   = 'DOMA' $DOM1   'VF' ;
 TDOM2   = 'DOMA' $DOM2   'VF' ;
 TDOM3   = 'DOMA' $DOM3   'VF' ;
 
 QDOMTOT = TDOMTOT  . 'QUAF' ;
 QDOM1   = TDOM1    . 'QUAF' ;
 QDOM2   = TDOM2    . 'QUAF' ;
 QDOM3   = TDOM3    . 'QUAF' ;
 
 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM1 ;
 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM2 ;
 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM3 ;
 
*
**** MOD1 = model (created just to display the CHAMELEMs)
*
 
 MOD1 =  'MODELISER'  ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
*
**** Line for graphics
*
 
 'OPTION' 'ELEM' 'SEG2' ;
 P1 = (DX '/' 2) (DX '/' 2);
 P2 = ((L1 '+' L2 '+' L3) '-' (DX '/' 2)) (DX '/' 2);
 LIGEVO = (P1 'DROIT' (N1 '+' N2 '+' N3 '-' 1) P2) 'COULEUR' 'ROUG'  ;
 'ELIMINATION' LIGEVO ('DOMA' $DOMTOT 'CENTRE') (0.001 '/' RAF) ;
 
* 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;
 
****************************************************************
****************************************************************
*****      Initial conditions and gas properties.          *****
****************************************************************
****************************************************************
 
*
**** 0 from a numerical point of view
*
 
 zero = 1.0d-9 ;
 
*
**** (Non-Uniform) initial conditions.(DOM1 et DOM2)
*    We compute the conservative variables.
*
*
* Contact discontinuity between DOM2 and DOM3 (same pressure,
* different temperature)
*
 
 PINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET' 
        ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 99700 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 99700
        'NATU' 'DISCRET') ;
 TINI =  ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
         ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 285. 'NATU' 'DISCRET') 'ET'
         ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 400
        'NATU' 'DISCRET') ; 
 WINI = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 2 'UX' 0.0 'UY' 0.0;
 
* We suppose that the combustion is complete everywhere, i.e. the
* progress variable is 1 after the detonation transition.
 
 CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0 ;
 
*
**** From initial conditions on molar fractions, we duduce
*    the data for the virtual mixture of burned/unburned
*    gases
 
*
**** Molar fractions (before combustion)
*
 
 XH21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'H2' 0.3 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'H2' 0.2
         'NATU' 'DISCRET') ;
 XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' zero
         'NATURE' 'DISCRET' ;
 XO21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'O2' 0.147 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'O2' 0.2
         'NATU' 'DISCRET') ;
 XN = XH21 'ET' XH2O1 'ET' XO21 ;
 CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0D0) ;
 XN21 =  CHUN '-' ('PSCAL' XN CHUN  ('MOTS' 'H2' 'O2' 'H2O')
         ('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
 XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
 XN = XN 'ET' XN21 ;
 
*
***** Gas properties
*
 
 PGAZ = 'TABLE' ;
 
* Polynomial degree of specific heats
 
 PGAZ . 'NORD' = 4 ;
 
* Species explicitly treated in the Euler Equations
 
 PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;
 
* Species non explicitly treated
 
 PGAZ . 'ESPNEULE' = 'N2  ';
 
* Single gases properties
 
 PGAZ .  'H2  ' = 'TABLE'  ;
 PGAZ .  'H2O ' = 'TABLE'  ;
 PGAZ .  'N2  ' = 'TABLE'  ;
 PGAZ .  'O2  ' = 'TABLE'  ;
 
* R (J/Kg/K)
 
 mH2 = 2. '*' 1.00797E-3 ;
 mo2 = 2. '*' 15.9994E-3 ;
 mH2O = mh2 '+' (0.5 '*' mo2) ;
 mN2 = 2 '*' 14.0067E-3 ;
 RGAS = 8.31441 ;
 
 PGAZ .  'H2  ' . 'R' = RGAS '/' mh2 ;
 PGAZ .  'H2O ' . 'R' = RGAS '/' mh2o ;
 PGAZ .  'N2  ' . 'R' = RGAS '/' mn2 ;
 PGAZ .  'O2  ' . 'R' = RGAS '/' mo2 ;
 
* Polynomials regressions coefficients
 
 PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
 PGAZ .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
                              -1.82753232E-08 2.44485692E-12 ;
 PGAZ .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
                              8.78233606E-09 -3.05514485E-13 ;
 PGAZ .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
                              2.33636971E-08 -1.53304905E-12;
 
* Formation enthalpies at 0K (J/Kg))
*
*     h_i(0K)  = h_i(T0) '-' \int_0^{T0} cp_i(x) dx 
*              = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
*     
 
 PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
 PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
 PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;
 
* Passive scalars names
 
 PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0  ';
 
*  
 
* 'OPTION' 'SAUVER' 'gas_properties.sauv' ;
* 'SAUVER' PGAZ ;
 
* 
**** We compute the mass fractions
*
 
 MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' Mh2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' Mo2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' Mh2o
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'N2' Mn2
         'NATURE' 'DISCRET') ;
 
 MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
 
 YH21 = (XH21 '*' MH2) '/' MTOT ;
 YO21 = (XO21 '*' MO2) '/' MTOT ;
 YN21 = (XN21 '*' MN2) '/' MTOT ;
 YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;
 
*
**** Burned gas
*    Here x refers to the mole producted per 1 mole of unburned gas
*    N.B. Number of moles changes after the combustion!
*
 
*
* H2MAG = 1 in the region where H2 is the majority species.
* Conversely H2MAG = 0
*
 
 H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
    ('EXCO' 'O2'  XN) ;
 O2MAG = CHUN '-' H2MAG ;
 
 XO22a = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' zero ;
 XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
 
 XO22a = XO22a ;
 XH22a = XH22a ;
 
 XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' zero ;
 XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (xh21 '-' xh22b))) ;
 
 XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
 XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;
 
 XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
 XN22 = XN21 ;
 
*
**** mtot = weight of a mole  of unburned gas
*    mtot2 = weight of a mole of unburned gas after the combustion
*    Of course mtot2 = mtot
*
 
 MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
   ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
 
 mtotcel = 'MAXIMUM' MTOT ;
 
 'SI' ( ('MAXIMUM'  (mtot2 '-' mtot) 'ABS') '>' (zero '*' mtotcel)) ;
    'MESSAGE' 'Probleme dan le calcul du gaz brule' ;
    'ERREUR' 21 ;
 'FINSI' ;
 
*
**** xtot2 = mole of burnet gas per mole of unburned gas
*
 
 xtot2 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
   CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
 
*
**** x = real molar fractions
*
 
 xh22 = xh22 '/' xtot2 ;
 xo22 = xo22 '/' xtot2 ;
 xh2o2 = xh2o2 '/' xtot2 ;
 xn22 = xn22 '/' xtot2 ;
 
 
 xtot3 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
   CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
 
 'SI' ( ('MAXIMUM' (xtot3 '-' 1.0D0) 'ABS' ) '>'
    zero) ;
    'MESSAGE' 'Probleme dan le calcul du gaz brule' ;
    'ERREUR' 21 ;
 'FINSI' ;
 
*
**** mtot2 = weight of 1 mole of burned gas
*
 
 MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
   ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
 
*
**** Mass fractions of burned gas after the complete combustion
*
 
 YH22 = (XH22 '*' MH2) '/' MTOT2 ;
 YO22 = (XO22 '*' MO2) '/' MTOT2 ;
 YN22 = (XN22 '*' MN2) '/' MTOT2 ;
 YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
 
*
**** The conservative variables
*
 
 RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
        (('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
        (('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
        (('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
   ;
 
 RN = PINI '/' (RINI '*' TINI) ;
 
* Internal energy
 
 eini = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 0.0 ;
 Tfois = Tini ;
 
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
 
    acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
            &BL1))  ;
 
    eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
 
    Tfois = Tfois * tini ;
 
 'FIN' BL1 ;
 
 RETN = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
   ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
 GN = RN '*' WINI ;
 YN = YH21 '+' YO21 '+' YH2O1 ;
 RYN = RN '*' YN ;
 
*
**** K0/DX = chemical source term
*
*    i.e. dcsi/dt  = K0 '/' DX . {criterium function}
*         where
*         {criterium function} = 1.0 or 0.0
*
*         d(rho etot)/dt = -rho (dcsi/dt) (h_f '-' h_i)
*
*         h_f reference enthalpy (0K) of the burned gas
*         h_i reference enthalpy (0K) of the unburned gas
*
**** Mass fraction of H2 before and after the combustion
*
*    CSIMAX = 1.0 -> complete combustion
*
 
 YININ = yh21 ;
 YFINN = YININ '+' ((yh22 '-' yh21) '*' csimax) ;
 K0N = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 500 ;
 
*
**** Passive scalars  
*
 
 RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
               ('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
 
 'SI' VRAI ;
*   Test  
    VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
    PCEL = 'MAXIMUM' PINI ;
    TCEL = 'MAXIMUM' TINI ;
    'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
       'MESSAGE' 'Probleme calcul energie totale' ;
       'ERREUR' 21 ;
    'FINSI' ;
    'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
       'MESSAGE' 'Probleme calcul energie totale' ;
       'ERREUR' 21 ;
    'FINSI' ;
 'FINSI' ;
 
****************************************************************
****************************************************************
*****           Parameters for the computations           ******
****************************************************************
****************************************************************
 
* Upwind scheme
* METO = 'VLH' ;
 METO = 'SS' ;
 
* Iterations
* Final time
* Safety Factor fot the time step
* Second order reconstruction?
 
 NITER = 10000 ;
 TFINAL = 10.0D-3 ;
 SAFFAC = 0.7 ;
 LOGSO = VRAI ;
 
*
* EPS_CC linked to CREBCOM criterion
*
 
 EPS_CC = 0.8 ;
 
* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER'  METO  NITER  TFINAL  SAFFAC EPS_CC  LOGSO ;
 
****************************************************************
****************************************************************
*****           The computation                           ******
****************************************************************
****************************************************************
 
* Mesh dimension (in 2D!)
 
 DELTAXN = ('DOMA' $DOMTOT 'VOLUME') '**' (1. '/' 2.)  ;
 YN = RYN '/' RN ;
 
*
**** We impose we have combustion in the ignition region (DOM1)
*    We can use the operator 'FLAM', with epsilon=0 and
*    DELTAT >> DELTATC
*    where DELTATC is the characteristic time step linked to
*    the source term.
*    Since DELTAXN is constant, DELTATC does not change during the
*    computation.
*
 
 SN = RSN '/' RN ;
 K0N = 'EXCO' 'K0' SN ;
 YININ = 'EXCO' 'H2IN' SN 'H2' ;
 YFINN = 'EXCO' 'H2FI' SN 'H2' ;
 
 DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
 LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
 LCOEF = 'PROG' 1.0 0.5 -1.0 ;
 
 DELTARE DELTARY = 'FLAM' 'CREBCOM2'  $DOM1 PGAZ LMOT1 LCOEF
   ('REDU' RN ('DOMA' $DOM1 'CENTRE'))
   ('REDU' YN ('DOMA' $DOM1 'CENTRE'))
   ('REDU' YININ ('DOMA' $DOM1 'CENTRE'))
   ('REDU' YFINN ('DOMA' $DOM1 'CENTRE'))
   ('REDU' K0N ('DOMA' $DOM1 'CENTRE'))
   ('REDU' DELTAXN ('DOMA' $DOM1 'CENTRE'))
   0.0 (1D4 '*' DELTATC)  0.3 ;
 
 RYN = RYN '+' DELTARY ;
 YN = RYN '/' RN ;
 RETN = RETN '+' DELTARE ;
 
 
 RN0 = 'COPIER' RN ;
 GN0 = 'COPIER' GN ;
 RETN0 = 'COPIER' RETN ;
 RYN0 = 'COPIER' RYN ;
 RSN0 = 'COPIER' RSN ;
 
*
**** Parameter for the time loop 
*
 
* Names of the residuum components
 
 LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' 'H2IN'
   'H2FI' 'K0' ;
 
*
**** Geometrical coefficient to compute gradients
*
 
 GRADR CACCA COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
          ('MOTS' 'SCAL') RN ;
 
 GRADV CACCA COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
          ('MOTS' 'UX' 'UY') GN ;
 
 TPS = 0.0 ;
 
*
**** The chemistry time step
*
*
 
 DT_CHEM = SAFFAC '*' DELTATC ;
 
*
**** Temporal loop
*
 
 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
 'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
 'MESSAGE' ;
 
  VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
 
  TNM1 = 'COPIER' TN ;
 
 'TEMPS' 'ZERO' ;
 'REPETER' BL1 NITER ;
 
*
**** Primitive variables
*
 
    VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN
      TNM1  ;
 
    TNM1 = 'COPIER' TN ;
 
    'SI' LOGSO ;
 
       GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') RN  'GRADGEO' COEFSCAL ;
 
       GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'SCAL') PN  'GRADGEO' COEFSCAL ;
 
       GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY') VN   'GRADGEO'  COEFVECT ;
 
       GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'H2' 'O2' 'H2O') YN  'GRADGEO' COEFSCAL ;
 
       GRADS ALS = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
          ('MOTS' 'H2IN' 'H2FI' 'K0')  SN  'GRADGEO' COEFSCAL ;
 
       ROF VITF PF YF SF = 'PRET' 'PERFTEMP'  2 1
                            $DOMTOT  PGAZ
                            RN    GRADR  ALR
                            VN    GRADV  ALV
                            PN    GRADP  ALP
                            YN    GRADY  ALY
                            SN    GRADS  ALS ;
 
    'SINON' ;
 
       ROF VITF PF YF SF = 'PRET' 'PERFTEMP'  1 1 
                              $DOMTOT  PGAZ
                              RN   
                              VN   
                              PN   
                              YN   
                              SN    ;
 
    'FINSI' ;
 
    RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
         $DOMTOT PGAZ LISTINCO  ROF VITF PF YF SF ;
 
    DT_CON = SAFFAC '*' DELTAT ;     
*
**** The time step linked to tps
*
 
    DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
 
*
**** Total time step
*
 
    DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM) ;
 
*
**** Increment of the variables (convection)
*
 
    RESIDU = DTMIN '*' RESIDU ;
 
    DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
    DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
    DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
    DRYN  = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
            ('MOTS' 'H2' 'O2' 'H2O') ;
    DRSN = 'EXCO' ('MOTS' 'H2IN' 'H2FI' 'K0') RESIDU
           ('MOTS' 'H2IN' 'H2FI' 'K0') ;
 
 
    TPS  = TPS '+' DTMIN ;
    RN   = RN '+' DRN ;
    GN   = GN '+' DGN ;
    RETN = RETN '+' DRETN ;
    RYN = RYN '+' DRYN ;
    RSN = RSN '+' DRSN ;
    YN = RYN '/' RN ;
    SN = RSN '/' RN ;
 
*
**** Increment of the variables (chemical source)
*
 
    K0N = 'EXCO' 'K0' SN 'SCAL' ;
    YININ = 'EXCO' 'H2IN' SN 'H2' ;
    YFINN = 'EXCO' 'H2FI' SN 'H2' ;
 
    DRETN DRYN = 'FLAM' 'CREBCOM2' $DOMTOT PGAZ LMOT1 LCOEF RN YN
    YININ YFINN K0N DELTAXN EPS_CC DTMIN  0.3 ;
 
    RYN = RYN '+' DRYN ;
    RETN = RETN '+' DRETN ;
 
 
    'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
        'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
    'FINSI' ;
 
    'SI' (TPS > TFINAL) ;
       'QUITTER' BL1 ;
    'FINSI' ;
 
 'FIN' BL1 ;
 'TEMPS' ;
 
* 'OPTION' 'SAUVER' ('CHAINE' 'result.sauv_' RAF 'tps_' TPS ) ;
* 'SAUVER'  ;
 
****************************************************************
****************************************************************
*****           The post-treatment                        ******
****************************************************************
****************************************************************
 
* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;
 
 GRAPH = VRAI ;
 GRAPH = FAUX ;
 
* 
**** The mesh
*
 
 'SI' GRAPH ;
    'TRACER' ('DOMA' $DOMTOT 'MAILLAGE')
         'TITR' ('CHAINE' 'Maillage: ' ('NBEL' DOMTOT) ' elts');
 'FINSI' ;
 
*
**** Initial conditions
*
 
  MOD1 =  'MODELISER'  ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
 'SI' GRAPH ;
 
    VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP'
       PGAZ RN0 GN0 RETN0 RYN0 RSN0 ;
 
    CHM_RN  =  'KCHA'  $DOMTOT  'CHAM'  RN0  ;
    CHM_VN  =  'KCHA'  $DOMTOT  'CHAM'  VN  ;
    CHM_PN  =  'KCHA'  $DOMTOT  'CHAM'  PN  ;
    CHM_TN  =  'KCHA'  $DOMTOT  'CHAM'  TN  ;
    CHM_YN  =  'KCHA'  $DOMTOT  'CHAM'  YN  ;
    CHM_SN  =  'KCHA'  $DOMTOT  'CHAM'  SN  ;
 
    'TRAC' CHM_RN  MOD1 
       'TITR' ('CHAINE' 'rho at t=' 0.0)   ;
    'TRAC' CHM_VN  MOD1 
       'TITR' ('CHAINE' 'v at t= '  0.0)   ;
    'TRAC' CHM_PN  MOD1 
       'TITR' ('CHAINE' 'p at t= '  0.0)   ;
    'TRAC' CHM_TN  MOD1 
       'TITR' ('CHAINE' 'T at t= '  0.0)   ;
    'TRAC' CHM_YN  MOD1 
       'TITR' ('CHAINE' 'y at t= '  0.0) ;
    'TRAC' CHM_SN  MOD1 
       'TITR' ('CHAINE' 's at t= '  0.0) ;
 
 'FINSI' ;
 
*
**** Results
*
 
  VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
 
 'SI' GRAPH ;
 
    CHM_RN  =  'KCHA'  $DOMTOT  'CHAM'  RN  ;
    CHM_VN  =  'KCHA'  $DOMTOT  'CHAM'  VN  ;
    CHM_PN  =  'KCHA'  $DOMTOT  'CHAM'  PN  ;
    CHM_TN  =  'KCHA'  $DOMTOT  'CHAM'  TN  ;
    CHM_YN  =  'KCHA'  $DOMTOT  'CHAM'  YN  ;
    CHM_SN  =  'KCHA'  $DOMTOT  'CHAM'  SN  ;
 
    'TRAC' CHM_RN  MOD1 
       'TITR' ('CHAINE' 'rho at t=' TPS)   ;
    'TRAC' CHM_VN  MOD1 
       'TITR' ('CHAINE' 'v at t= '  TPS)   ;
    'TRAC' CHM_PN  MOD1 
       'TITR' ('CHAINE' 'p at t= '  TPS)   ;
    'TRAC' CHM_TN  MOD1 
       'TITR' ('CHAINE' 'T at t= '  TPS)   ;
    'TRAC' CHM_YN  MOD1 
       'TITR' ('CHAINE' 'y at t= '  TPS)   ;
    'TRAC' CHM_SN  MOD1 
       'TITR' ('CHAINE' 's at t= '  TPS)   ;
 
 'FINSI' ;
 
*
*** Evolution Objects
*
 
 SS = 'COORDONNEE' 1 LIGEVO;
 
 evxx = 'EVOL' 'CHPO' ss LIGEVO ;
 lxx = 'EXTRAIRE'  evxx 'ORDO' ;
 
 x0 = 'MINIMUM' lxx;
 x1 = 'MAXIMUM' lxx ;
 
 titolo  = 'CHAINE'  ' at t = ' TPS '  Methode = ' METO
   ' SF = ' SAFFAC ;
 
 VNN = 'EXCO' 'UX' VN ;
 VNT = 'EXCO' 'UY' VN ;
 
* rho
 
 evro = 'EVOL'  'CHPO' RN LIGEVO ;
 lro  = 'EXTRAIRE'  evro 'ORDO' ;
 evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro ;
 tro  = 'CHAINE' 'ro ' titolo ;
 
* u
 
 evu  = 'EVOL' 'CHPO' VNN LIGEVO ;
 lu   = 'EXTRAIRE'  evu 'ORDO' ;
 evu  = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
 tu   = 'CHAINE' 'u ' titolo ;
 
* yini
 
 YININ = 'EXCO' 'H2IN' SN  ;
 evyini  = 'EVOL' 'CHPO' YININ LIGEVO ;
 lyini   = 'EXTRAIRE'  evyini 'ORDO' ;
 evyini  = 'EVOL' 'MANU' 'x' lxx 'yini' lyini ;
 tyini   = 'CHAINE' 'yini ' titolo ;
 
* yfin
 
 YFINN = 'EXCO' 'H2FI' SN  ;
 evyfin  = 'EVOL' 'CHPO' YFINN LIGEVO ;
 lyfin   = 'EXTRAIRE'  evyfin 'ORDO' ;
 evyfin  = 'EVOL' 'MANU' 'x' lxx 'yfin' lyfin ;
 tyfin   = 'CHAINE' 'yfin ' titolo ;
 
* K0
 
 K0N = 'EXCO' 'K0' SN  ;
 evk0  = 'EVOL' 'CHPO' K0N LIGEVO ;
 lk0   = 'EXTRAIRE'  evk0 'ORDO' ;
 evk0  = 'EVOL' 'MANU' 'x' lxx 'k0' lk0 ;
 tk0   = 'CHAINE' 'k0 ' titolo ;
 
* csi
 
 CSIN = (('EXCO' 'H2' YN) '-' YININ) '/'
   (YFINN '-' YININ) ;
 evcsi  = 'EVOL' 'CHPO' CSIN LIGEVO ;
 lcsi   = 'EXTRAIRE'  evcsi 'ORDO' ;
 evcsi  = 'EVOL' 'MANU' 'x' lxx 'csi' lcsi ;
 tcsi   = 'CHAINE' 'csi ' titolo ;
 
* v
 
 evv  = 'EVOL' 'CHPO' VNT LIGEVO ;
 lv   = 'EXTRAIRE'  evv 'ORDO' ;
 evv  = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
 tv   = 'CHAINE' 'v ' titolo ;
 
 
* p
 
 evp  = 'EVOL' 'CHPO' PN LIGEVO ;
 lp   = 'EXTRAIRE'  evp 'ORDO' ;
 evp  = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
 tp   = 'CHAINE' 'p  ' titolo ;
 
 
* T
 
 evT  = 'EVOL' 'CHPO' TN LIGEVO ;
 lT   = 'EXTRAIRE'  evT 'ORDO' ;
 evT  = 'EVOL' 'MANU' 'x' lxx 'T' lT ;
 tT   = 'CHAINE' 'T  ' titolo ;
 
 
 'SI' GRAPH ;
    'DESSIN' evro 'TITRE' tro 'XBOR' x0 x1 'MIMA' ;
    'DESSIN' evp  'TITRE' tp 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evu 'TITRE'  tu 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evv 'TITRE'  tv 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evT  'TITRE' tT 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evcsi 'TITRE' tcsi 'XBOR' x0 x1 'MIMA' ;
    'DESSIN' evk0 'TITRE' tk0 'XBOR' x0 x1 'MIMA' ;
    'DESSIN' evyini 'TITRE' tyini 'XBOR' x0 x1 'MIMA' ;
    'DESSIN' evyfin 'TITRE' tyfin 'XBOR' x0 x1 'MIMA' ;
 'FINSI' ;
 
 
* Test
 
 PVN = 3.2E6 ;
 PCJ = 1.1E6 ;
 RHOVN = 3.7 ;
 RHOCJ = 1.6 ;
 
 PMAX = 'MAXIMUM' PN ;
 'SI' ((PMAX '>' PVN) 'OU' (PMAX '<' PCJ)) ;
    'MESSAGE' 'Pression = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 RHOMAX = 'MAXIMUM' RN ;
 'SI' ((RHOMAX '>' RHOVN) 'OU' (RHOMAX '<' RHOCJ)) ;
    'MESSAGE' 'Densite = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 'FIN' ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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