Télécharger flamarrh.dgibi

Retour à la liste

* fichier :  flamarrh.dgibi
************************************************************************
************************************************************************
**********************************************************************
****                                                              ****
**** OPERATEUR FLAM                                               ****
**** Combustion Hydrogene-Oxygene, cinetique de type Arrhenius    ****
****                                                              ****
**** A. BECCANTINI DRN/DMT/SEMT/LTMF        FEVRIER 1999          ****
**********************************************************************
 
'OPTION'  'DIME' 2 ;
'OPTION'  'ELEM' QUA4 ;
'OPTION'  'TRAC' 'X' ;
'OPTION'  'ECHO'  0  ;
 
*
*** GRAPH
*
 
* GRAPH = VRAI ;
GRAPH = FAUX ;
 
***********************************************************************
***************** CALCUL DU TAUX DE REACTION  *************************
***********************************************************************
 
'DEBPROC' TAUX ;
*
**** La meme chose que FLAM ARRHENIUS
*    
 
ARGUMENT TS*CHPOINT A*FLOTTANT b*FLOTTANT con*FLOTTANT Ta*FLOTTANT
H0H2*FLOTTANT H0O2*FLOTTANT H0H2O*FLOTTANT DT*FLOTTANT
RYH2*CHPOINT RYO2*CHPOINT RYH2O*CHPOINT T*CHPOINT ;
 
 
MH2 = 2.016D-3 ;
MO2 = 31.998D-3 ;
MH2O = 18.015D-3 ;
 
MAIL = 'EXTRAIRE' TS 'MAILLAGE' ;
NPOI = 'NBEL' MAIL ;
 
LIRYH2  = 'PROG' NPOI * 0.0D0 ;
LIRYO2  = 'PROG' NPOI * 0.0D0 ;
LIRYH2O = 'PROG' NPOI * 0.0D0 ;
LIDE  = 'PROG' NPOI * 0.0D0 ; 
 
'REPETER' BL1 NPOI ;
    P1 = 'POIN' &BL1 MAIL ;
    RYH2I = 'EXTRAIRE' RYH2 'H2  ' P1 ;
    RYO2I = 'EXTRAIRE' RYO2 'O2  ' P1 ;
    RYH2OI = 'EXTRAIRE'RYH2O  'H2O ' P1 ;
    TS0 = 'EXTRAIRE' TS 'SCAL' P1 ; 
    T0  = 'EXTRAIRE' T  'SCAL' P1 ; 
    XH2 = RYH2I '/' MH2 ;
    XO2 = RYO2I '/' MO2 ;
    'SI' (T0 > TS0) ;
       'SI' (XH2 > (2.0D0 * XO2)) ;
*       H2 majoritaire
           LAMBDA = EXP(-1.0D0*TA/T0);
           LAMBDA = LAMBDA * (T0 ** (-1.0D0*B));
           LAMBDA = LAMBDA * (RYH2I ** CON);
           LAMBDA = LAMBDA * MO2 * A;
           NRYO2I = RYO2I / (1.0D0 + (LAMBDA * DT));
           DRYO2 = NRYO2I '-' RYO2I;
           DMOL = DRYO2 / MO2;
           DRYH2 = 2.0D0 * DMOL * MH2;
           DRYH2O = -2.0D0 * DMOL * MH2O;
           NRYH2I = RYH2I '+' DRYH2 ;
           NRYH2OI = RYH2OI '+' DRYH2O ;
           'MESSAGE' ('CHAINE' 'POINT ' &BL1) ;
           'MESSAGE' 'O2  minoritaire' ;
           'MESSAGE' ('CHAINE' 'DTc = ' (1.0 '/' LAMBDA) '  DT = ' DT) ;
       'SINON' ;
*       O2 majoritaire
           LAMBDA = EXP(-1.0D0*TA/T0);
           LAMBDA = LAMBDA * (T0 ** (-1.0D0*B));
           LAMBDA = LAMBDA * (RYH2I ** (CON-1.0D0)) * RYO2I;
           LAMBDA = 2.0D0 * LAMBDA * MH2 * A; 
           NRYH2I = RYH2I / (1.0D0 + (LAMBDA * DT)) ;
           DRYH2 =  NRYH2I '-' RYH2I;
           DMOL = DRYH2 / MH2 / 2.0D0;
           DRYO2 = DMOL * MO2 ;
           DRYH2O = -2.0D0 * DMOL * MH2O;
           NRYO2I = RYO2I '+' DRYO2 ;
           NRYH2OI = RYH2OI '+' DRYH2O ;
           'MESSAGE' ('CHAINE' 'POINT ' &BL1) ;
           'MESSAGE' 'H2  minoritaire ' ;
           'MESSAGE' ('CHAINE' 'DTc = ' (1.0 '/' LAMBDA) '  DT = ' DT) ;
        'FINSI' ;
        DEI = (DRYH2 * H0H2) + (DRYO2 * H0O2) + (DRYH2O * H0H2O) ;
        DEI = -1.0D0 * DEI ;
     'SINON' ;
        NRYO2I = RYO2I ;
        NRYH2I = RYH2I ;
        NRYH2OI = RYH2OI ;
        DEI = 0.0D0 ;
     'FINSI' ;
     'REMPLACER' LIRYH2  &BL1 NRYH2I ;
     'REMPLACER' LIRYO2  &BL1 NRYO2I ;
     'REMPLACER' LIRYH2O &BL1 NRYH2OI ;
     'REMPLACER' LIDE    &BL1 DEI ;
'FIN' BL1 ;
 
NRYH2  = 'MANUEL' 'CHPO' MAIL 1 'H2  ' LIRYH2  ;
NRYO2  = 'MANUEL' 'CHPO' MAIL 1 'O2  ' LIRYO2  ;
NRYH2O = 'MANUEL' 'CHPO' MAIL 1 'H2O ' LIRYH2O ;
DE = 'MANUEL' 'CHPO' MAIL 1 'SCAL' LIDE ;
 
 
'FINPROC' NRYH2 NRYO2 NRYH2O DE ;
 
***********************
**** LA TABLE PGAZ ****
***********************
 
PGAZ = 'TABLE' ;
 
*
**** Ordre des polynoms
*
 
PGAZ . 'NORD' = 4 ;
 
*
**** Especes qui sont dans les equations d'Euler
*
 
PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;
 
*
**** Espece qui n'y est pas
*
 
 
PGAZ . 'ESPNEULE' = 'N2  ';
 
*
 
PGAZ .  'H2  ' = 'TABLE'  ;
PGAZ .  'H2O ' = 'TABLE'  ;
PGAZ .  'N2  ' = 'TABLE'  ;
PGAZ .  'O2  ' = 'TABLE'  ;
 
*
**** R (J/Kg/K)
*
 
PGAZ .  'H2  ' . 'R' = 4130.0 ;
PGAZ .  'H2O ' . 'R' = 461.4 ;
PGAZ .  'N2  ' . 'R' = 296.8 ;
PGAZ .  'O2  ' . 'R' = 259.8 ;
 
 
*
**** Regressions polynomials pour le cvs
*
 
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;
 
 
 
'SI' FAUX ;
 
*
**** Calcul des enthalpies de formation a 0K
*
*    A 298.15K
*
*    HOH2 = 0 ;
*    HOO2 = 0 ;
*    HON2 = 0 ;
*    HOH2O = -1.3428D7 J/kg ;
*
 
   RH2  = PGAZ .  'H2  ' . 'R' ;
   RO2  = PGAZ .  'O2  ' . 'R' ;
   RN2  = PGAZ .  'N2  ' . 'R' ;
   RH2O = PGAZ .  'H2O ' . 'R' ;
 
   A0H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 1 ;
   A1H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 2 ;
   A2H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 3 ;
   A3H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 4 ;
   A4H2 = 'EXTRAIRE' (PGAZ . 'H2  ' . 'A') 5 ;
 
   A0H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 1 ;
   A1H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 2 ;
   A2H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 3 ;
   A3H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 4 ;
   A4H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 5 ;
 
   A0N2 = 'EXTRAIRE' (PGAZ . 'N2  ' . 'A') 1 ;
   A1N2 = 'EXTRAIRE' (PGAZ . 'N2  ' . 'A') 2 ;
   A2N2 = 'EXTRAIRE' (PGAZ . 'N2  ' . 'A') 3 ;
   A3N2 = 'EXTRAIRE' (PGAZ . 'N2  ' . 'A') 4 ;
   A4N2 = 'EXTRAIRE' (PGAZ . 'N2  ' . 'A') 5 ;
 
   A0O2 = 'EXTRAIRE' (PGAZ . 'O2  ' . 'A') 1 ;
   A1O2 = 'EXTRAIRE' (PGAZ . 'O2  ' . 'A') 2 ;
   A2O2 = 'EXTRAIRE' (PGAZ . 'O2  ' . 'A') 3 ;
   A3O2 = 'EXTRAIRE' (PGAZ . 'O2  ' . 'A') 4 ;
   A4O2 = 'EXTRAIRE' (PGAZ . 'O2  ' . 'A') 5 ;
 
*  On calcule \int_0^298.15 cv(y) dy pour chaque especes
 
   T = 298.15  ;
   T2 =  T * T ;
   T3 = T2 * T ;
   T4 = T3 * T ;
   T5 = T4 * T ;
 
   ETHERH2 = ((A0H2 * T) '+' ((A1H2 '/' 2.0) * T2) '+'
             ((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
             ((A4H2 '/' 5.0) * T5)) ; 
 
   ETHERO2 = ((A0O2 * T) '+' ((A1O2 '/' 2.0) * T2) '+'
             ((A2O2 '/' 3.0) * T3) '+' ((A3O2 '/' 4.0) * T4) '+'
             ((A4O2 '/' 5.0) * T5)) ; 
 
   ETHERH2O= ((A0H2O * T) '+' ((A1H2O '/' 2.0) * T2) '+'
             ((A2H2O '/' 3.0) * T3) '+' ((A3H2O '/' 4.0) * T4) '+'
             ((A4H2O '/' 5.0) * T5)) ; 
 
   ETHERN2 = ((A0N2 * T) '+' ((A1N2 '/' 2.0) * T2) '+'
             ((A2N2 '/' 3.0) * T3) '+' ((A3N2 '/' 4.0) * T4) '+'
             ((A4N2 '/' 5.0) * T5)) ;
*
 
   ENTHH2 = -1.0D0 * (ETHERH2 '+' (RH2 '*' 298.15)) ;
 
   ENTHO2 = -1.0D0 * (ETHERO2 '+' (RO2 '*' 298.15)) ;
 
   ENTHN2 = -1.0D0 * (ETHERN2 '+' (RN2 '*' 298.15)) ;
 
   ENTHH2O = -1.0D0 * (1.3428D7  '+' ETHERH2O '+' (RH2O '*' 298.15)) ;
 
'FINSI' ;
 
 
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
*
 
PGAZ .  'H2  ' . 'H0K' = -4.194D6 ;
PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;
 
 
*************************** 
***** DOMAINE SPATIAL  ****
***************************
 
*
**** Trois points
*
 
A1 = 0.0D0 0.0D0 ;
DOM1 = 'MANUEL' 'POI1' A1 ;
A2 = 1.0D0 0.0D0 ;
DOM2 = 'MANUEL' 'POI1' A2 ;
A3 = 2.0D0 0.0D0 ;
DOM3 = 'MANUEL' 'POI1' A3 ;
 
DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
 
*
**** FLAM, 'ARRHENIU'
*
*    2 H2 + O2 -> 2 H2O
*
*    w = A H(T - Ts) RH2^c RO2 * T^(-b) exp(-Ta/T)
*                      
*    ou
*
*    w = taux de reaction MOLAIRE pour la reaction 
*        chimie, qui est suppose' etre elementaire
*                      
*    H(T - Ts) est la fontion d'Heavyside
*  
*    N.B. L'energie thermique de chaque epsece est definie par
*
*          e = \int_{0}^T cv(y) dy
*
 
Ta = 8310.0;
A = 1.1725D14;
b = 0.91;
con = 2.0;
DT = 0.1;
TS = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' 1000 'NATURE' 'DISCRET' ;      
 
********************************************
************* Donnes initiales *************
********************************************
 
RH2  = PGAZ .  'H2  ' . 'R' ;
RO2  = PGAZ .  'O2  ' . 'R' ;
RN2  = PGAZ .  'N2  ' . 'R' ;
RH2O = PGAZ .  'H2O ' . 'R' ;
 
H0H2  = PGAZ .  'H2  ' . 'H0K' ;
H0O2  = PGAZ .  'O2  ' . 'H0K' ;
H0H2O = PGAZ .  'H2O ' . 'H0K' ;
* H0N2  = PGAZ .  'N2  ' . 'H0K' ; 
 
MH2 = 2.016D-3 ;
MO2 = 31.998D-3 ;
MH2O = 18.015D-3 ;
MN2 = 28.013D-3 ;
 
XH2 = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 0.3 0.25 0.3) ;
 
XO2 = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 0.147 0.147 0.147) ;
 
XH2O = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 1.0D-6 1.0D-6 1.0D-6) ;
 
XN2 = ('MANUEL' 'CHPO' DOMTOT 1 'SCAL' 1.0 'NATURE' 'DISCRET' ) '-'
      (XH2 '+' XO2 '+' XH2O) ;
 
* Mass des especes dans une mole
 
MH2SMOL = XH2 * MH2 ;
MO2SMOL = XO2 * MO2 ;
MH2OSMOL = XH2O * MH2O ;
MN2SMOL = XN2 * MN2 ;
MTOT = MH2SMOL '+' MO2SMOL '+' MH2OSMOL '+' MN2SMOL ;
 
YH2  = MH2SMOL  '/' MTOT ;
YO2  = MO2SMOL  '/' MTOT ;
YH2O = MH2OSMOL '/' MTOT ;
YN2  = MN2SMOL  '/' MTOT ;
 
RTOT = (YH2 * RH2) '+' (YO2 * RO2) '+' (YH2O * RH2O) '+' (YN2 * RN2) ;
 
* Pression et temperature
 
P =  'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 2.83D6 2.0D6 2.5D6) ; 
 
T =  'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 1527D0 2500 1650 ) ; 
 
RO = P '/' (RTOT '*' T) ;
 
 
RYH2 = 'NOMC'  'H2  ' (RO * YH2) ;
 
RYO2 = 'NOMC'  'O2  ' (RO * YO2) ;
 
RYH2O = 'NOMC'  'H2O ' (RO * YH2O) ;
 
 
* Dans le cas de detonation dans un tube ou C1 = zone von Neumann et
* C3 = reste du tube, on a
* u ~ 1000 m/s
* gamma ~ 1.25
* donc
* ason ~ 875 m/s -> (u + cson) ~ 1875 m/s
*
* DT ~ DX '/' (u + cson)
*
* DX ~ 0.1 m  ->  DT ~ 5.00000E-05
*
 
DT = 0.5D-5 ;
 
 
NRYH2 NRYO2 NRYH2O DE = 'FLAM' 'ARRHENIU' TS A b con Ta H0H2 H0O2
                        H0H2O DT RYH2 RYO2 RYH2O T ;
 
*
**** Conservativite'
*
 
ERRRO = (('EXCO' 'H2  ' NRYH2) '+' ('EXCO' 'O2  'NRYO2)  '+'
         ('EXCO' 'H2O ' NRYH2O)) '-'
         (('EXCO' 'H2  ' RYH2) '+' ('EXCO' 'O2  ' RYO2)  '+'
         ('EXCO' 'H2O ' RYH2O));
 
ERRRO = ERRRO '/' ('MAXIMUM' RO) ;         
 
'SI' (('MAXIMUM' ERRRO 'ABS') > 1.0D-10) ;
    'ERREUR' 5 ;
'FINSI' ;    
 
NRYH21 NRYO21 NRYH2O1 DE1 = TAUX TS A b con Ta H0H2 H0O2
                        H0H2O DT RYH2 RYO2 RYH2O T ;
 
 
ERRH2  = ('MAXIMUM' (NRYH2  '-' NRYH21)   'ABS')
         '/' ('MAXIMUM' NRYH2  ) ;
 
ERRO2  = 'MAXIMUM' (NRYO2  '-' NRYO21)   'ABS' 
         '/' ('MAXIMUM' NRYO2  ) ;
 
ERRH2O = 'MAXIMUM' (NRYH2O  '-' NRYH2O1) 'ABS'
         '/' ('MAXIMUM' NRYH2O  ) ;
 
 
ERRDE = 'MAXIMUM' (DE '-' DE1) 'ABS' '/'
        ('MAXIMUM' DE) ;
 
 
'SI' (('MAXIMUM' ('PROG' ERRH2 ERRO2 ERRH2O ERRDE)) > 1.0D-10) ;
     'ERREUR' 5;
'FINSI' ;     
 
*
**** Option HEAVYSIDE: les especes minoritaires brule completement
*
 
NRYH2H NRYO2H NRYH2OH DEH = 'FLAM' 'HEAVYSID'  TS  H0H2 H0O2
                             H0H2O  RYH2 RYO2 RYH2O T ;
 
*
**** Conservativite'
*
 
ERRRO = (('EXCO' 'H2  ' NRYH2H) '+' ('EXCO' 'O2  'NRYO2H)  '+'
         ('EXCO' 'H2O ' NRYH2OH)) '-'
         (('EXCO' 'H2  ' RYH2) '+' ('EXCO' 'O2  ' RYO2)  '+'
         ('EXCO' 'H2O ' RYH2O));
 
ERRRO = ERRRO '/' ('MAXIMUM' RO) ;         
 
'SI' (('MAXIMUM' ERRRO 'ABS') > 1.0D-10) ;
    'ERREUR' 5 ;
'FINSI' ;    
 
'FIN' ;
 
 
 
 
 
 
 

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