* fichier : rdem_surf1Daxi.dgibi
****************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la ****
**** combustion. ****
**** MODELE RDEM ****
**** ****
**** OPERATEURS 'PRIM', PRET, KONV ****
**** PROPAGATION D'UNE FLAMME DANS UN TUBE ****
**** We verify that S_flame = S_tube ****
**** ****
**** S. KUDRIAKOV SFME/LTMF ****
****************************************************************
*
'OPTION' 'ECHO' 1 'DIME' 2
'ELEM' 'QUA4' 'TRAC' 'X' ;
'OPTION' 'MODE' 'AXIS' ;
RAF = 1 ;
************************************************
* Constructing mesh *
************************************************
*---------------------------------
******** refinement level ********
*---------------------------------
EPSI = 0.154 '/' 10000.0D0 ;
*----------------------------------
** radius ************************
*----------------------------------
R = 0.15 '/' 2.0D0 ;
*-----length of the obstacle--------
LD = 0.03015D0 ;
*--------------------------------------------------------
**** number of points along the 'REST' of the radius
*--------------------------------------------------------
NRAD = 4 '*' RAF ;
******* base of the tube
*!!!!!!!!!!!!!!!!!!!!
A0 = 0.0 0.0 ;
A1 = R 0.0 ;
'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
A0 = EPSI 0.0 ;
'FINSI' ;
*********************************
DX = R '/' NRAD ;
****************************************************
****************************************************
A0A1 = 'DROIT' NRAD A0 A1 ;
*----------------------------------------
IBASE = A0A1 ;
* 'TRACER' IBASE ;
*-------------------------------------------------------
****** Length of the first section (without obstacles)
*-------------------------------------------------------
LFIRST = 3.0 ;
NFIRST = 'ENTIER' ((LFIRST '/' R) '*' NRAD) ;
*----------------------------
B0 = A0 'PLUS' (0.0 LFIRST) ;
B1 = A1 'PLUS' (0.0 LFIRST) ;
*--------------------------
A0B0 = 'DROIT' NFIRST A0 B0 ;
B0B1 = 'DROIT' NRAD B0 B1 ;
B1A1 = 'DROIT' NFIRST B1 A1 ;
*------------------------------
DOMTOT = 'DALLER' A0B0 B0B1 B1A1 ('INVERSE' A0A1) 'PLANE' ;
*********************************************************
***** Extracting ignition region ***********************
*********************************************************
LRIG = R ;
*----------------
YYAXI = 'COORDONNEE' 2 DOMTOT ;
YY_I = YYAXI 'POIN' 'EGINFE' LRIG ;
TIGN = DOMTOT 'ELEM' 'APPUYE' 'STRICTEMENT' YY_I ;
* 'TRACER' TIGN ;
'OUBLIER' YYAXI ;
'OUBLIER' YY_I ;
*************************************************
***** adapting the names ************************
*************************************************
DOM1 = TIGN ;
DOM2 = 'DIFF' DOMTOT DOM1 ;
DOMLIM = 'CONTOUR' DOMTOT ;
DOMTOT = DOMTOT 'COULEUR' 'VERT' ;
*********************************************
**** Extracting a line for posttreatment ****
*********************************************
'OPTION' 'ELEM' 'SEG2' ;
vev1 = (0.3 '*' DX) (0.25 '*' DX) ;
vev2 = (0.3 '*' DX) (LFIRST '-'(0.25 '*' DX)) ;
vev3 = (1.0 '*' DX) (LFIRST '-'(0.25 '*' DX)) ;
vev4 = (1.0 '*' DX) (0.25 '*' DX) ;
LIMCON = 'MANUEL' 'QUA4'
(A0 'PLUS' vev1)
(A0 'PLUS' vev2)
(A0 'PLUS' vev3)
(A0 'PLUS' vev4) ;
LIMCON = LIMCON 'COULEUR' 'ROUG' ;
***************************************************
****** Making models !!!!!!!!!!!!!!!!! ************
***************************************************
MDNS = 'EULER' ;
*---------------------------------
$DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
$DOMLIM = 'MODELISER' DOMLIM 'EULER' ;
$DOM1 = 'MODELISER' DOM1 'EULER' ;
$DOM2 = 'MODELISER' DOM2 'EULER' ;
*
TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
TDOMLIM = 'DOMA' $DOMLIM 'VF' ;
TDOM1 = 'DOMA' $DOM1 'VF' ;
TDOM2 = 'DOMA' $DOM2 'VF' ;
*
QDOMTOT = TDOMTOT . 'QUAF' ;
QDOMLIM = TDOMLIM . 'QUAF' ;
QDOM1 = TDOM1 . 'QUAF' ;
QDOM2 = TDOM2 . 'QUAF' ;
*
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOMLIM ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM1 ;
'ELIMINATION' QDOMTOT (0.001 '*' DX) QDOM2 ;
*--------------------------------------------
*--------------------------------------------
LIGEVO = (TDOMTOT . 'CENTRE') 'INCL' LIMCON ;
P1 = 'POIN' (COOR 1 LIGEVO) 'MINIMUM' ;
'ORDONNER' LIGEVO P1 ;
LIGEVO = 'QUELCONQUE' 'SEG2' LIGEVO 'COULEUR' 'VERT';
*************************************************
******* Volumes !!!!!!!!!!! *****************
*************************************************
***--- Volume of the total domain
VTOT = 'DOMA' $DOMTOT 'VOLUME' ;
CENDOM = 'DOMA' $DOMTOT 'CENTRE' ;
CHUN = ('MANUEL' 'CHPO' CENDOM 1 'SCAL' 1.0 'NATU' 'DISCRET') ;
VTOT = ('XTY' CHUN VTOT ('MOTS' 'SCAL')
('MOTS' 'SCAL')) ;
**************************************************
******* Pressure transducers *********************
**************************************************
PBAS = A0 ;
PTOP = B0 ;
*--------------------
PC1 = ('DOMA' $DOMTOT 'CENTRE') 'POIN' 'PROC' PBAS ;
PC2 = ('DOMA' $DOMTOT 'CENTRE') 'POIN' 'PROC' PTOP ;
PCTOT = PC1 'ET' PC2 ;
*-------------------------------------
GRAPH = VRAI ;
* GRAPH = FAUX ;
* Upwind scheme
METO = 'SS' ;
*****************************************
**** vitese fondamentale **************
*****************************************
K0 = 10.0D0 ;
***********************************************
**** Computational parameters *****************
***********************************************
NITER = 10000000 ;
TFINAL = 1.0D-3 ;
SAFFAC = 1.0 ;
LOGSO = VRAI ;
*************************************
**** table for pressure transducers *
*************************************
PCEN = 'TABLE' ;
PCEN.'TIME' = 'PROG' ;
PCEN.'PRESSION' = 'TABLE' ;
NPP = 'NBEL' PCTOT ;
***** two pressure transducers
PCEN.'PRESSION'. 1 = 'PROG' ;
PCEN.'PRESSION'. 2 = 'PROG' ;
*******************************************
*** table for velocity along centerline **
*******************************************
VCEN = 'TABLE' ;
VCEN . 'TIME' = 'PROG' ;
IPOST = 0 ;
**** 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 ;
*
* 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' ;
***************************************************
********* end of personal subroutines *************
***************************************************
****************************************************************
****************************************************************
***** Initial conditions and gas properties. *****
****************************************************************
****************************************************************
*
*************************************************
**** 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.
* meme que 'species' YH2_ini, YH2_fin, YO2_fin, YH2O_ini
PGAS . 'MASSFRA' = 'PROG' 1.03333E-02 0.964039E-11 1.48501E-01
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. * 1.00797E-3 ;
PGAS . 'O2 ' . 'W' = 2. * 15.9994E-3 ;
PGAS . 'H2O ' . 'W' = (PGAS . 'H2 ' . 'W' ) '+'
(0.5 * (PGAS . 'O2 ' . 'W' )) ;
PGAS . 'N2 ' . 'W' = 2 * 14.0067E-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 initial conditions *********************
*************************************************
*
eps = 1.0D-64 ;
*
T1 = 300.0D0 ;
* alph11 = 1.0 '-' 1.0D-5 ;
alph11 = 1.0D-5 ;
alph12 = 1.0 '-' 1.0D-5 ;
alph21 = 1.0 '-' alph11 ;
alph22 = 1.0 - alph12 ;
T2 = 1325.1D0 ;
un1 = 0.0 ;
un2 = 0.0 ;
ut1 = 0.0 ;
ut2 = 0.0 ;
pre1 = 1.D5 ;
pre2 = 1.D5 ;
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*** these lists are needed for flame evolution with time
*** along the centerline
DFLAM = 'PROG' LRIG ;
DFLAM1 = 'PROG' LRIG ;
TFLAM = 'PROG' 0.0 ;
TFLAM1 = 'PROG' 0.0 ;
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*
CHVN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un1
'UY' ut1) ;
CHVN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' un2
'UY' ut2) ;
*
CHTN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T1) ;
CHTN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' T2) ;
*
CHPN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre1) ;
CHPN2 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' pre2) ;
********************************************************
********************************************************
CHAL1 = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
alph11) '+'
('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
alph12) ;
CHAL2 = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL'
alph21) '+'
('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1 'SCAL'
alph22) ;
********************************************************
********************************************************
CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS
PGAS CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS 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 (TDOM1 . 'CENTRE'))) 'ABS' ;
ERRO = ERRO '/' (pre1) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (PRE2 '-' ('REDU' PN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
ERRO = ERRO '/' (pre2) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (T1 '-' ('REDU' TN1 (TDOM1 . 'CENTRE'))) 'ABS' ;
ERRO = ERRO '/' (T1) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (T2 '-' ('REDU' TN2 (TDOM2 . 'CENTRE'))) 'ABS' ;
ERRO = ERRO '/' (T2) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
**************************************************
**** initial conditions **************************
**************************************************
*** density
* DENSM = (CHAL1 * CHRN1) '+' (CHAL2 * CHRN2) ;
* CHM_RN = 'KCHA' $DOMTOT 'CHAM' DENSM ;
* 'TRACER' CHM_RN $DOMTOT
* 'TITR' ('CHAINE' 'RN at t=' 0.0) ;
**************************************************
**** temperat
* TEMPM = (CHAL1 * TN1) '+' (CHAL2 * TN2) ;
* CHM_TN = 'KCHA' $DOMTOT 'CHAM' TEMPM ;
* 'TRACER' CHM_TN $DOMTOT
* 'TITR' ('CHAINE' 'TN at t=' 0.0) ;
* 'OPTION' DONN 5 ;
*
**** Parameter for the time loop
*
* Names of the residuum components
*
LISTINCO = ('MOTS' 'ALF1' 'RN1' 'RNX1' 'RNY1' 'RET1'
'ALF2' 'RN2' 'RNX2' 'RNY2' 'RET2') ;
*
**** BC
*
KONLIM = 'DIFF' (TDOMLIM . 'CENTRE') (TDOMLIM . 'CENTRE') ;
CHPVID CHMVID = 'KOPS' 'MATRIK' ;
RESLIM = 'COPIER' CHPVID ;
************************************************************
* mass of burnt gase
************************************************************
MBURB = (CHAL2 * CHRN2) '*' ('DOMA' $DOMTOT 'VOLUME') ;
FMBURB = 'MAXIMUM' ('RESULT' MBURB) ;
TPSM1 = 0.0 ;
DENUNB = (CHAL1 * CHRN1) ;
LSFLAM = 'PROG' R ;
'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
LSFLAM = 'PROG' ((R '**' 2.0) '/' 2.0D0) ;
'FINSI' ;
* 'LISTE' FMBURB ;
* 'LISTE' ((CHAL2 * CHRN2)) ;
*
****************************************************************
****************************************************************
****************************************************************
***** The computation ******
****************************************************************
****************************************************************
AL1 = CHAL1 ;
AL2 = CHAL2 ;
ARN1 = CHAL1 * CHRN1 ;
ARN2 = CHAL2 * CHRN2 ;
AGN1 = CHAL1 * CHGN1 ;
AGN2 = CHAL2 * CHGN2 ;
AREN1 = CHAL1 * CHRET1 ;
AREN2 = CHAL2 * CHRET2 ;
*
AL10 = 'COPIER' AL1 ;
AL20 = 'COPIER' AL2 ;
ARN10 = 'COPIER' ARN1 ;
ARN20 = 'COPIER' ARN2 ;
AGN10 = 'COPIER' AGN1 ;
AGN20 = 'COPIER' AGN2 ;
AREN10 = 'COPIER' AREN1 ;
AREN20 = 'COPIER' AREN2 ;
TN10 = 'COPIER' TN1 ;
TN20 = 'COPIER' TN2 ;
*
**** Geometrical coefficient to compute grad(alpha)/|grad(alpha)|
*
EPSGR = 1.0D-6 ;
CHPL1 = CHPVID ;
CHPL2 = 'MANUEL' 'CHPO' (TDOMLIM . 'CENTRE') 2
'P1DX' 0.0 'P1DY' 0.0 ;
GRALP1 GRAL = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1 CHPL1 CHPL2 ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
*
**** Geometrical coefficient to compute gradients
*
GRADR0 ALR0 COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') AL10 ;
GRADR0 = GRADR0 * 0.0 ;
ALR0 = ALR0 * 0.0 ;
GRADV0 ALV0 COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
('MOTS' 'UX' 'UY') AGN10 ;
GRADV0 = GRADV0 * 0.0 ;
ALV0 = ALV0 * 0.0 ;
*
**** Temporal loop
*
TPS = 0.0 ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
'TEMPS' 'ZERO' ;
'REPETER' BLITER NITER ;
*
**** Primitive variables
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 TN1M TN2M EPS ;
TN1M = 'COPIER' TN1 ;
TN2M = 'COPIER' TN2 ;
**** Gradient of alpha
GRALP1 = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1 CHPL1 CHPL2 'GRADGEO' GRAL ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
'SI' LOGSO ;
GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL1 'GRADGEO' COEFSCAL ;
GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN1 'GRADGEO' COEFSCAL ;
GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN1 'GRADGEO' COEFSCAL ;
GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN1 'GRADGEO' COEFVECT ;
GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL2 'GRADGEO' COEFSCAL ;
GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN2 'GRADGEO' COEFSCAL ;
GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN2 'GRADGEO' COEFSCAL ;
GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN2 'GRADGEO' COEFVECT ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1 GRADA1 ALA1
AL2 GRADA2 ALA2
* AL1 GRADA1 ALR0
* AL2 GRADA2 ALR0
RN1 GRADR1 ALR1
RN2 GRADR2 ALR2
VN1 GRADV1 ALV1
VN2 GRADV2 ALV2
PN1 GRADP1 ALP1
PN2 GRADP2 ALP2 ;
'SINON' ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1 GRADR0 ALR0
AL2 GRADR0 ALR0
RN1 GRADR0 ALR0
RN2 GRADR0 ALR0
VN1 GRADV0 ALV0
VN2 GRADV0 ALV0
PN1 GRADR0 ALR0
PN2 GRADR0 ALR0 ;
'FINSI' ;
RESIDU DELTAT SURFF = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1 AL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
RESIDU = RESIDU '+' RESLIM ;
*********************************************************
*** axisymmetric case ****
* Contribution of pressure ****
* Warning: at first order, ALP initialised to 0 ****
*********************************************************
'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
******** first order for while
ALP1 = 0.0 '*' ALP1 ;
ALP2 = 0.0 '*' ALP2 ;
RESIP1 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
('MOTS' 'RN1' 'RNX1' 'RNY1' 'RET1')
(AL1 '*' PN1) GRADP1 ALP1 ;
RESIP2 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
('MOTS' 'RN2' 'RNX2' 'RNY2' 'RET2')
(AL2 '*' PN2) GRADP2 ALP2 ;
RESIDU = RESIDU 'ET' RESIP1 'ET' RESIP2 ;
'FINSI' ;
********************************************************
********************************************************
* 'REPETER' BL1 ('DIME' LISTINCO) ;
* mot1 = 'EXTRAIRE' LISTINCO &BL1 ;
* valre1 = 'MAXIMUM' ('EXCO' RESLIM mot1) 'ABS' ;
* tit1 = 'CHAINE' 'Component ' mot1 ', reference value ' valre1 ;
* evaa = 'EVOL' 'CHPO' RESIDU mot1 LIGEVO ;
* 'DESSIN' evaa 'TITRE' tit1 ;
* 'FIN' BL1 ;
* 'OPTION' DONN 5 ;
* Problem with RNY1, RNY2
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' 1.0D-6) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
*
**** Increment of the variables (predictor)
*
* RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1 AL2
* ARN1 ARN2 AGN1 AGN2 AREN1 AREN2 CHTN1 CHTN2 EPS ;
RESIDU = (0.5 * DTMIN) '*' RESIDU ;
DALP1 = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
DRN1 = 'EXCO' 'RN1' RESIDU 'SCAL' ;
DGN1 = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
DRET1 = 'EXCO' 'RET1' RESIDU 'SCAL' ;
DALP2 = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
DRN2 = 'EXCO' 'RN2' RESIDU 'SCAL' ;
DGN2 = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
DRET2 = 'EXCO' 'RET2' RESIDU 'SCAL' ;
AL1B = AL1 '+' DALP1 ;
ARN1B = ARN1 '+' DRN1 ;
AGN1B = AGN1 '+' DGN1 ;
AREN1B = AREN1 '+' DRET1 ;
AL2B = AL2 '+' DALP2 ;
ARN2B = ARN2 '+' DRN2 ;
AGN2B = AGN2 '+' DGN2 ;
AREN2B = AREN2 '+' DRET2 ;
*
**** Corrector
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS AL1B AL2B
ARN1B ARN2B AGN1B AGN2B AREN1B AREN2B TN1M TN2M EPS ;
**** Gradient of alpha
GRALP1 = 'PENT' $DOMTOT 'FACE'
'DIAMAN2' ('MOTS' 'SCAL') ('MOTS' 'P1DX' 'P1DY')
AL1B CHPL1 CHPL2 'GRADGEO' GRAL ;
* 'SI' VRAI ;
* V1 = 'VECTEUR'
* ('NOMC' ('MOTS' 'P1DX' 'P1DY') GRALP1 ('MOTS' 'UX' 'UY')) ;
* 'TRACER' DOMTOT V1 DOMLIM ;
* 'FINSI' ;
GRALP1 = GRALP1 '+' EPSGR ;
GRALP1 = GRALP1 '/' (('PSCAL' GRALP1 GRALP1 ('MOTS' 'P1DX' 'P1DY')
('MOTS' 'P1DX' 'P1DY')) '**' 0.5) ;
'SI' LOGSO ;
GRADA1 ALA1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL1B 'GRADGEO' COEFSCAL ;
GRADR1 ALR1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN1 'GRADGEO' COEFSCAL ;
GRADP1 ALP1 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN1 'GRADGEO' COEFSCAL ;
GRADV1 ALV1 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN1 'GRADGEO' COEFVECT ;
GRADA2 ALA2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') AL2B 'GRADGEO' COEFSCAL ;
GRADR2 ALR2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN2 'GRADGEO' COEFSCAL ;
GRADP2 ALP2 = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN2 'GRADGEO' COEFSCAL ;
GRADV2 ALV2 = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN2 'GRADGEO' COEFVECT ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1B GRADA1 ALA1
AL2B GRADA2 ALA2
RN1 GRADR1 ALR1
RN2 GRADR2 ALR2
VN1 GRADV1 ALV1
VN2 GRADV2 ALV2
PN1 GRADP1 ALP1
PN2 GRADP2 ALP2 ;
'SINON' ;
CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 =
'PRET' 'DEM' $DOMTOT
AL1B GRADR0 ALR0
AL2B GRADR0 ALR0
RN1 GRADR0 ALR0
RN2 GRADR0 ALR0
VN1 GRADV0 ALV0
VN2 GRADV0 ALV0
PN1 GRADR0 ALR0
PN2 GRADR0 ALR0 ;
'FINSI' ;
RESIDU DELTAT SURFP = 'KONV' 'VF' 'DEM' 'RESI' METO 'CONS'
$DOMTOT PGAS LISTINCO AL1B AL2B CHFAL1 CHFAL2 CHFRN1 CHFRN2
CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS KONLIM ;
RESIDU = RESIDU '+' RESLIM ;
*********************************************************
*** axisymmetric case ****
* Contribution of pressure ****
* Warning: at first order, ALP initialised to 0 ****
*********************************************************
'SI' ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
******** first order for while
ALP1 = 0.0 '*' ALP1 ;
ALP2 = 0.0 '*' ALP2 ;
RESIP1 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
('MOTS' 'RN1' 'RNX1' 'RNY1' 'RET1')
(AL1B '*' PN1) GRADP1 ALP1 ;
RESIP2 = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
('MOTS' 'RN2' 'RNX2' 'RNY2' 'RET2')
(AL2B '*' PN2) GRADP2 ALP2 ;
RESIDU = RESIDU 'ET' RESIP1 'ET' RESIP2 ;
'FINSI' ;
*
RESIDU = DTMIN '*' RESIDU ;
DALP1 = 'EXCO' 'ALF1' RESIDU 'SCAL' ;
DRN1 = 'EXCO' 'RN1' RESIDU 'SCAL' ;
DGN1 = 'EXCO' ('MOTS' 'RNX1' 'RNY1') RESIDU ('MOTS' 'UX' 'UY') ;
DRET1 = 'EXCO' 'RET1' RESIDU 'SCAL' ;
DALP2 = 'EXCO' 'ALF2' RESIDU 'SCAL' ;
DRN2 = 'EXCO' 'RN2' RESIDU 'SCAL' ;
DGN2 = 'EXCO' ('MOTS' 'RNX2' 'RNY2') RESIDU ('MOTS' 'UX' 'UY') ;
DRET2 = 'EXCO' 'RET2' RESIDU 'SCAL' ;
TPS = TPS '+' DTMIN ;
AL1 = AL1 '+' DALP1 ;
ARN1 = ARN1 '+' DRN1 ;
AGN1 = AGN1 '+' DGN1 ;
AREN1 = AREN1 '+' DRET1 ;
AL2 = AL2 '+' DALP2 ;
ARN2 = ARN2 '+' DRN2 ;
AGN2 = AGN2 '+' DGN2 ;
AREN2 = AREN2 '+' DRET2 ;
'SI' (((&BLITER '/' 10) '*' 10) 'EGA' &BLITER) ;
'MESSAGE' ('CHAINE' 'ITER =' &BLITER ' TPS =' TPS) ;
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*!!!!! HERE !!!!!!!
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IPOST = IPOST '+' 1 ;
* 'LISTE' DRN2 ;
************************************************
******* pressure transducers *******************
************************************************
PN = (AL1 '*' PN1) '+' (AL2 '*' PN2) ;
PPC = 'REDU' PN PCTOT ;
PCEN.'TIME' = PCEN.'TIME' 'ET' ('PROG' TPS) ;
'REPETER' BOUC NPP ;
PP1 = 'EXTRAIRE' PPC 'SCAL' ('POIN' PCTOT &BOUC);
PCEN.'PRESSION'. &BOUC = (PCEN.'PRESSION'. &BOUC) 'ET'
('PROG' PP1) ;
'FIN' BOUC ;
*************************************************
***** velocity along centerline *****************
*************************************************
VN = (AL1 '*' VN1) '+' (AL2 '*' VN2) ;
VYV = 'EXCO' 'UY' VN ;
EVVYV = ('EVOL' 'CHPO' VYV LIGEVO) 'COULEUR' 'ROUG' ;
* 'DESSIN' (EVVYV) 'GRILL' 'NCLK' ;
VCEN . IPOST = ('EVOL' 'CHPO' VYV LIGEVO) ;
* 'LISTE' (VCEN . IPOST) ;
*******************************************************
****** making variable ksi = \alpha_1 * (1-\alpha_1)
*******************************************************
EVALP = 'EVOL' 'CHPO' AL1 LIGEVO ;
* 'DESSIN' EVALP 'NCLK' ;
LCOOR = 'EXTRAIRE' EVALP 'ABSC' ;
LALP1 = 'EXTRAIRE' EVALP 'ORDO' ;
LUNI = 'PROG' ('DIME' LALP1) * 1.0 ;
VKSI = (-1.0D0 '*' LALP1) '+' LUNI ;
VKSI = LALP1 '*' VKSI ;
EVKSI = 'EVOL' 'MANUEL' LCOOR VKSI ;
************************************************
****** new way of computing flame location *****
************************************************
EVKIS = 'EVOL' 'MANUEL' VKSI LCOOR ;
EVKI = 'EXTRAIRE' EVKIS 'APRE' 0.23D0 ;
EVK = 'EXTRAIRE' EVKI 'ORDO' ;
'SI' (('DIME' EVK) '>' 1) ;
X1S = 'EXTRAIRE' EVK 1 ;
X2S = 'EXTRAIRE' EVK ('DIME' EVK) ;
DISTFL1 = (X1S '+' X2S) '/' 2.0D0 ;
DFLAM1 = DFLAM1 'ET' ('PROG' DISTFL1) ;
TFLAM1 = TFLAM1 'ET' ('PROG' TPS) ;
'FINSI' ;
**************************************************
********* the flame is where the maximum of ksi
**************************************************
OBJET2 DISTFL OBJET4 = 'MAXIMUM' EVKSI ;
DFLAM = DFLAM 'ET' ('PROG' DISTFL) ;
TFLAM = TFLAM 'ET' ('PROG' TPS) ;
EVDFLAM = 'EVOL' 'MANUEL' TFLAM DFLAM ;
*********************************************************
********* Flame surface storage *********************
*********************************************************
LSFLAM = LSFLAM 'ET' ('PROG' SURFF) ;
EVSFLAM = 'EVOL' 'MANUEL' TFLAM LSFLAM ;
* 'DESSIN' EVSFLAM 'GRILL' 'NCLK' ;
*-------------------------------------
* CHM_TN = 'KCHA' $DOMTOT 'CHAM' AL1 ;
* CNT1 = 'CONTOUR' ('DOMA' $DOMTOT 'MAILLAGE') ;
* 'TRACER' CHM_TN $DOMTOT CNT1
* 'TITR' ('CHAINE' 'TN at t=' TPS) 'NCLK' ;
* 'FINSI' ;
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
'FINSI' ;
'SI' (TPS > TFINAL) ;
'QUITTER' BLITER ;
'FINSI' ;
'FIN' BLITER ;
****************************************************
***** Test de non-regression ******************
****************************************************
***** exact flame surface
SEXACT = (R '**' 2.0) '/' 2.0D0 ;
***** maximum of computed flame surface
SMCOM = 'MAXIMUM' LSFLAM ;
* 'LISTE' SMCOM ;
****** relative error
ERREL = (SEXACT '-' SMCOM) '/' SEXACT ;
ERREL = ('ABS' ERREL) '*' 100.0D0 ;
* 'LISTE' ERREL ;
'SI' (ERREL '>' 0.02 ) ;
'MESSAGE' 'Problem 1';
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer