Numérotation des lignes :

* fichier :  prim_ther_dem3D.dgibi*************************************************************************************************************************************************************************************************************** FV "Cell-Centred Formulation" for the solution of ******** the Euler equations for a thermally perfect gas.  ******** Discrete Equation Method for the propagation of   ******** infinitely thin flames in initially homogeneous   ******** media.                                            ******** PRIM operator                                     ******** It computes primitive variables from conservative ******** ones.                                             ******** 3D case                                           ********                                                   ******** A. BECCANTINI DEN/DM2S/SFME/LTMF  DECEMBER 2010   *****************************************************************   alpha1 = 1 -> unburnt phase*   alpha2 = 1 -> burnt phase*   In exact algebra, alpha2 = 1 - alpha1*   alpha2 can be seen as a progress variable* ******************************** PRODECURE PRIMCONS ********************************* 'DEBPROC' PRIMCONS ; 'ARGUMENT' PGAS*TABLE TN1*'CHPOINT ' TN2*'CHPOINT '   PN1*'CHPOINT ' PN2*'CHPOINT ' VN1*'CHPOINT ' VN2*'CHPOINT ' ;** ETHER = int_0^T cv(T') dT'           T <  TMAX*       = int_0^TMAX cv(T') dT' '+'*         cv(TMAX)                     T >= TMAX ESP1 = 'EXTRAIRE' (PGAS . 'SPECIES') 1 ;* DY1 = y_i - y_f for the species 1  DY1 = (('EXTRAIRE' (PGAS . 'MASSFRA') 1) '-'   ('EXTRAIRE' (PGAS . 'MASSFRA') 2)) ; COEF1 = ('EXTRAIRE' (PGAS . 'CHEMCOEF') 1) '*'    (PGAS . ESP1 . 'W') ; YFINPH1 = 1.0 ; YFINPH2 = 1.0 ; 'SI' (COEF1 > 0) ;    YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;    YPH1 = YPH2 '+' DY1 ; 'SINON' ;    YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;    YPH2 = YPH1 '-' DY1 ; 'FINSI' ; YFINPH1 = YFINPH1 '-' YPH1 ; YFINPH2 = YFINPH2 '-' YPH2 ; PRYPH1 = 'PROG' YPH1 ; PRYPH2 = 'PROG' YPH2 ; 'REPETER' BLESP (('DIME' (PGAS . 'SPECIES')) '-' 2) ;    ESP = 'EXTRAIRE' (PGAS . 'SPECIES') (&BLESP '+' 1)  ;    COEF = ('EXTRAIRE' (PGAS . 'CHEMCOEF') (&BLESP '+' 1))        '*' (PGAS . ESP . 'W') ;    DY = (DY1 * (COEF '/' COEF1)) ;    'SI' (COEF > 0) ;       YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;       YPH1 = YPH2 '+' DY ;    'SINON' ;       YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;       YPH2 = YPH1 '-' DY ;    'FINSI' ;    PRYPH1 = PRYPH1 'ET' ('PROG' YPH1) ;    PRYPH2 = PRYPH2 'ET' ('PROG' YPH2) ;    YFINPH1 = YFINPH1 '-' YPH1 ;    YFINPH2 = YFINPH2 '-' YPH2 ; 'FIN' BLESP ; PRYPH1 = PRYPH1 'ET' ('PROG' YFINPH1) ;  PRYPH2 = PRYPH2 'ET' ('PROG' YFINPH2) ; 'LISTE' PRYPH1 ; 'LISTE' PRYPH2 ; * TMAX = (PGAS . 'TMAX') ;* TCAL1 = MIN TN1, TMAX  TCAL1 = 0.5D0 '*' ((TMAX '+' TN1) '-' ('ABS' (TN1 '-' TMAX))) ; DTN1 = TN1 '-' TCAL1 ;* TCAL1 = MIN TN1, TMAX  TCAL2 = 0.5D0 '*' ((TMAX '+' TN2) '-' ('ABS' (TN2 '-' TMAX))) ; DTN2 = TN2 '-' TCAL2 ;** Internal energy (J/kg in SI)* ETHER1 = 0.0 ; CV1 = 0.0 ; ETHER2 = 0.0 ; CV2 = 0.0 ; FUNTN1 = 1.0 ; FUNTN2 = 1.0 ; 'REPETER' BLPO ((PGAS . 'NORD') '+' 1) ;    'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;       ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;       YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;       YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;       AA = 'EXTRAIRE' (PGAS . ESP . 'A') &BLPO ;       DCV1 = (AA * YCEL1 * FUNTN1) ;       DCV2 = (AA * YCEL2 * FUNTN2) ;       CV1 = CV1 '+' DCV1 ;       CV2 = CV2 '+' DCV2 ;       ETHER1 = ETHER1 '+' (DCV1 * TCAL1 '/' (&BLPO)) ;       ETHER2 = ETHER2 '+' (DCV2 * TCAL2 '/' (&BLPO)) ;    'FIN' BLESP ;    FUNTN1 = FUNTN1 '*' TCAL1 ;    FUNTN2 = FUNTN2 '*' TCAL2 ; 'FIN' BLPO ; ETHER1 = ETHER1 '+' (CV1 '*' DTN1) ; ETHER2 = ETHER2 '+' (CV2 '*' DTN2) ;** Formation energy/enthalpy (J/kg in SI) and gas constant (J/kg/K)* EFORM1 = 0.0 ; EFORM2 = 0.0 ; RGAS1 = 0.0 ; RGAS2 = 0.0 ; 'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;    ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;    YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;    YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;    EFORM1 = EFORM1 '+' (YCEL1 * (PGAS . ESP . 'H0K')) ;    EFORM2 = EFORM2 '+' (YCEL2 * (PGAS . ESP . 'H0K')) ;    RGAS1 = RGAS1 '+' (YCEL1 *  (PGAS . 'RUNIV') '/'       (PGAS . ESP . 'W')) ;    RGAS2 = RGAS2 '+' (YCEL2 *  (PGAS . 'RUNIV') '/'       (PGAS . ESP . 'W')) ; 'FIN' BLESP ; 'LISTE' RGAS1 ; 'LISTE' RGAS2 ;** Computation of the conservative variables* RN1 = PN1 '/' (RGAS1 '*' TN1) ; RN2 = PN2 '/' (RGAS2 '*' TN2) ; GN1 = RN1 * VN1 ; GN2 = RN2 * VN2 ; LVEL = 'MOTS' 'UX' 'UY' 'UZ' ; ECIN1 = 0.5D0 '*' ('PSCAL' VN1 VN1 LVEL LVEL) ; ECIN2 = 0.5D0 '*' ('PSCAL' VN2 VN2 LVEL LVEL) ; RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ;  RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ;* 'RESPRO' RN1 RN2 GN1 GN2 RETN1 RETN2 ; 'FINPROC' ;* ************************************ FIN PRODECURE PRIMCONS ************************************** 'OPTION'  'DIME'  3 ; 'OPTION'  'ELEM' 'CUB8' ; 'OPTION'  'ECHO'  1 ; 'OPTION'  'TRAC' 'X';***** GRAPH* GRAPH = FAUX ;* GRAPH = VRAI ;**** We consider a mixture of H2, O2, H2O, N2****************************************************** The table for the properties of the gas ****************************************************** PGAS = 'TABLE' ;***** Order of the polynomial order for cv = cv(T) *    For T > TMAX, cv(T) = cv(Tmax)*  PGAS . 'TMAX' = 6000.0 ; PGAS . 'NORD' = 4 ;***** Species involved in the mixture (before or after*    the chemical reaction)* PGAS . 'SPECIES' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' 'N2  ' ;****** Coefficient of the chemical reaction.*    Note that for the first species this coefficient should be positive*    Normal, we take it equal to 1.**    H2 '+' 0.5 O2 ---> H2O* PGAS . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 0.0 ;* **** Mass fraction of the first species before and after the combustion*    Final mass fractions of the species with positive coefficients.*    Final mass fractions of the species with non-positive coefficient.*    The mass fraction of the last species is not given.* PGAS . 'MASSFRA' = 'PROG' 0.285219E-01  0.964039E-11  0.765104E-10                           0.127442E-10 ;***** Coef with the gas properties* PGAS .  'H2  ' = 'TABLE'  ; PGAS .  'H2O ' = 'TABLE'  ; PGAS .  'N2  ' = 'TABLE'  ; PGAS .  'O2  ' = 'TABLE'  ;***** Runiv (J/mole/K)* PGAS . 'RUNIV' = 8.31441 ;***** W (kg/mole). Gas constant (J/kg/K = Runiv/W)* PGAS .  'H2  ' . 'W' = 2.016E-3 ; PGAS .  'O2  ' . 'W' = 31.999E-3 ; PGAS .  'H2O ' . 'W' = 18.0155E-3 ; PGAS .  'N2  ' . 'W' = 28.013E-3 ;***** Polynomial coefficients * PGAS .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836                               -2.37281455E-07 1.84701105E-11 ; PGAS .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05                              -1.82753232E-08 2.44485692E-12 ; PGAS .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05                              8.78233606E-09 -3.05514485E-13 ; PGAS .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865                              2.33636971E-08 -1.53304905E-12;***** Formation enthalpies (energies) at 0K (J/Kg)* PGAS .  'H2  ' . 'H0K' = -4.195D6 ; PGAS .  'H2O ' . 'H0K' = -1.395D7 ; PGAS .  'N2  ' . 'H0K' = -2.953D5 ; PGAS .  'O2  ' . 'H0K' = -2.634D5 ;  ***************************************************** The mesh                                *****************************************************  A1 = 0.0D0 0.0D0 0.0D0 ; A2 = 0.0D0 1.0D0 0.0D0 ; A3 = 1.0D0 0.0D0 0.0D0 ; A4 = 1.0D0 1.0D0 0.0D0 ; A5 = 2.0D0 0.0D0 0.0D0 ; A6 = 2.0D0 1.0D0 0.0D0 ;  MAIL1  = 'MANUEL' 'POI1' A1 ; MAIL1 = (MAIL1 'ET' A2) 'ET' A3 ; MAIL2  = 'MANUEL' 'POI1' A4 ; MAIL2 = (MAIL2 'ET' A5) 'ET' A6 ; MAIL = MAIL1 'ET' MAIL2 ;  'SI' GRAPH ;   'TRACER' MAIL 'TITRE' 'Domaine' ; 'FINSI' ;  ***************************************************** Initial condition on primitive variables ****************************************************  TN1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 293.16) '+'   ('MANUEL' 'CHPO' MAIL2 1 'SCAL' 283.16) ; TN2 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 2800.15) '+'    ('MANUEL' 'CHPO' MAIL2 1 'SCAL' 2910.15) ; PN1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 1.013E5) '+'   ('MANUEL' 'CHPO' MAIL2 1 'SCAL' 1.011E5) ; PN2 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 1.023E5) '+'   ('MANUEL' 'CHPO' MAIL2 1 'SCAL' 1.033E5) ; VN1 = ('MANUEL' 'CHPO' MAIL1 3 'UX' 122.    'UY' 101. 'UZ' 0.0 ) '+'     ('MANUEL' 'CHPO' MAIL2 3 'UX' 120.     'UY' 102. 'UZ' 0.0 ) ; VN2 =  ('MANUEL' 'CHPO' MAIL1 3 'UX' 150.    'UY' 160. 'UZ' 0.0 ) '+'    ('MANUEL' 'CHPO' MAIL2 3 'UX' 140.     'UY' 160. 'UZ' 0.0 )  ; ALPHA1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 0.9) '+'   ('MANUEL' 'CHPO' MAIL2 1 'SCAL' 0.1)  ; ALPHA2 = -1 '*' (ALPHA1 '-' 1.0) ;** Computation of the conservative variables* RN1 RN2 GN1 GN2 RETN1 RETN2 = PRIMCONS   PGAS TN1 TN2 PN1 PN2 VN1 VN2 ; * ARN1 = ALPHA1 * RN1 ; ARN2 = ALPHA2 * RN2 ; AGN1 = ALPHA1 * GN1 ; AGN2 = ALPHA2 * GN2 ; ARETN1 = ALPHA1 * RETN1 ; ARETN2 = ALPHA2 * RETN2 ;** K0 = 22.5 ; EPS = 1.0D-16 ;* TG1 = TN1 '*' 10. ; TG2 = TN2 '/' 10. ; R1 R2 V1 V2 P1 P2 T1 T2       = 'PRIM' 'DEM' PGAS ALPHA1  ALPHA2          ARN1 ARN2 AGN1 AGN2           ARETN1 ARETN2 TG1 TG2 EPS ;* ERRO1 = 'MAXIMUM' ((TN1 '-' T1) '/' TN1) 'ABS' ; ERRO2 = 'MAXIMUM' ((TN2 '-' T2) '/' TN2) 'ABS' ; ERRO3 = 'MAXIMUM' ((PN2 '-' P2) '/' PN2) 'ABS' ; ERRO4 = 'MAXIMUM' ((PN1 '-' P1) '/' PN1) 'ABS' ; ERRO5 = 'MAXIMUM' ((RN1 '-' R1) '/' RN1) 'ABS' ; ERRO6 = 'MAXIMUM' ((RN2 '-' R2) '/' RN2) 'ABS' ; ERRO = 'MAXIMUM' ('PROG' ERRO1 ERRO2 ERRO3 ERRO4   ERRO5 ERRO6) ; 'SI' (ERRO > 1.0D-12) ;     'ERREUR' 5 ; 'FINSI' ;**** We consider a mixture of H2, O2, H2O****************************************************** The table for the properties of the gas ****************************************************** PGAS2 = 'TABLE' ;***** Order of the polynomial order for cv = cv(T) *    For T > TMAX, cv(T) = cv(Tmax)*  PGAS2 . 'TMAX' = 6000.0 ; PGAS2 . 'NORD' = 4 ;***** Species involved in the mixture (before or after*    the chemical reaction)* PGAS2 . 'SPECIES' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;****** Coefficient of the chemical reaction.*    Note that for the first species this coefficient should be positive*    Normal, we take it equal to 1.**    H2 '+' 0.5 O2 ---> H2O* PGAS2 . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 ;***** Coef with the gas properties* PGAS2 .  'H2  ' = 'TABLE'  ; PGAS2 .  'H2O ' = 'TABLE'  ; PGAS2 .  'O2  ' = 'TABLE'  ;***** Runiv (J/mole/K)* PGAS2 . 'RUNIV' = 8.31441 ;***** W (kg/mole). Gas constant (J/kg/K = Runiv/W)* PGAS2 .  'H2  ' . 'W' =  2. * 1.00797E-3 ; PGAS2 .  'O2  ' . 'W' =  2. * 15.9994E-3 ; PGAS2 .  'H2O ' . 'W' = (PGAS2 .  'H2  ' . 'W' ) '+'    (0.5 * (PGAS2 .  'O2  ' . 'W' )) ;* **** Mass fraction of the first species before and after the combustion*    Final mass fractions of the species with positive coefficients.*    Final mass fractions of the species with non-positive coefficient.*    The mass fraction of the last species is not given.* NH2 = 1.; NO2 = 0.4999999999 * NH2 ; NH2O = 1.0D-12 ; MH2 = NH2 *  (PGAS2 .  'H2  ' . 'W') ; MO2 = NO2 *  (PGAS2 .  'O2  ' . 'W') ; MH2O = NH2O *  (PGAS2 .  'H2O ' . 'W') ; MTOT = MH2 '+' MO2 '+' MH2O ; YH2 = MH2 '/' MTOT ; YO2 = MO2 '/' MTOT ; YH2O = MH2O '/' MTOT ; SI FAUX ;***** Complete combustion*    SI (NH2 > (2.0 * NO2)) ;       YO2F = 1.0D-12 ;       YH2F = YH2 '-' ((YO2 '-' YO2F) * ((PGAS2 . 'H2  ' . 'W') '/'          (0.5 * (PGAS2 . 'O2  ' . 'W')))) ;    'SINON' ;       YH2F = 1.0D-12 ;       YO2F = YO2 '-' (0.5 * (YH2 '-' YH2F) * ((PGAS2 . 'O2  ' . 'W')          '/'(PGAS2 . 'H2  ' . 'W'))) ;    'FINSI' ; 'SINON' ;    YH2F = 1. '/' 3. * YH2 ;    YO2F = YO2 '-' (0.5 * (YH2 '-' YH2F) * ((PGAS2 . 'O2  ' . 'W') '/'       (PGAS2 . 'H2  ' . 'W'))) ; 'FINSI' ; YH2OF = YH2O '+' (YH2 '-' YH2F) '+' (YO2 '-' YO2F) ; ERRO = ((YH2F '+' YO2F '+' YH2OF) '-' 1.0) ; SI (ERRO > 1.0D-12) ;    'ERREUR' 21 ; 'FINSI' ;*  PGAS2 . 'MASSFRA' = 'PROG' YH2 YH2F YO2F ;***** Polynomial coefficients * PGAS2 .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836                               -2.37281455E-07 1.84701105E-11 ; PGAS2 .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05                              -1.82753232E-08 2.44485692E-12 ; PGAS2 .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865                              2.33636971E-08 -1.53304905E-12;***** Formation enthalpies (energies) at 0K (J/Kg)* PGAS2 .  'H2  ' . 'H0K' = -4.195D6 ; PGAS2 .  'H2O ' . 'H0K' = -1.395D7 ; PGAS2 .  'O2  ' . 'H0K' = -2.634D5 ;******************************************************* The initial conditions *********************************************************************** eps = 1.0D-64 ; K0 = 3000. ;*  T1   = 473. ;* alph11 = 1.0 '-' 1.0D-5 ; alph11 = 1.0 '-' 5.0D-3 ;* alph12 = 1.0D-5 ; alph12 = 5.0D-3 ; alph21 = 1.0 '-' alph11 ; alph22 = 1.0 - alph12 ; T2   = 4000. ; un1  = 1900. ; un2  = 1200. ; ut1  = 1221. ; ut2  = 1432. ; utt1 = 1231. ; utt2 = 1439. ; pre1 = 6.909D5 ; pre2 = 6.909D5 ;* CHVN1  = ('MANUEL' 'CHPO' MAIL 3 'UX' un1   'UY' ut1 'UZ' utt1) ; CHVN2  = ('MANUEL' 'CHPO' MAIL 3 'UX' un2   'UY' ut2 'UZ' utt2) ;* CHTN1  = ('MANUEL' 'CHPO' MAIL 1 'SCAL' T1) ; CHTN2  = ('MANUEL' 'CHPO' MAIL 1 'SCAL' T2) ;* CHPN1  = ('MANUEL' 'CHPO' MAIL 1 'SCAL' pre1) ; CHPN2  = ('MANUEL' 'CHPO' MAIL 1 'SCAL' pre2) ;* DOMAL = MAIL1 ; DOMRES = 'DIFF' MAIL DOMAL ; CHAL1  = ('MANUEL' 'CHPO' DOMRES 1 'SCAL'   alph11) '+'   ('MANUEL' 'CHPO' DOMAL 1 'SCAL'   alph12) ; CHAL2  = ('MANUEL' 'CHPO' DOMRES 1 'SCAL'   alph21) '+'   ('MANUEL' 'CHPO' DOMAL 1 'SCAL'   alph22) ;* CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS   PGAS2 CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;* RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS2 CHAL1 CHAL2   (CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1)   (CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2)    CHTN1 CHTN2 EPS ; TN1M = COPIER TN1 ; TN2M = COPIER TN2 ;* ERRO = 'MAXIMUM' (PRE1 '-' ('REDU' PN1 DOMAL)) 'ABS' ; ERRO = ERRO '/' (pre1) ; 'SI' (ERRO > 1.0D-6) ;     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);     'ERREUR' 21 ; 'FINSI' ;  ERRO = 'MAXIMUM' (PRE2 '-' ('REDU' PN2 DOMRES)) 'ABS' ; ERRO = ERRO '/' (pre2) ; 'SI' (ERRO > 1.0D-6) ;     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);     'ERREUR' 21 ; 'FINSI' ; ERRO = 'MAXIMUM' (T1 '-' ('REDU' TN1 DOMAL)) 'ABS' ; ERRO = ERRO '/' (T1) ; 'SI' (ERRO > 1.0D-6) ;     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);     'ERREUR' 21 ; 'FINSI' ; ERRO = 'MAXIMUM' (T2 '-' ('REDU' TN2 DOMRES)) 'ABS' ; ERRO = ERRO '/' (T2) ; 'SI' (ERRO > 1.0D-6) ;     'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);     'ERREUR' 21 ; 'FINSI' ;* 'FIN' ;

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