* fichier : rut_tg_2.dgibi
************************************************************************
**************** CAS-TEST RUT *************************************
************************************************************************
* Date de creation : 16 decembre 2002
* Auteur : Alexandre Bleyer
* Commentaire : ce cas-test est de coupe en plusieurs
* parties : procedures, definition du maillage, de la table PGAZ
* et definition des conditons initiales
************************************************************************
* Liste des procedures : CONS, DETO, EXEX, PERSO, PREPADET, VALP
************************************************************************
'OPTI' 'ECHO' 0;
*$$$$ CONS
***************************************************************
***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
***************************************************************
'DEBPROC' CONS RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT'
YN*'CHPOINT' $DM*'MMODEL' PGAZ*'TABLE' ;
*
A0 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A1 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A2 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A3 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
A4 = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
RTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
YTOT = 'KCHT' $DM 'SCAL' 'CENTRE' 0.0 ;
'REPETER' BL1 ('DIME' (PGAZ . 'ESPEULE')) ;
NOMCOM = 'EXTRAIRE' &BL1 (PGAZ . 'ESPEULE');
YCEL = 'EXCO' NOMCOM YN ;
RYCEL = 'NOMC' (YCEL * RN) NOMCOM 'NATU' 'DISCRET' ;
YTOT = YTOT '+' YCEL ;
RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') * YCEL) ;
A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
'SI' ('EGA' &BL1 1) ;
RYN = 'KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL ;
'SINON' ;
RYN = RYN 'ET'
('KCHT' $DM 'SCAL' 'CENTRE' 'COMP' NOMCOM RYCEL) ;
'FINSI' ;
'FIN' BL1 ;
YCEL = ('KCHT' $DM 'SCAL' 'CENTRE' 1.0) '-' YTOT ;
NOMCOM = PGAZ . 'ESPNEULE' ;
A0 = A0 '+' (('EXTRAIRE' 1 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A1 = A1 '+' (('EXTRAIRE' 2 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A2 = A2 '+' (('EXTRAIRE' 3 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A3 = A3 '+' (('EXTRAIRE' 4 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
A4 = A4 '+' (('EXTRAIRE' 5 (PGAZ . NOMCOM . 'A')) '*' YCEL) ;
RTOT = RTOT '+' ((PGAZ . NOMCOM . 'R') '*' YCEL);
TN = PN '/' (RN '*' RTOT) ;
T2 = TN '*' TN ;
T3 = T2 '*' TN ;
T4 = T3 '*' TN ;
T5 = T4 '*' TN ;
ETHER = (A0 * TN) '+' ((A1 '/' 2.0) * T2) '+'
((A2 '/' 3.0) * T3) '+' ((A3 '/' 4.0) * T4) '+'
((A4 '/' 5.0) * T5) ;
GN = RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX ' 'UY ')
('MOTS' 'UX ' 'UY ') ;
GN = 'KCHT' $DM 'VECT' 'CENTRE' GN ;
ECIN = 0.5D0 '*' ('PSCAL' GN VN ('MOTS' 'UX ' 'UY ')
('MOTS' 'UX ' 'UY ')) ;
REN = (RN '*' ETHER) '+' ECIN ;
REN = 'KCHT' $DM 'SCAL' 'CENTRE' REN ;
'FINPROC' GN REN RYN ;
*$$$$ DETO
'DEBPROC' DETO ;
'ARGUMENT' PGAZ*'TABLE' TN*'CHPOINT' RYN*'CHPOINT' DTI*'FLOTTANT' ;
************************
**** PROCEDURE DETO ****
************************
*
* dRYO2/dT = MO2 * A * (T ** -B) EXP( - TA / T)
* (RYH2 ** C) * (RYO2)
*
**** Les parametres physiques
*
* Ts = temperature de seuil
*
TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' 600. ;
TA = 8310.;
A = 1.1725D14;
B = 0.91;
CON = 2;
H0H2 = PGAZ . 'H2 ' . 'H0K' ;
H0O2 = PGAZ . 'O2 ' . 'H0K' ;
H0H2O = PGAZ . 'H2O ' . 'H0K' ;
RYH2 = 'EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET' ;
RYO2 = 'EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET' ;
RYH2O = 'EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET' ;
RYH2 RYO2 RYH2O DELTAE = 'FLAM' 'ARRHENIU' TS A B CON TA
H0H2 H0O2 H0H2O DTI RYH2 RYO2 RYH2O
TN ;
RY = RYH2 'ET' RYO2 'ET' RYH2O ;
'FINPROC' RY DELTAE ;
****************************
**** FIN PROCEDURE DETO ****
****************************
*$$$$ EXEX
'DEBPROC' EXEX RV*'TABLE' ;
****************************************************************
****************************************************************
***** Parameters for the computations ******
****************************************************************
****************************************************************
* Upwind scheme
METO = RV.'CI'.'METO' ;
* Iterations
* Final time
* Safety Factor for the time step
* Second order reconstruction?
TPS = RV . 'PASDETPS' . 'TINI' ;
NTPS = 'DIME' RV.'PASDETPS'.'LTPS' ;
'SI' (NTPS '>' 1) ;
TPS = 'EXTR' (RV.'PASDETPS'.'LTPS') NTPS ;
'FINSI';
NITER = RV . 'PASDETPS' . 'NITMA' ;
TFINAL = RV . 'PASDETPS' . 'TFINA' ;
SAFFAC = RV . 'PASDETPS' . 'SAFFAC' ;
LOGSO = 'EGA' (RV.'ORDREESP') 2 ;
GEO = RV . 'GEO' ;
INCO = RV . 'INCO' ;
PGAZ = RV . 'PGAZ' ;
$DOMTOT= GEO . 'MODTOT' ;
ZERO = 1.0d-9 ;
*
**** Chemical time step
*
* We compute d(RHY2)/dT at the von Neumann state
*
* Molar reaction rate
*
TVN = 2000 ;
TA = 8310.;
A = 1.1725D14;
B = 0.91;
CON = 2;
H0H2 = PGAZ . 'H2 ' . 'H0K' ;
H0O2 = PGAZ . 'O2 ' . 'H0K' ;
H0H2O = PGAZ . 'H2O ' . 'H0K' ;
RYN = INCO . 'RYN' ;
GN = INCO . 'GN' ;
RN = INCO . 'RN' ;
RETN = INCO . 'RETN' ;
IGNI1 = 2 ;
MOD_I = 'CHAIN' 'MOD' IGNI1;
$DOM2 = GEO . MOD_I;
RYNRED = 'REDU' RYN ('DOMA' $DOM2 'CENTRE') ;
RYH2 = 'MAXIMUM' ('EXCO' 'H2' RYNRED 'SCAL') ;
RYO2 = 'MAXIMUM' ('EXCO' 'O2' RYNRED 'SCAL') ;
OMEGA = A '*' (RYH2 '**' CON) '*' RYO2 '*' (TVN '**' (-1.0 '*' B))
'*' ('EXP' (-1.0 '*' TA '*' (TVN '**' -1))) ;
* OMEGA in (mol/m3/s)
* OMEGAH2 in (kg/m3/s)
OMEGAH2 = OMEGA '*' 2 '*' 2E-3 ;
DELTATC = RYH2 '/' OMEGAH2 ;
* NB DELTATC << DELTAT_CON
*
**** Zone d activation = etat AICC
*
IGNI1 = RV . 'MODETO' .'IGNIREG' ;
MOD_I = 'CHAIN' 'MOD' IGNI1;
$DOM1 = GEO . MOD_I;
VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
RYH2 = ('EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET') ;
RYO2 = ('EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET') ;
RYH2O = ('EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET') ;
TN = TN '-' ('REDU' TN ('DOMA' $DOM1 'CENTRE')) '+'
('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' TVN) ;
TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' (0.9 '*' TVN)
;
RYH2 RYO2 RYH2O DELTARE = 'FLAM' 'ARRHENIU' TS A B CON TA
H0H2 H0O2 H0H2O (1.0D5 '*' DELTATC)
RYH2 RYO2 RYH2O TN ;
RYN = RYH2 'ET' RYO2 'ET' RYH2O ;
RETN = RETN '+' DELTARE ;
RN0 = 'COPIER' RN ;
GN0 = 'COPIER' GN ;
RETN0 = 'COPIER' RETN ;
RYN0 = 'COPIER' RYN ;
*
**** Parameter for the time loop
*
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ;
*
**** 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 ;
*
**** The chemistry time step is too little
* We take DT_CHEM > DELTATC
*
*
DT_CHEM = 50 '*' 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 GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
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 ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 2 1
$DOMTOT PGAZ
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
YN GRADY ALY ;
'SINON' ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 1 1
$DOMTOT PGAZ
RN
VN
PN
YN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
$DOMTOT PGAZ LISTINCO ROF VITF PF YF ;
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
*
**** Total time step
*
DTIMP = RV.'PASDETPS'.'DTIMP';
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM DTIMP) ;
RV.'PASDETPS'.'DT' = RV.'PASDETPS'.'DT' 'ET' ('PROG' DTMIN) ;
*
**** 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') ;
TPS = TPS '+' DTMIN ;
RV.'PASDETPS'.'LTPS' = RV.'PASDETPS'.'LTPS' 'ET'
('PROG' TPS) ;
RV.'PASDETPS'.'TPS' = TPS ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
RYN = RYN '+' DRYN ;
*
**** Increment of the variables (chemical source)
*
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
TNM1 ;
RYN DRETN = DETO PGAZ TN RYN DTMIN ;
RETN = RETN '+' DRETN ;
RV.'INCO'.'PN' = PN ;
RV.'INCO'.'TN' = TN ;
RV.'INCO'.'RN' = RN ;
RV.'INCO'.'VN' = VN ;
RV.'INCO'.'YN' = YN ;
RV.'INCO'.'GAMMAN' = GAMMAN ;
RV.'INCO'.'GN' = GN ;
RV.'INCO'.'RETN' = RETN ;
RV.'INCO'.'RYN' = RYN ;
'SI' ('EXISTE' RV 'TABC' );
TABC = VALP RV ;
RV . 'TABC' = TABC ;
'FINSI' ;
'SI' ('EXISTE' RV 'PRCPERSO' );
FREQ0 = RV.'PRCPERSO'.'ARG3' ;
'SI' (((&BL1 '/' FREQ0) '*' FREQ0) 'EGA' &BL1) ;
('TEXT' (RV.'PRCPERSO'.'NOMPROC')) RV (RV.'PRCPERSO'.'ARG1')
(RV.'PRCPERSO'.'ARG2') ;
'FINSI' ;
'FINSI' ;
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS '>EG' TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
RV.'PASDETPS'.'TINI' = TPS ;
'FINPROC' RV;
*$$$$ PERSO
'DEBP' PERSO RV*'TABLE' $DOMTRA*'MMODEL' TRACE*'ENTIER';
INCO = RV . 'INCO' ;
GEO = RV . 'GEO' ;
TPS = RV . 'PASDETPS' . 'TPS' ;
*
DOMT1 = 'DOMA' $DOMTRA 'MAILLAGE' ;
X0 Y0 = 'COORD' DOMT1 ;
Xmin = 'MINI' X0 ;
Xmax = 'MAXI' X0 ;
Ymin = 'MINI' Y0 ;
Ymax = 'MAXI' Y0 ;
DTX = (Xmax - Xmin)/5. ;
DTY = (Ymax - Ymin)/3. ;
Pt1 = (Xmin - DTX) (Ymin - DTY) ;
Pt2 = (Xmax + DTX) (Ymin - DTY) ;
Pt3 = (Xmax + DTX) (Ymax + DTY) ;
Pt4 = (Xmin - DTX) (Ymax + DTY) ;
CONT1 = CONT DOMT1 ;
TRACHA = FAUX ;
TRACAP = FAUX ;
'SI' (TRACE 'EGA' 1);
TRACHA = VRAI ;
'FINSI';
'SI' (TRACE 'EGA' 2);
TRACAP = VRAI ;
'FINSI';
'SI' (TRACE 'EGA' 3);
TRACHA = VRAI ;
TRACAP = VRAI ;
'FINSI';
'SI' TRACHA ;
AB = KCHA $DOMTRA INCO.'PN' 'CHAM' ;
TIT1 = 'CHAIN' 'PRESSION' TPS ;
LIST1= 'PROG' 1.E5 'PAS' 0.4E5 13.E5 ;
TRAC AB $DOMTRA CONT1 LIST1 'TITR' TIT1 ;
*
AB = KCHA $DOMTRA INCO.'TN' 'CHAM' ;
TIT1 = 'CHAIN' 'TEMPERATURE' TPS ;
LIST1= 'PROG' 363. 'PAS' 50. 3000. ;
TRAC AB $DOMTRA CONT1 LIST1 'TITR' TIT1 ;
*
*AB = KCHA $DOMT1 INCO.'RN' 'CHAM' ;
*TIT1 = 'CHAIN' 'DENSITE' TPS ;
*TRAC AB $DOMT1 cont1 'TITR' TIT1 ;
*
*YN = INCO.'YN' ;
*AA = 'EXCO' 'H2 ' YN ;
*AB = KCHA $DOMT1 AA 'CHAM' ;
*TIT1 = 'CHAIN' 'YH2' TPS ;
*TRAC AB $DOMT1 cont1 'TITR' TIT1 ;
*
*AA = 'EXCO' 'O2 ' YN ;
*AB = KCHA $DOMT1 AA 'CHAM' ;
*TIT1 = 'CHAIN' 'YO2' TPS ;
*TRAC AB $DOMT1 cont1 'TITR' TIT1 ;
*
*AA = 'EXCO' 'H2O ' YN ;
*AB = KCHA $DOMT1 AA 'CHAM' ;
*TIT1 = 'CHAIN' 'YH2O' TPS ;
*TRAC AB $DOMT1 cont1 'TITR' TIT1 ;
*
V2 = KCHT $DOMTRA VECT CENTRE INCO.'VN' ;
V1 = VECT V2 0.003 ;
TIT1 = 'CHAIN' 'VITESSE' TPS ;
TRAC V1 (CONT1 'ET' Pt1 'ET' Pt2 'ET' Pt3 'ET' Pt4) 'TITR' TIT1 ;
*
'FINSI';
'SI' TRACAP ;
TABC = RV . 'TABC' ;
NBEVOL = 'NBEL' TABC.'LPOINTS';
NBCOMP = 'DIME' TABC.'COMP' ;
LOG1 = 'EXISTE' TABC 'NOMC' ;
*
*- Initialisations
*
MARQLIST = 'MOTS' 'CARR' 'LOSA' 'TRIA' 'TRIB' 'CROI' 'PLUS' 'ETOI' ;
COULLIST = 'MOTS' 'TURQ' 'BLEU' 'ROUG' 'VERT' ;
LEGEND = 'TABLE' ;
TEVOL = 'TABLE' ;
'REPETER' BLO10 NBCOMP ;
LCOMP1 = 'EXTR' (TABC.'COMP') &BLO10 ;
IPOS2 = 1 ;
IPOS0 = 0 ;
I100 = 0 ;
'REPETER' BLO100 NBEVOL ;
I100 = I100 + 1 ;
'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' NBEVOL )) ;
SUPP1 = 'EXTR' (TABC.'NOMC') I100;
'SINON' ;
SUPP1 = 'CHAINE' 'CAPT' I100;
'FINSI';
IPOS1 = I100 - (7 * (I100 / 7)) ;
'SI' (IPOS1 EGA 0) ;
IPOS0 = 1 ;
IPOS1 = 7 ;
'FINSI' ;
IPOS3 = (IPOS1 + 1) - (2 * ((IPOS1 + 1) / 2)) ;
'SI' (IPOS3 EGA 0) ;
IPOS2 = IPOS2 + 1 ;
'SI' (IPOS2 EGA 5) ;
IPOS2 = 1 ;
IPOS3 = IPOS1 ;
'FINSI' ;
'FINSI' ;
IMARQ1 = 'EXTR' MARQLIST IPOS1 ;
ICOUL1 = 'EXTR' COULLIST IPOS2 ;
LEGEN1 = 'CHAIN' 'MARQ ' IMARQ1 ;
LEGEN2 = 'CHAIN' SUPP1 ;
*
EVOL1 = EVOL ICOUL1 'MANU' 'TIME' TABC.'LTPS' TABC.SUPP1.LCOMP1 ;
'SI' (IPOS1 EGA 1) ;
TEVOL . LCOMP1 = EVOL1 ;
LEGEND . LCOMP1 = TABLE ;
LEGEND . LCOMP1 . 'TITRE' = TABLE ;
LEGEND . LCOMP1 . 1 = 'MOT' LEGEN1 ;
LEGEND . LCOMP1 . 'TITRE' . 1 = 'MOT' LEGEN2 ;
'SINON' ;
TEVOL . LCOMP1 = (TEVOL . LCOMP1) ET EVOL1 ;
NDIML1 = 'DIME' (LEGEND . LCOMP1 . 'TITRE') ;
LEGEND . LCOMP1 . (NDIML1 + 1) = 'MOT' LEGEN1 ;
LEGEND . LCOMP1 . 'TITRE' . (NDIML1 + 1) = 'MOT' LEGEN2 ;
'FINSI' ;
'SI' ((IPOS0 EGA 1) 'OU' (I100 'EGA' NBEVOL)) ;
IPOS0 = 0 ;
'DESS' (TEVOL . LCOMP1) 'LEGE' 'MIMA' 'GRIL' (LEGEND .LCOMP1)
'TITR' ('CHAIN' 'Tracé de ' LCOMP1) 'DATE' ;
'FINSI' ;
*
'FIN' BLO100 ;
'FIN' BLO10 ;
'FINSI' ;
*
FINP ;
*** Debut proc perso
*$$$$ PREPADET
'DEBPROC' PREPADET CONDINI*'TABLE' PGAZ*'TABLE' TABC/'TABLE';
ERROR = 0 ;
*** DEFINITION DU PARAMETRE EPSI POUR ELIM DU MAILLAGE
'SI' ('EXISTE' CONDINI 'EPSI');
EPSI = CONDINI . 'EPSI' ;
'SINON';
EPSI = 0.01 ;
'FINSI';
*** DEFINITION DU PARAMETRE EPSI POUR ELIM DU MAILLAGE
'SI' ('NON' ('EXISTE' CONDINI 'GRAPH'));
CONDINI . 'GRAPH' = FAUX ;
'FINSI';
*** INITILISATION DE TABLES TEMPORAIRES OU NON
INI = 'TABLE' ;
INCO = 'TABLE' 'INCO' ;
GEO = 'TABLE' 'GEO';
FRACY = 'TABLE' ;
*** CREATION DES OBJETS MODELES
*** ils sont sauvegardes dans la table GEO de RV
MAILTOT = CONDINI . 'MAILTOT' ;
$MAILTOT = 'MODE' MAILTOT 'EULER' ;
MMAILTOT = 'DOMA' $MAILTOT 'QUAF' ;
'ELIM' MMAILTOT EPSI ;
*** Parametres verifiant si l utilisateur a defini
*** - soit des fractions massiques
*** - soit des fractions volumiques
*** Pour le calcul on n utilise uniquement des fractions
*** massiques.
*** Si l utilisateur a defini des fractions volumiques
*** il faut se ramener a des fractions massiques (pre-calculs)
LOGX = FAUX ;
LOGY = FAUX ;
'REPETER' BCL1 CONDINI . 'NBDOM';
MODE_I = 'CHAIN' 'MOD' &BCL1;
MQUA_I = 'CHAIN' 'MQUA' &BCL1;
DOM_I = 'CHAIN' 'DOM' &BCL1;
M_I = 'CHAIN' 'MAIL' &BCL1;
MAIL_I = CONDINI.DOM_I.'MAIL' ;
$MAIL_I = 'MODE' MAIL_I 'EULER' ;
MMAIL_I = 'DOMA' $MAIL_I 'QUAF' ;
'ELIM' (MMAILTOT 'ET' MMAIL_I) EPSI ;
GEO.MODE_I = $MAIL_I ;
GEO.MQUA_I = MMAIL_I ;
GEO.M_I = 'DOMA' $MAIL_I 'MAILLAGE';
'SI' (&BCL1 'EGA' 1);
LOGX = 'EXISTE' CONDINI.DOM_I 'XH2' ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XO2') ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XN2') ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2O');
LOGY = 'EXISTE' CONDINI.DOM_I 'YH2' ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YO2') ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YN2') ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2O');
'SINON' ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2') ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XO2') ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XN2') ;
LOGX = LOGX 'ET' ('EXISTE' CONDINI.DOM_I 'XH2O');
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2') ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YO2') ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YN2') ;
LOGY = LOGY 'ET' ('EXISTE' CONDINI.DOM_I 'YH2O');
'FINSI';
'FIN' BCL1;
GEO.'MODTOT' = $MAILTOT ;
GEO.'MQUATOT' = MMAILTOT ;
GEO.'MAILTOT' = 'DOMA' $MAILTOT 'MAILLAGE';
GEO.'EPSI' = EPSI ;
*** Verification du bon remplissage de la table de donnees
'SI' (LOGX 'ET' ('NON' LOGY));
'MESS' 'C est bien ! Vous n avez defini que des fractions volumiques';
'SINON';
'SI' (LOGY 'ET' ('NON' LOGX));
'MESS' 'C est bien ! Vous n avez defini que des fractions massiques';
'SINON';
'MESS' 'Vous ne devez definir uniquement des fractions massiques';
'MESS' 'ou uniquement des fractions volumiques !!!';
'MESS' 'En aucun cas, un melange des deux !!!';
ERROR = ERROR + 1 ;
'FINSI';
'FINSI';
'SI' LOGX ;
FRACY.'MH2' = 2. '*' 1.00797E-3 ;
FRACY.'MO2' = 2. '*' 15.9994E-3 ;
FRACY.'MH2O' = (FRACY.'MH2') '+' (0.5 '*' (FRACY.'MO2')) ;
FRACY.'MN2' = 2. '*' 14.0067E-3 ;
'FINSI';
'REPETER' BCL1 CONDINI . 'NBDOM';
DOM_I = 'CHAIN' 'DOM' &BCL1;
FRACY.DOM_I = 'TABLE';
'SI' LOGX;
'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'T'));
'MESS' 'Comme vous avez defini les fractions volumiques,';
'MESS' 'vous devez definir la temperature !';
ERROR = ERROR + 1 ;
'FINSI';
*** CALCUL DES MASSES TOTALES!!!
FRACY.DOM_I.'MTOT' = ((CONDINI.DOM_I.'XH2') '*' (FRACY.'MH2')) '+'
((CONDINI.DOM_I.'XO2') '*' (FRACY.'MO2')) '+'
((CONDINI.DOM_I.'XH2O') '*' (FRACY.'MH2O')) '+'
((CONDINI.DOM_I.'XN2') '*' (FRACY.'MN2')) ;
*** FIN CALCUL DES MASSES TOTALES
'FINSI';
*****************
'SI' LOGY;
'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'RHO'));
'MESS' 'Comme vous avez defini les fractions massiques,';
'MESS' 'vous devez definir la densite du melange !';
ERROR = ERROR + 1 ;
'FINSI';
'FINSI';
*****************
'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'P'));
'MESS' 'Il manque l initialisation de la pression';
'MESS' 'dans le domaine' DOM_I;
ERROR = ERROR + 1 ;
'FINSI';
*****************
'SI' ('NON' ('EXISTE' CONDINI.DOM_I 'V'));
'MESS' 'Il manque l initialisation de la vitesse';
'MESS' 'dans le domaine' DOM_I;
ERROR = ERROR + 1 ;
'FINSI';
'FIN' BCL1;
*****************
'SI' (ERROR 'NEG' 0);
'MESS' 'Il y a des problèmes d initialisation !';
'ERRE' 5;
'FINSI';
****** TRANSFORMATION DES FRACTIONS VOLUMIQUES EN MASSIQUES
'SI' LOGX ;
'REPETER' BCL1 CONDINI . 'NBDOM';
DOM_I = 'CHAIN' 'DOM' &BCL1;
FRACY.DOM_I.'YH2' = ((CONDINI.DOM_I.'XH2') '*' (FRACY.'MH2')) '/'
(FRACY.DOM_I.'MTOT') ;
FRACY.DOM_I.'YO2' = ((CONDINI.DOM_I.'XO2') '*' FRACY.'MO2') '/'
(FRACY.DOM_I.'MTOT') ;
FRACY.DOM_I.'YH2O'= ((CONDINI.DOM_I.'XH2O') '*' FRACY.'MH2O') '/'
(FRACY.DOM_I.'MTOT') ;
FRACY.DOM_I.'YN2' = ((CONDINI.DOM_I.'XN2') '*' FRACY.'MN2') '/'
(FRACY.DOM_I.'MTOT') ;
FRACY.DOM_I.'RTOT'= ((FRACY.DOM_I.'YH2') '*' (PGAZ.'H2'.'R')) '+'
((FRACY.DOM_I.'YO2') '*' (PGAZ.'O2'.'R')) '+'
((FRACY.DOM_I.'YH2O') '*' (PGAZ.'H2O'.'R')) '+'
((FRACY.DOM_I.'YN2') '*' (PGAZ.'N2'.'R')) ;
FRACY.DOM_I.'RHO' = (CONDINI.DOM_I.'P') '/'
((FRACY.DOM_I.'RTOT') '*' (CONDINI.DOM_I.'T'));
'FIN' BCL1;
'FINSI';
**** INITIALISATION DES CHAMPOINTS POUR LE CALCUL
'REPETER' BCL1 CONDINI . 'NBDOM';
MODE_I = 'CHAIN' 'MOD' &BCL1;
DOM_I = 'CHAIN' 'DOM' &BCL1;
RHO_I = 'CHAIN' 'RHO' &BCL1;
P_I = 'CHAIN' 'P' &BCL1;
V_I = 'CHAIN' 'V' &BCL1;
YH2_I = 'CHAIN' 'YH2' &BCL1;
YO2_I = 'CHAIN' 'YO2' &BCL1;
YH2O_I = 'CHAIN' 'YH2O' &BCL1;
INI . P_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
CONDINI.DOM_I.'P';
INI . V_I = 'KCHT' GEO.MODE_I 'VECT' 'CENTRE'
CONDINI.DOM_I.'V';
'SI' LOGY;
INI . RHO_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
CONDINI.DOM_I.'RHO';
INI . YH2_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'H2 ' CONDINI.DOM_I.'YH2';
INI . YO2_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'O2 ' CONDINI.DOM_I.'YO2';
INI . YH2O_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'H2O ' CONDINI.DOM_I.'YH2O';
'SINON';
'SI' LOGX;
INI . RHO_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
FRACY.DOM_I.'RHO';
INI . YH2_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'H2 ' FRACY.DOM_I.'YH2';
INI . YO2_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'O2 ' FRACY.DOM_I.'YO2';
INI . YH2O_I = 'KCHT' GEO.MODE_I 'SCAL' 'CENTRE'
'COMP' 'H2O ' FRACY.DOM_I.'YH2O';
'FINSI';
'FINSI';
'SI' (&BCL1 'EGA' 1);
RN = 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' INI . RHO_I ;
PN = 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' INI . P_I ;
VN = 'KCHT' GEO.'MODTOT' 'VECT' 'CENTRE' INI . V_I ;
YH2= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'H2 ' INI . YH2_I ;
YO2= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'O2 ' INI . YO2_I ;
YH2O= 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'H2O ' INI . YH2O_I ;
'SINON';
RN = RN 'ET' ('KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' INI . RHO_I);
PN = PN 'ET' ('KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' INI . P_I );
VN = VN 'ET' ('KCHT' GEO.'MODTOT' 'VECT' 'CENTRE' INI . V_I );
YH2=YH2 'ET' ('KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'H2 ' INI . YH2_I );
YO2=YO2 'ET' ('KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'O2 ' INI . YO2_I );
YH2O=YH2O 'ET' ('KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE'
'COMP' 'H2O ' INI . YH2O_I );
'FINSI';
'FIN' BCL1;
*** On supprime les tables qui ne nous sont plus necessaires
'OUBLIER' INI ;
'OUBLIER' FRACY ;
***
INCO . 'RN' = 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' RN ;
INCO . 'PN' = 'KCHT' GEO.'MODTOT' 'SCAL' 'CENTRE' PN ;
INCO . 'VN' = 'KCHT' GEO.'MODTOT' 'VECT' 'CENTRE' VN ;
YN = YH2 'ET' YO2 'ET' YH2O ;
INCO . 'YN' = YN ;
'SI' CONDINI.'GRAPH';
'SI' (('VALE' DIME) 'EGA' 2);
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' RN 'CHAM')
('CONT' GEO.'MAILTOT') 'CACH' 'TITR' 'RO' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YH2 'CHAM')
('CONT' GEO.'MAILTOT') 'CACH' 'TITR' 'YH2' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YO2 'CHAM')
('CONT' GEO.'MAILTOT') 'CACH' 'TITR' 'YO2' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YH2O 'CHAM')
('CONT' GEO.'MAILTOT') 'CACH' 'TITR' 'YH2O';
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' PN 'CHAM')
('CONT' GEO.'MAILTOT') 'CACH' 'TITR' 'P' ;
'SINON';
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' RN 'CHAM')
'CACH' 'TITR' 'RO' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YH2 'CHAM')
'CACH' 'TITR' 'YH2' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YO2 'CHAM')
'CACH' 'TITR' 'YO2' ;
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' YH2O 'CHAM')
'CACH' 'TITR' 'YH2O';
'TRAC' GEO.'MODTOT' (KCHA GEO.'MODTOT' PN 'CHAM')
'CACH' 'TITR' 'P' ;
'FINSI';
'FINSI';
GN REN RYN = CONS RN VN PN YN GEO.'MODTOT' PGAZ ;
***************************
**** La table PASDETPS ****
***************************
PASDETPS = 'TABLE' 'PASDETPS' ;
'SI' ('EXISTE' CONDINI 'CFL');
PASDETPS.'CFL' = CONDINI . 'CFL' ;
'SINON';
PASDETPS.'CFL' = 1.0D0 ;
'MESS' 'On prend CFL par défaut = ' CFL ;
'FINSI';
'SI' ('EXISTE' CONDINI 'NITMA');
PASDETPS.'NITMA' = CONDINI . 'NITMA' ;
'SINON';
PASDETPS.'NITMA' = 10000 ;
'MESS' 'On prend ITMA par défaut = ' NITMA ;
'FINSI';
'SI' ('EXISTE' CONDINI 'TINI');
PASDETPS.'TINI' = CONDINI . 'TINI' ;
'SINON';
PASDETPS.'TINI' = 0.0 ;
'MESS' 'Temps initial = ' PASDETPS.'TINI';
'FINSI';
PASDETPS.'LTPS' = 'PROG' PASDETPS.'TINI' ;
PASDETPS.'TPS' = PASDETPS.'TINI' ;
'SI' ('EXISTE' CONDINI 'TFINA');
PASDETPS.'TFINA' = CONDINI . 'TFINA' ;
'SINON';
'MESS' 'Vous devez entrer un temps final';
ERROR = ERROR + 1 ;
'FINSI';
'SI' ('EXISTE' CONDINI 'DTIMP');
PASDETPS.'DTIMP' = CONDINI . 'DTIMP' ;
'SINON';
'MESS' 'Vous devez entrer un pas de temps';
ERROR = ERROR + 1 ;
'FINSI';
PASDETPS.'DT' = 'PROG' ;
*
****
*
'SI' ('EXISTE' CONDINI 'SAFFAC');
SAFFAC = CONDINI . 'SAFFAC' ;
'SINON';
SAFFAC = 0.7 ;
'MESS' 'On prend le facteur de securite pour le pas de temps';
'MESS' ' par défaut = ' SAFFAC ;
'FINSI';
PASDETPS.'SAFFAC' = SAFFAC ;
'SI' ('EXISTE' CONDINI 'METO');
METO = CONDINI . 'METO' ;
'SINON';
METO = 'VLH ';
'MESS' 'Methode VLH par defaut !!!';
'FINSI';
CONDINI.'METO' = METO ;
RV = 'TABLE' ;
RV . 'INCO' = INCO ;
'SI' ('EXISTE' 'TABC');
RV . 'TABC' = TABC ;
'FINSI';
RV . 'CI' = CONDINI ;
RV . 'GEO' = GEO ;
RV . 'PASDETPS'= PASDETPS ;
RV . 'INCO' . 'RN' = 'COPIER' RN ;
RV . 'INCO' . 'GN' = 'COPIER' GN ;
RV . 'INCO' . 'RETN' = 'COPIER' REN ;
RV . 'INCO' . 'RYN' = 'COPIER' RYN ;
*
*** Le gaz
*
RV . 'PGAZ' = PGAZ ;
*
**** Ordre en espace
*
'SI' ('EXISTE' CONDINI 'ORDREESP');
IE = CONDINI . 'ORDREESP' ;
'SINON';
IE = 1 ;
'MESS' 'On prend ORDREESP par défaut = ' IE ;
'FINSI';
RV . 'ORDREESP' = IE ;
*
**** Ordre en temps
*
'SI' ('EXISTE' CONDINI 'ORDRETPS');
IT = CONDINI . 'ORDRETPS' ;
'SINON';
IT = 1 ;
'MESS' 'On prend ORDRETPS par défaut = ' IT ;
'FINSI';
RV . 'ORDRETPS' = IT ;
*
**** Definition du modele de combustion/detonation
*
'SI' ('EXISTE' CONDINI 'MODETO');
MODETO = CONDINI . 'MODETO' ;
'SINON';
MODETO = 'ARRHENIU' ;
'MESS' 'On prend comme modele de combustion/deto ';
'MESS' ' par défaut = ' MODETO ;
'FINSI';
RV . 'MODETO' = 'TABLE' 'DETO' ;
RV . 'MODETO' .'NOM' = MODETO ;
**** Si le modele de combustion est CREBCOM :
'SI' ('EGA' MODETO 'CREBCOM');
* numero du domaine d'ignition
'SI' ('EXISTE' CONDINI 'IGNIREG');
RV . 'MODETO' .'IGNIREG' = CONDINI . 'IGNIREG' ;
'SINON';
RV . 'MODETO' .'IGNIREG' = 1 ;
'MESS' 'On prend le 1er domaine comme domaine d ignition';
'FINSI';
* critere Crebcom (chimie)
'SI' ('EXISTE' CONDINI 'EPS_CC');
RV . 'MODETO' .'EPS_CC' = CONDINI . 'EPS_CC' ;
'SINON';
RV . 'MODETO' .'EPS_CC' = 0.7 ;
'MESS' 'On prend par défaut EPS_CC = ' RV . 'MODETO' .'EPS_CC';
'FINSI';
* critere Crebcom (vitesse)
'SI' ('EXISTE' CONDINI 'K0');
RV . 'MODETO' .'K0' = CONDINI . 'K0' ;
'SINON';
RV . 'MODETO' .'K0' = 500.0 ;
'MESS' 'On prend par défaut EPS_CC = ' RV . 'MODETO' .'K0';
'FINSI';
'FINSI';
**** Si le modele de combustion est ARRHENIUS :
'SI' ('EGA' MODETO 'ARRHENIU');
* numero du domaine d'ignition
'SI' ('EXISTE' CONDINI 'IGNIREG');
RV . 'MODETO' .'IGNIREG' = CONDINI . 'IGNIREG' ;
'SINON';
RV . 'MODETO' .'IGNIREG' = 1 ;
'MESS' 'On prend le 1er domaine comme domaine d ignition';
'FINSI';
'FINSI';
*
****
*
'SI' (ERROR 'NEG' 0);
'MESS' 'Il y a des problèmes d initialisation !';
'ERRE' 5;
'FINSI';
'FINPROC' RV;
*$$$$ VALP
*****************************************************
*****************************************************
** PROCEDURE VALP POUR CALCULER LES VALEURS DE **
** PRESSION EN DES POINTS DONNES ET A DES TEMPS **
** DONNES **
*****************************************************
*****************************************************
'DEBPROC' VALP RV*'TABLE' ;
TABC = RV . 'TABC' ;
MODTOT = RV . 'GEO' . 'MODTOT' ;
CENTOT = 'DOMA' MODTOT 'CENTRE';
PTS = TABC . 'LPOINTS' ;
M1 = 'NBEL' PTS ;
LTPS = TABC . 'LTPS' ;
NTPS = DIME LTPS ;
L1 = TABC . 'COMP' ;
N1 = 'DIME' L1 ;
* On etudie uniquement la pression
* peut evoluer par la suite !!!
LOG1 = 'EXISTE' TABC 'NOMC' ;
* Initialisation des tables
'SI' (NTPS EGA 0 );
'REPETER' BOU1 M1;
PT1 = PTS 'POIN' &BOU1;
'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
'SINON' ;
CAPI = 'CHAINE' 'CAPT' &BOU1;
'FINSI';
PT2 = CENTOT 'POIN' 'PROCHE' PT1;
'SI' (&BOU1 'EGA' 1);
TABC.'RPOINTS' = PT2 ;
'SINON';
TABC.'RPOINTS' = RV.'TABC'.'RPOINTS' 'ET' PT2;
'FINSI';
TABC . CAPI = 'TABLE';
'REPETER' BOU2 N1;
MOI = EXTR L1 &BOU2 ;
TABC . CAPI . MOI = 'PROG';
'FIN' BOU2;
'FIN' BOU1;
'SI' RV.'CI'.'GRAPH' ;
'SI' (('VALE' DIME) 'EGA' 2) ;
CONT1 = 'CONT' RV.'GEO'.'MAILTOT' ;
TRACE (('COUL' TABC.'RPOINTS' 'ROUGE') 'ET' CONT1 ) ;
'SINON';
ENV1 = 'ENVE' RV.'GEO'.'MAILTOT' ;
TRACE (('COUL' TABC.'RPOINTS' 'ROUGE') 'ET' ENV1 ) ;
'FINSI';
'FINSI';
'FINSI';
'REPETER' BOU1 M1;
PT1 = TABC.'RPOINTS' 'POIN' &BOU1;
'SI' (LOG1 'ET' (('DIME' (TABC.'NOMC')) 'EGA' M1 )) ;
CAPI = 'EXTR' (TABC.'NOMC') &BOU1;
'SINON' ;
CAPI = 'CHAINE' 'CAPT' &BOU1;
'FINSI';
'REPETER' BOU2 N1;
MOI = EXTR L1 &BOU2 ;
*** captage de la pression
'SI' ('EGA' MOI 'PN');
X1 = 'EXTR' (RV . 'INCO' . 'PN') 'SCAL' PT1 ;
TABC . CAPI . MOI = TABC . CAPI . MOI 'ET' ('PROG' X1);
'FINSI';
*** captage de temperature
'SI' ('EGA' MOI 'TN');
X1 = 'EXTR' (RV . 'INCO' . 'TN') 'SCAL' PT1 ;
TABC . CAPI . MOI = TABC . CAPI . MOI 'ET' ('PROG' X1);
'FINSI';
'FIN' BOU2;
'FIN' BOU1;
TABC . 'LTPS' = TABC . 'LTPS' 'ET' ('PROG' RV . 'PASDETPS' . 'TPS');
'FINPROC' TABC;
************************************************************************
***************** FIN DEFINTION DES PROCEDURES *************************
************************************************************************
***************** REMPLISSAGE DE LA TABLE PGAZ *************************
*
***** 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 ;
*'OPTION' 'SAUVER' 'FORMAT' 'pgaz.sauv' ;
*'SAUVER' 'FORMAT' PGAZ ;
*'FIN';
*** ------------------------------------------------- ***
*** Fichier correspondant a la definition du maillage ***
*** ------------------------------------------------- ***
************
* MAILLAGE *
************
raf = 1 ;
'OPTI' 'DIME' 2 'ELEM' 'QUA4' ;
*
*
* Maillage du canal RUT sans obstacles
*
*
*NX1 = 10 ; NX2 = 20 ; NX3 = 106 ; NX4 = 180 ;
*DINI1 = 0.1 ; DFIN1 = 1. ;
*NY1 = 15 ; NY2 = 25 ; NY3 = 10 ; NY4 = 13 ;
X0 = 31.6 ;
X1 = X0 ;
Y0 = 0.0 ;
Y1 = Y0 ;
IGNI1 = 1. ;
LC1 = 3. ;
LK1 = 10.6 ;
LC2 = 18. ;
LC3 = 54. ;
HK1 = 4.;
HC = 2.3 ;
*YJH = Y1 + 1. - ((HC - 1.)/(2. * NY4 - 1.));repositionnement
*YJB = 1.5 - HK1 + ((1.5 - Y1) / ( 2. * NY1 - 1.)) ;repositionnement
YJH = Y1 + 1. ;
YJB = 1.5 - HK1 ;
*---------------------------------------> Points du contour
P1 = X1 Y1 ;
P2 = (X1 + LC1) Y1 ;
P3 = (X1 + LC1) (Y1 - HK1) ;
P4 = (X1 + LC1 + LK1) (Y1 - HK1) ;
P5 = (X1 + LC1 + LK1) Y1 ;
PF1 = (X1 + LC1 + LK1 + LC2) Y1 ;
PF2 = (X1 + LC1 + LK1 + LC2) (Y1 + HC) ;
P6 = (X1 + LC1 + LK1 + LC2 + LC3) Y1 ;
P7 = (X1 + LC1 + LK1 + LC2 + LC3) (Y1 + HC) ;
P8 = (X1 + LC1 + LK1) (Y1 + HC) ;
P9 = (X1 + LC1) (Y1 + HC) ;
P10 = X1 (Y1 + HC) ;
*---------------------------------------> Points pour les jauges
J1 = X1 YJH ;
J2 = (X1 + IGNI1) YJH ;
J2B = (X1 + IGNI1) Y1 ;
J2H = (X1 + IGNI1) (Y1 + HC) ;
J3 = (X1 + LC1) YJH ;
J4 = (X1 + LC1 + LK1) YJH ;
J5 = (X1 + LC1 + LK1 + LC2) YJH ;
J6 = (X1 + LC1 + LK1 + LC2 + LC3) YJH ;
J7 = (X1 + LC1) YJB ;
J8 = (X1 + LC1 + LK1) YJB ;
*
*- Lignes pour les sous-domaines
*
* MAILLAGE FIN
*NX1 = 10 ; NX2 = 20 ; NX3 = 106 ; NX4 = 180 ;
*DINI1 = 0.1 ; DFIN1 = 1. ;
*NY1 = 15 ; NY2 = 25 ; NY3 = 10 ; NY4 = 13 ;
* MAILLAGE TRES FIN
*NX1 = 15 ; NX2 = 30 ; NX3 = 150 ; NX4 = 270 ;
*DINI1 = 0.1 ; DFIN1 = 1. ;
*NY1 = 22 ; NY2 = 37 ; NY3 = 15 ; NY4 = 19 ;
* MAILLAGE TRES GROSSIER
'MESS' 'Maillage tres grossier !' ;
NX1 = 2 ; NX2 = 4 ; NX3 = 25 ; NX4 = 36 ;
DINI1 = 0.5 ; DFIN1 = 5. ;
NY1 = 3 ; NY2 = 5 ; NY3 = 2 ; NY4 = 3 ;
* MAILLAGE GROSSIER
*NX1 = 4 ; NX2 = 8 ; NX3 = 50 ; NX4 = 72 ;
*DINI1 = 0.25 ; DFIN1 = 2.5 ;
*NY1 = 6 ; NY2 = 10 ; NY3 = 4 ; NY4 = 6 ;
*---------------------------------------> Contour externe
P1J2B = P1 'DROI' NX1 J2B ;
J2BP2 = J2B 'DROI' NX2 P2 ;
P2J7 = P2 'DROI' NY2 J7 ;
J7P3 = J7 'DROI' NY1 P3 ;
P3P4 = P3 'DROI' NX3 P4 ;
P4J8 = P4 'DROI' NY1 J8 ;
J8P5 = J8 'DROI' NY2 P5 ;
P5PF1 = P5 'DROI' NX4 PF1 ;
PF1P6 = PF1 'DROI' P6 'DINI' DINI1 'DFIN' DFIN1 ;
P6J6 = P6 'DROI' NY3 J6 ;
J6P7 = J6 'DROI' NY4 P7 ;
P7PF2 = P7 'DROI' PF2 'DINI' DFIN1 'DFIN' DINI1 ;
PF2P8 = PF2 'DROI' NX4 P8 ;
P8P9 = P8 'DROI' NX3 P9 ;
P9J2H = P9 'DROI' NX2 J2H ;
J2HP10 = J2H 'DROI' NX1 P10 ;
P10J1 = P10 'DROI' NY4 J1 ;
J1P1 = J1 'DROI' NY3 P1 ;
*----------------------------------------> Jauges
J1J2 = J1 'DROI' NX1 J2 ;
J2J3 = J2 'DROI' NX2 J3 ;
J3J4 = J3 'DROI' NX3 J4 ;
J4J5 = J4 'DROI' NX4 J5 ;
J5J6 = J5 'DROI' J6 'DINI' DINI1 'DFIN' DFIN1 ;
J7J8 = J7 'DROI' NX3 J8 ;
*----------------------------------------> Autres lignes internes
J2BJ2 = J2B 'DROI' NY3 J2 ;
J2J2H = J2 'DROI' NY4 J2H ;
P2J3 = P2 'DROI' NY3 J3 ;
J3P9 = J3 'DROI' NY4 P9 ;
P5J4 = P5 'DROI' NY3 J4 ;
J4P8 = J4 'DROI' NY4 P8 ;
PF1J5 = PF1 'DROI' NY3 J5 ;
J5PF2 = J5 'DROI' NY4 PF2 ;
P2P5 = P2 'DROI' NX3 P5 ;
J7J8 = J7 'DROI' NX3 J8 ;
*
*- Pavage de chacune de zones identifiées
*
*---------------------------------------> Premier canal
DC1 = 'DALL' P1J2B J2BJ2 ('INVE' J1J2) J1P1 ;
DC2 = 'DALL' J1J2 J2J2H J2HP10 P10J1 ;
DOM1 = DC1 ET DC2 ;
DC3 = 'DALL' J2BP2 P2J3 ('INVE' J2J3) ('INVE' J2BJ2) ;
DC4 = 'DALL' J2J3 J3P9 P9J2H ('INVE' J2J2H) ;
DOM2 = DC3 ET DC4 ;
DOMC1 = DOM1 ET DOM2 ;
*---------------------------------------> Deuxieme canal
DC5 = 'DALL' P5PF1 PF1J5 ('INVE' J4J5) ('INVE' P5J4) ;
DC6 = 'DALL' J4J5 J5PF2 PF2P8 ('INVE' J4P8) ;
DOM3 = DC5 ET DC6 ;
DC7 = 'DALL' PF1P6 P6J6 ('INVE' J5J6) ('INVE' PF1J5) ;
DC8 = 'DALL' J5J6 J6P7 P7PF2 ('INVE' J5PF2) ;
DOM4 = DC7 ET DC8 ;
DOMC2 = DOM3 ET DOM4 ;
*---------------------------------------> Canyon
DK1 = 'DALL' P3P4 P4J8 ('INVE' J7J8) J7P3 ;
DK2 = 'DALL' J7J8 J8P5 ('INVE' P2P5) P2J7 ;
DK3 = 'DALL' P2P5 P5J4 ('INVE' J3J4) ('INVE' P2J3) ;
DK4 = 'DALL' J3J4 J4P8 P8P9 ('INVE' J3P9) ;
DOMK1 = DK1 ET DK2 ET DK3 ET DK4 ;
*---------------------------------------> Maillage global
DOMTOT = DOMC1 ET DOMC2 ET DOMK1 ;
DDOM1 = DOM1 ;
DDOM2 = DOM2 ET DOMC2 ET DOMK1 ;
*
*** PCEN1 Pour determiner PRES(PCEN,t)
*
XINIT0 = 0.33 ;
*repositionnement
*XINIT0 = 1.41 ;
YINIT0 = 0.79 ;
PCEL2 = (XINIT0 + 34.78) YJB ;
PCEL4 = (XINIT0 + 41.29) YJB ;
PCEL5 = (XINIT0 + 42.97) YJB ;
PCEL6 = (XINIT0 + 44.84) YJB ;
PCEL7 = (XINIT0 + 34.54) YJH ;
PCEL8 = (XINIT0 + 36.79) YJH ;
PCEL9 = (XINIT0 + 39.04) YJH ;
PCEL10= (XINIT0 + 41.74) YJH ;
PCEL11= (XINIT0 + 45.04) YJH ;
PCEL12= (XINIT0 + 50.02) YJH ;
PCEL13= (XINIT0 + 55.04) YJH ;
PCEL14= (XINIT0 + 64.00) YJH ;
*
*PCEL2 = 35.11 YJB ;
*PCEL4 = 41.62 YJB ;
*PCEL5 = 43.30 YJB ;
*PCEL6 = 45.17 YJB ;
*PCEL7 = 34.87 YJH ;
*PCEL8 = 37.12 YJH ;
*PCEL9 = 39.37 YJH ;
*PCEL10= 42.07 YJH ;
*PCEL11= 45.37 YJH ;
*PCEL12= 50.35 YJH ;
*PCEL13= 55.37 YJH ;
*PCEL14= 64.33 YJH ;
LPOINTS = PCEL2 'ET' PCEL4 'ET' PCEL5 'ET' PCEL6
'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 'ET' PCEL10 'ET' PCEL11
'ET' PCEL12 'ET' PCEL13 'ET' PCEL14 ;
*'OPTI' 'SAUV' 'FORMAT' 'maillage_rut_fin.save';
*'OPTI' 'SAUV' 'FORMAT' 'maillage_rut_gros.save';
*'SAUV' 'FORMAT' ;
*'FIN';
******* FICHIER MAITRE *************************************
*** RECUPERATION DES OBJETS MODELES + POINTS DES CAPTEURS ***
*'OPTI' 'REST' 'FORMAT' 'maillage_rut_gros.save';
*'OPTI' 'REST' 'FORMAT' 'maillage_rut_fin.save';
*'REST' 'FORMAT' ;
*
TABC = 'TABLE' 'CAPTEUR' ;
TABC . 'LPOINTS' = PCEL2 'ET' PCEL4 'ET' PCEL5 'ET' PCEL6
'ET' PCEL7 'ET' PCEL8 'ET' PCEL9 'ET' PCEL10 'ET' PCEL11
'ET' PCEL12 'ET' PCEL13 'ET' PCEL14 ;
TABC . 'LTPS' = 'PROG';
TABC . 'NOMC' = 'MOTS' 'C2' 'C4' 'C5' 'C6' 'C7'
'C8' 'C9' 'C10' 'C11' 'C12'
'C13' 'C14' ;
TABC . 'COMP' = 'MOTS' 'PN' 'TN';
TABC . 'FREQ' = 1 ;
*** DEFINITION DES GAZ EN PRESENCE + LEURS PROPRIETES ********
***********************
**** LA TABLE PGAZ ****
***********************
*'OPTI' 'REST' 'FORMAT' 'pgaz.sauv';
*'REST' 'FORMAT' ;
*
*** Fin PGAZ
*
*** DEFINITION DES CONDITIONS INITIALES *************************
CONDINI = 'TABLE' ;
CONDINI . 'NBDOM' = 2 ;
* 1.1) DOM1 = gauche
* 1.2) DOM2 = droit
CONDINI . 'DOM1' = 'TABLE' ;
CONDINI . 'DOM2' = 'TABLE' ;
* ZONE GAUCHE
tg = 3000.0 ;
pg =20.00D5 ;
xh2sec = 0.250 ;
xh2o = 0.250 ;
xh2 = xh2sec * (1.D0 - xh2o) ;
xo2 = 0.21 * (1.D0 - xh2 - xh2o) ;
xn2 = 1.0D0 - (xh2 + xo2 + xh2O) ;
uxg = 0.0 ;
uyg = 0.0 ;
CONDINI . 'DOM1' .'MAIL' = DDOM1;
CONDINI . 'DOM1' .'T' = tg ;
CONDINI . 'DOM1' .'P' = pg ;
CONDINI . 'DOM1' .'V' = (uxg uyg) ;
CONDINI . 'DOM1' .'XH2' = xh2 ;
CONDINI . 'DOM1' .'XO2' = xo2 ;
CONDINI . 'DOM1' .'XN2' = xn2 ;
CONDINI . 'DOM1' .'XH2O' = xh2o ;
* ZONE DROITE
td = 363. ;
und = 0.0D0 ;
utd = 0.0D0 ;
pd = 1.0D5 ;
CONDINI . 'DOM2' .'MAIL' = DDOM2;
CONDINI . 'DOM2' .'T' = td ;
CONDINI . 'DOM2' .'P' = pd ;
CONDINI . 'DOM2' .'V' = (und utd) ;
CONDINI . 'DOM2' .'XH2' = xh2 ;
CONDINI . 'DOM2' .'XO2' = xo2 ;
CONDINI . 'DOM2' .'XN2' = xn2 ;
CONDINI . 'DOM2' .'XH2O' = xh2o ;
*** Methodes possibles :
*
* 'VLH '
CONDINI . 'METO' = 'VLH ' ;
CONDINI . 'MAILTOT' = DOMTOT ;
CONDINI . 'TINI' = 0.0 ;
CONDINI . 'TFINA' = 8.D-3 ;
CONDINI . 'ORDREESP' = 2 ;
* Maillage fin
**CONDINI . 'DTIMP' = 1.D-5 ;
* Maillage tres grossier
CONDINI . 'DTIMP' = 3.D-5 ;
*CONDINI . 'GRAPH' = VRAI ;
CONDINI . 'GRAPH' = FAUX ;
'OPTI' 'TRACE' 'PSC' ;
*** FIN remplissage de la table definissant les conditions initiales
*** preparation des equations et initialisation des champoints
*** recuperation de la table de calculs
RV = PREPADET CONDINI PGAZ TABC ;
MDOMT = RV.'GEO'.'MQUATOT' ;
DOMT1 = RV.'GEO'.'MAILTOT';
X1 = 'COORD' 1 DOMT1;
X2 = 'CHAN' 'CHAM' X1 DOMT1;
X3 = 'ELEM' X2 'INFERIEUR' 48. ;
X4 = 'CHAN' X3 'QUAF' ;
'ELIM' (X4 ET MDOMT) RV.'GEO'.'EPSI';
$DOMT1 = 'MODE' X4 'NAVIER_STOKES' 'LINE' ;
TDOMT1 = 'DOMA' $DOMT1 'TABLE' ;
TDOMT1 . 'PRECONDI' = 1 ;
* Table perso pour le trace de courbes au cours du temps
*RV.'PRCPERSO' = 'TABLE' 'PERSO' ;
*RV.'PRCPERSO'.'NOMPROC' = 'MOT' PERSO ;
*RV.'PRCPERSO'.'ARG1' = $DOMT1 ;
*RV.'PRCPERSO'.'ARG2' = 1 ;
* Maillage fin
**RV.'PRCPERSO'.'ARG3' = 15 ;
* Maillage grossier
*RV.'PRCPERSO'.'ARG3' = 5 ;
*** execution du calcul
EXEX RV ;
*** relancement d un calcul
*RV . 'PASDETPS' . 'TFINA' = 6.D-3 ;
* EXEX RV ;
* Maillage fin
**'OPTION' 'SAUV' 'main_rut_fin.sauv' ;
* Maillage grossier
*'OPTION' 'SAUV' 'main_rut_gros.sauv' ;
*'SAUV' ;
*************** TESTS de non-regression ********************************
TPS = RV.'TABC'.'LTPS';
CA = 'TABLE' ;
CA . 1 = RV.'TABC'.'C7'.'PN';
CA . 2 = RV.'TABC'.'C8'.'PN';
CA . 3 = RV.'TABC'.'C9'.'PN';
CA . 4 = RV.'TABC'.'C10'.'PN';
N1 = 'DIME' TPS;
TREF = 'PROG' 1.86E-3 3.00E-3 4.53E-3 6.12E-3 ;
PREF = 'PROG' 1.12503E+06 1.02158E+06 9.41526E+05 8.81071E+05 ;
PCAL = 'PROG' ;
TOL1 = 1.E-6 ;
NN1 = 1 ;
'REPETER' BLO1 N1 ;
T1 = EXTR TPS &BLO1 ;
T2 = EXTR TREF NN1 ;
'SI' (T1 'EGA' T2 TOL1);
P1 = 'EXTR' (CA . NN1) &BLO1 ;
PCAL= PCAL 'ET' ('PROG' P1) ;
NN1 = NN1 + 1 ;
'MESS' 'Tps/pression : ' T1 P1 ;
'FINSI';
'SI' (NN1 '>' ('DIME' PREF));
'QUITTER' BLO1 ;
'FINSI';
'FIN' BLO1 ;
'MESS' 'Pression de reference' ;
'LIST' PREF;
'MESS' 'Pression calculee' ;
'LIST' PCAL;
DIFF1 = 'MAXI' (PREF - PCAL) 'ABS' ;
* AB: Je definis PREFS pour calculer une erreur relative.
PREFS = 'MINIMUM' PREF ;
EPS_P = 1.E-5 ;
ERROR = 0 ;
**** TEST sur la pression
'MESS' 'Max difference : ' (DIFF1 '/' PREFS) ;
'SI' ((DIFF1 '/' PREFS) '>' EPS_P) ;
'MESS' 'Probleme sur la pression ' ;
ERROR = ERROR + 1 ;
'FINSI';
**** Calcul de la vitesse de flamme
X1 Y1 = COORD (RV.'TABC'.'RPOINTS');
P0 = EXTR X1 MAILLAGE;
P7 = POINT P0 5;
P8 = POINT P0 6;
P9 = POINT P0 7;
P10 = POINT P0 8;
X7 = EXTR X1 'SCAL' P7 ;
X8 = EXTR X1 'SCAL' P8 ;
X9 = EXTR X1 'SCAL' P9 ;
X10 = EXTR X1 'SCAL' P10 ;
DX8 = X8 - X7 ;
DX9 = X9 - X8 ;
DX10 = X10 - X9 ;
T2 = 'ENLE' TREF 1 ;
T1 = 'ENLE' TREF ('DIME' TREF) ;
DT1 = T2 - T1 ;
DX1 = 'PROG' DX8 DX9 DX10 ;
US1 = DX1 / DT1 ;
'MESS' 'Vitesse de flamme : '; 'LIST' US1 ;
'SI' (ERROR '>' 0) ;
'MESS' 'Il y a ' ERROR ' probleme(s) ' ;
'ERRE' 5 ;
'FINSI';
'FIN';
© Cast3M 2003 - All rights reserved.
Disclaimer