*****************************************************
************************************************************************
* Section : Chimie Combustion
************************************************************************
* fichier :  flamcat.dgibi                          *
** modifie le 15/06/2014 passage EQPR -> EQEX       *
*****************************************************
*****************************************
* COMBUSTION CATALYTIQUE DE L'HYDROGENE *
* on teste l'operateur FLAM             *
* H.PAILLERE/TTMF 97                    *
*****************************************
   DISCR  = 'MACRO'    ;
   KPRESS = 'CENTRE'   ;
   BETAP  = 1.         ;

 DEBP FILTREKE ;
 ARGU RX*TABLE ;
 rv=rx.'EQEX' ;
 rvp=rv.'PROJ' ;

* Filtre sur K et Epsilon
* K > 0  Epsilon tel que l'echelle de longueur reste inférieure
* à L0 diamètre enceinte

 rv=rx.'EQEX' ;
 iarg=rx.'IARG' ;
 si( non ( ega iarg 3)) ;
 mess 'Procedure FILTREKE : nombre d arguments incorrect ' iarg ;
 quitter FILTREKE  ;
 finsi ;
 si ( ega ('TYPE' rx.'ARG1') 'MOT     ') ;
 U1=rv.'INCO'.(rx.'ARG1') ;
 sinon ;
 si ( ega ('TYPE' (rx.'ARG1')) 'CHPOINT') ;
 U1=rx.'ARG1' ;
 sinon ;
 mess 'Procedure FILTREKE : type argument 1 invalide ' ;
 quitter FILTREKE ;
 finsi ;
 finsi ;

 si ( ega ('TYPE' rx.'ARG2') 'MOT     ') ;
 L0=rv.'INCO'.(rx.'ARG2') ;
 sinon ;
 si ( ega ('TYPE' (rx.'ARG2')) 'FLOTTANT') ;
 L0=rx.'ARG2' ;
 sinon ;
 mess 'Procedure FILTREKE : type argument 2 invalide ' ;
 quitter FILTREKE ;
 finsi ;
 finsi ;

 si ( ega ('TYPE' rx.'ARG3') 'MOT     ') ;
 NU=rv.'INCO'.(rx.'ARG3') ;
 sinon ;
 si ( ega ('TYPE' (rx.'ARG3')) 'FLOTTANT') ;
 NU=rx.'ARG3' ;
 sinon ;
 mess 'Procedure FILTREKE : type argument 3 invalide ' ;
 quitter FILTREKE ;
 finsi ;
 finsi ;

 nic=dime (rx.'LISTINCO') ;
 si( non ( ega nic 2)) ;
 mess 'Procedure FILTREKE : nombre d inconnues incorrect ' nic ;
 quitter FILTREKE ;
 finsi ;
 nomi1=extr 1 (rx.'LISTINCO');
 nomi2=extr 2 (rx.'LISTINCO');
 nom1= mot (text (chai nomi1)) ;
 nom2= mot (text (chai nomi2)) ;

 en=rv.'INCO'.nom2 ;
 kn=rv.'INCO'.nom1 ;
 Rec=100.;
 k0 = 1.e-10 ;

 lcu1=extr u1 'COMP'       ;
 mdu=u1 lcu1 'PSCA' u1 lcu1;
 mdu=mdu ** 0.5 ;
 Re=kops (kops (kops mdu '*' L0) '/' nu) '+' (Rec / 10.) ;
 a= exp (kops Rec '/' Re ) ;
 kn=kops kn '|<' k0 ;
 E0= kops (kops kn '**' 1.5) '*' (a / L0) ;
 en=kops en '|<' E0 ;

 rv.'INCO'.'EN'=en ;
 rv.'INCO'.'KN'=kn ;

 as2 ama1  = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1;


*********************************************************************
* PROCEDURE RECOMBINEUR
*********************************************************************

DEBPROC RECOMB ;
ARGU RVX*'TABLE' ;

RV = RVX.'EQEX' ;

Dt = RV.PASDETPS.'DELTAT' ;
Dt = Dt * RV.'ALFA' ;

YH2e  = KCHT $ENTREE SCAL SOMMET RV.INCO.'YH2' ;
YO2e  = KCHT $ENTREE SCAL SOMMET RV.INCO.'YO2' ;
YN2e  = KCHT $ENTREE SCAL SOMMET RV.INCO.'YN2' ;
YH2Oe = KCHT $ENTREE SCAL SOMMET RV.INCO.'YH2O' ;

SENTR1 = 'SOMT' ('DOMA' $ENTREE 'XXDIAGSI') ;
YH2e  = 'SOMT' (YH2e * ('DOMA' $ENTREE 'XXDIAGSI')) / SENTR1 ;
YO2e  = 'SOMT' (YO2e * ('DOMA' $ENTREE 'XXDIAGSI')) / SENTR1 ;
YN2e  = 'SOMT' (YN2e * ('DOMA' $ENTREE 'XXDIAGSI')) / SENTR1 ;
YH2Oe = 'SOMT' (YH2Oe* ('DOMA' $ENTREE 'XXDIAGSI')) / SENTR1 ;
YHEe  = 0.D0 ;
YCO2e = 0.D0 ;
YCOe  = 0.D0 ;
Te = KCHT $ENTREE SCAL SOMMET RV.INCO.'TN' ;
Te = 'SOMT' (Te* ('DOMA' $ENTREE 'XXDIAGSI')) / SENTR1 ;

Pe = RV.INCO.'PM' ;

Tp = RV.INCO.'TPLAQUE' ;

AA1    = RV . 'INCO' . 'A'      ;
BB1    = RV . 'INCO' . 'B'      ;
XH2ON  = RV . 'INCO' . 'XH2ON'  ;
XH2OFF = RV . 'INCO' . 'XH2OFF' ;
BETAI  = RV . 'INCO' . 'BETA'   ;
*IKALP = 1 ;
     YH2s YO2s YN2s YH2Os YHEs YCO2s YCOs Ts TP DEB EFF XH2ONF
     = 'FLAM' 'RECOMBIN'
     YH2e YO2e YN2e YH2Oe YHEe YCO2e YCOe
     Te TP Pe Dt 1 'SIEMENS' AA1 BB1 XH2ON XH2OFF BETAI ;
*     Te TP Pe Dt IKALP 'SIEMENS' AA1 BB1 XH2ON XH2OFF BETAI ;
*YH2s YO2s YN2s YH2Os TS TP DEB EFF XH2ONF
*=  'FLAM' 'RECOMBIN'
*YH2e YO2e YN2e YH2Oe Te TP Pe Dt 1
*'SIEMENS' AA1 BB1 XH2ON XH2OFF BETAI ;

RV . 'INCO' . 'XH2ON' = XH2ONF ;

DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/20 ;
LO = (DD-(20*NN)) EGA 0 ;
SI ( LO ) ;
P_TP = PROG TP ;
RV.INCO.'TP' = RV.INCO.'TP' ET P_TP ;
P_DEB = PROG DEB ;
RV.INCO.'DEBIT' = RV.INCO.'DEBIT' ET P_DEB ;
P_EFF = PROG EFF ;
RV.INCO.'EFFIC' = RV.INCO.'EFFIC' ET P_EFF ;
P_TS = PROG TS ;
RV.INCO.'TS' = RV.INCO.'TS' ET P_TS ;
P_PM = PROG Pe ;
RV.INCO.'PRESSIO' = RV.INCO.'PRESSIO' ET P_PM ;
MAST = SOMT (KOPS RV.INCO.'RHOM' '*' VOL) ;
P_MAST = PROG MAST ;
RV.INCO.'MASSETOT' = RV.INCO.'MASSETOT' ET P_MAST ;

FINSI ;


R = (YH2s*RH2)+(YO2s*RO2)+(YN2s*RN2)+(YH2Os*RH2O);
PSIs = R*TS ;

RV.INCO.'TPLAQUE' = TP ;

DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/20 ;
LO = (DD-(20*NN)) EGA 0 ;
SI ( LO ) ;
MESSAGE 'TEMPERATURE PLAQUE EN K = ' TP ;
MESSAGE 'TEMPERATURE GAZ SORTIE EN K = ' TS ;
MESSAGE 'EFFICACITE RECOMBINEUR EN % = ' EFF ;
MESSAGE 'DEBIT EN KG/S = ' DEB ;
FINSI ;

* MODIFICATION DES CONDITIONS LIMITES

cl_YH2 = 'EXCO' 'YH2' RV.'CLIM' 'YH2' ;
cl_YO2 = 'EXCO' 'YO2' RV.'CLIM' 'YO2' ;
cl_YH2O= 'EXCO' 'YH2O' RV.'CLIM' 'YH2O' ;
cl_GNX = 'EXCO' '1GN' RV.'CLIM' '1GN' ;
cl_GNY = 'EXCO' '2GN' RV.'CLIM' '2GN' ;
cl_KN  = 'EXCO' 'KN' RV.'CLIM' 'KN' ;
cl_EN  = 'EXCO' 'EN' RV.'CLIM' 'EN' ;
cl_PSI = 'EXCO' 'PSIN' RV.'CLIM' 'PSIN' ;

cl_YH2 = 0.D0 * cl_YH2 + YH2s ;
cl_YO2 = 0.D0 * cl_YO2 + YO2s ;
cl_YH2O = 0.D0 * cl_YH2O + YH2Os ;
cl_PSI = 0.D0 * cl_PSI + PSIs ;

cl_GNYm = MANU 'CHPO' (MUR_HOR) 1 '2GN' 0. ;
cl_GNXm = MANU 'CHPO' (MUR_VER) 1 '1GN' 0. ;
cl_GNYe = MANU 'CHPO' (ENTREE) 1 '2GN' (DEB*length) ;
cl_GNXs = MANU 'CHPO' (SORTIE) 1 '1GN' (DEB*length*SIGN) ;

cl_GNX = 0.D0 * cl_GNX + cl_GNXm + cl_GNXs ;
cl_GNY = 0.D0 * cl_GNY + cl_GNYm + cl_GNYe ;

RV.'CLIM' = cl_YH2 + cl_YO2 + cl_YH2O + cl_GNX + cl_GNY +
        cl_KN + cl_EN + cl_PSI ;

detr cl_YH2 ;
detr cl_YO2 ;
detr cl_YH2O ;
detr cl_GNX ;
detr cl_GNY ;
detr cl_KN ;
detr cl_EN ;
detr cl_PSI ;
detr cl_GNYm ;
detr cl_GNXm ;
detr cl_GNYe ;
detr cl_GNXs ;

 as2 ama1  = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1;


*********************************************************************
* PROCEDURE CALCUL1
*********************************************************************
DEBPROC CALCUL1 ;
ARGU RVX*'TABLE' ;

RV = RVX.'EQEX' ;
RVP = RV.'PROJ' ;

RV.PASDETPS.'DELTAT' = 0.5 ;

* filtres sur fractions massiques et calcul de la fraction massique d'azote

YH2  = RV.INCO.'YH2' ;
YO2  = RV.INCO.'YO2' ;
YH2O = RV.INCO.'YH2O' ;
YH2  = KOPS YH2 '|<' 0. ;
YO2  = KOPS YO2 '|<' 0. ;
YH2O = KOPS YH2O '|<' 0. ;

RV.INCO.'YH2' =  YH2 ;
RV.INCO.'YO2'  = YO2 ;
RV.INCO.'YH2O' = YH2O ;
RV.INCO.'YN2'  = KOPS 1.0 '-' RV.INCO.'YH2O' ;
RV.INCO.'YN2'  = KOPS RV.INCO.'YN2' '-' RV.INCO.'YH2' ;
RV.INCO.'YN2'  = KOPS RV.INCO.'YN2' '-' RV.INCO.'YO2' ;

* calcul de la constante des gaz et de la temperature

R = KOPS (KOPS YH2 '*' RH2) '+' (KOPS YO2 '*' RO2) ;
R = KOPS R '+' (KOPS RV.INCO.'YN2' '*' RN2) ;
R = KOPS R '+' (KOPS YH2O '*' RH2O) ;
RV.INCO.'R'  = KCHT $MT 'SCAL' 'SOMMET' R ;
T = KOPS RV.INCO.'PSIN' '/' R ;
RV.INCO.'TN' = KCHT $MT 'SCAL' 'SOMMET' T ;

* calcul des concentrations max et min d'hydrogène

DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/20 ;
LO = (DD-(20*NN)) EGA 0 ;
SI ( LO ) ;
M = KOPS Rg '/' RV.INCO.'R' ;
XH2 = KOPS RV.INCO.'YH2' '*' M ;
XH2 = XH2/MH2 ;
XH2 = XH2*100. ;
XH2MAX = MAXI XH2 ;
XH2MIN = MINI XH2 ;
P_XH2MAX = PROG XH2MAX ;
P_XH2MIN = PROG XH2MIN ;
RV.INCO.'XH2MAX' = RV.INCO.'XH2MAX' ET P_XH2MAX ;
RV.INCO.'XH2MIN' = RV.INCO.'XH2MIN' ET P_XH2MIN ;
FINSI ;

* calcul de rho à partir de la loi d'état des gaz parfaits

RV.INCO.'RHO' = KOPS (RV.INCO.'PM') '/' (RV.INCO.'PSIN') ;
RV.INCO.'RHOC' = NOEL $MT RV.INCO.'RHO' ;
RV.INCO.'RHOC' = KCHT $MT 'SCAL' 'CENTRE' (RV.INCO.'RHOC') ;
RHOM = (SOMT (KOPS RV.INCO.'RHOC' '*' VOL))/VOLT ;
RV.INCO.'RHOM' = RHOM ;

* calcul du terme (rho-rhom)*g de la QDM

DRHO = KOPS (RV.INCO.'RHOC') '-' (RV.INCO.'RHOM') ;
rogx = KCHT $MT SCAL CENTRE 0. ;
rogy = KCHT $MT SCAL CENTRE (-9.81*DRHO) ;
rogx = nomc 'UX' rogx 'NATU' 'DISCRET' ;
rogy = nomc 'UY' rogy 'NATU' 'DISCRET' ;
RV.INCO.'ROG' = KCHT $MT VECT CENTRE (rogx ET rogy) ;

* calcul la vitesse à partir de la qdm

GN = RV.INCO.'GN' ;
gnx= kcht $MT scal sommet (exco 'UX' gn) ;
gny= kcht $MT scal sommet (exco 'UY' gn) ;
unx = kops gnx '/' RV.INCO.'RHO' ;
uny = kops gny '/' RV.INCO.'RHO' ;
unx = nomc 'UX' unx 'NATU' 'DISCRET' ;
uny = nomc 'UY' uny 'NATU' 'DISCRET' ;
un = kcht $MT vect sommet (unx et uny) ;
RV.INCO.'UN' = un ;

* calcul de nut à partir de mut et de rho

RV.INCO.'NUT' = KOPS (RV.INCO.'MUT') '/' (RV.INCO.'RHOC') ;

 as2 ama1  = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1;

*********************************************************************
* PROCEDURE CALCUL2
*********************************************************************
DEBPROC CALCUL2 ;
ARGU RVX*'TABLE' ;

RV = RVX.'EQEX' ;
RVP = RV.'PROJ' ;

* on récupère le pas de temps

Dt = RV.PASDETPS.'DELTAT' ;
Dt = Dt * RV.'ALFA' ;

DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/20 ;
LO = (DD-(20*NN)) EGA 0 ;
SI ( LO ) ;
P_TEMPS = PROG RV.PASDETPS.'TPS' ;
RV.INCO.'TEMPS' = RV.INCO.'TEMPS' ET P_TEMPS ;
FINSI ;

* calcul l'intégrale se trouvant dans le terme Dp/Dt.

ZT = KOPS (RV.INCO.'PSIN') '-' (RV.INCO.'PSI1') ;
ZT = KOPS ZT '/' Dt ;
ZT = KOPS ZT '/' (RV.INCO.'PSIN') ;

ZZ = KOPS ZT '*' RV.INCO.'RHO' ;
ZZ = NOEL $MT ZZ ;
ZZ = KCHT $MT 'SCAL' 'CENTRE' ZZ ;
ZZT = SOMT ( KOPS ZZ '*' VOL ) ;

RV.INCO.'PSI1' = KCHT $MT 'SCAL' 'SOMMET' (RV.INCO.'PSIN') ;

* calcul de dPdt

q0 = 'DBIT' RV.INCO.'GN' $SORTIE ;
q1 = 'DBIT' RV.INCO.'GN' $ENTREE ;
q0 = q0+q1 ;
dPdt = ZZT + q0 ;
dPdt = dPdt*(RV.INCO.'PM')/VOLT/RV.INCO.'RHOM' ;

* terme source de l'équation d'énergie

zuzu = (gamma-1.)/gamma ;

QC = NOEL $MT RV.INCO.'Q' ;
QC = KCHT $MT 'SCAL' 'CENTRE' QC ;
TOTO = KOPS QC '+' dPdt ;
TOTO = KOPS TOTO '/' RV.INCO.'RHOC' ;
TOTO = KOPS TOTO '*' zuzu ;
RV.INCO.'S' = KCHT $MT 'SCAL' 'CENTRE' TOTO ;

* terme source équation de Poisson

Z1 = dPdt/RV.INCO.'PM' ;
ZP = KOPS RV.INCO.'RHOC' '*' Z1 ;
ZS = KOPS ZP '-' ZZ ;
ZZS = KOPS ZS '*' VOL ;
ZZS = ZZS * (-1.0) ;
RV.INCO.'ZS' = KCHT $MT 'SCAL' 'CENTRE' ZZS ;

* calcul de l'évolution de la pression 

RV.INCO.'PM' = RV.INCO.'PM' + (dPdt*dt) ;

* affichage informations toutes les 20 iterations

DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/20 ;
LO = (DD-(20*NN)) EGA 0 ;
SI ( LO ) ;
MESSAGE '==========================================================' ;
MESSAGE 'TEMPS = ' RV.PASDETPS.'TPS' ' PRESSION = ' RV.INCO.'PM' ;
FINSI ;

 as2 ama1  = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1;

***************************************************************
* JEU DE DONNEES NURETH-8 : CALCUL DE DISTRIBUTION/COMBUSTION *
* ALGORITHME 'COMPRESSIBLE FAIBLE MACH' SEMI-IMPLICITE        *
*                    === DISTRIBUTION ===                     *
***************************************************************

opti dime 2 ;
opti elem qua8 ;

*========================================================*
*               GENERATION DU MAILLAGE                   *
*========================================================*

L = 1. ;
H = 2. ;
height= 1. ;
depth = 0.3 ;
length= 0.2 ;

* position du recombineur
XR = 0.2 ;
YR = 0.2 ;
epsi = 0.01 ;
SREC = 'D' ;
DINJ = depth ;

A1 = 0.0 0.0 ;
A2 = L   0.0 ;
A3 = L   H   ;
A4 = 0.0 H   ;

R1 = XR YR ;
R2 = (XR+depth) YR ;
R3 = (XR+depth) (YR+height-depth) ;
R4 = (XR+depth) (YR+height) ;
R5 = XR (YR+height) ;
R6 = XR (YR+height-depth) ;

C1 = XR 0.0 ;
C2 = (XR+depth) 0.0 ;
C3 = L YR ;
C4 = L (YR+height-depth) ;
C5 = L (YR+height) ;
C6 = (XR+depth) H ;
C7 = XR H ;
C8 = 0.0 (YR+height) ;
C9 = 0.0 (YR+height-depth) ;
C10 = 0.0 YR ;

DE1 = 1.2 ;
DE2 = 0.8 ;
DE3 = 0.2 ;

SI (EGA SREC 'G') ;
DE4 = DE3 ;
DE5 = DE3 ;
SINON ;
DE4 = DE3 ;
DE5 = DE3 ;
FINSI ;

A1C1 = 'DROI' A1 C1 'DINI' DE1 'DFIN' DE4 ;
C1R1 = 'DROI' C1 R1 'DINI' DE2 'DFIN' DE3 ;
R1C10 = 'DROI' R1 C10 'DINI' DE4 'DFIN' DE1 ;
C10A1 = 'DROI' C10 A1 'DINI' DE3 'DFIN' DE2 ;

C1C2 = 'DROI' C1 C2 'DINI' DE3 'DFIN' DE3 ;
C2R2 = 'DROI' C2 R2 'DINI' DE2 'DFIN' DE3 ;
R2R1 = 'DROI' R2 R1 'DINI' DE3 'DFIN' DE3 ;
R1C1 = 'DROI' R1 C1 'DINI' DE3 'DFIN' DE2 ;

C2A2 = 'DROI' C2 A2 'DINI' DE5 'DFIN' DE1 ;
A2C3 = 'DROI' A2 C3 'DINI' DE2 'DFIN' DE3 ;
C3R2 = 'DROI' C3 R2 'DINI' DE1 'DFIN' DE5 ;
R2C2 = 'DROI' R2 C2 'DINI' DE3 'DFIN' DE2 ;

R2C3 = 'DROI' R2 C3 'DINI' DE5 'DFIN' DE1 ;
C3C4 = 'DROI' C3 C4 'DINI' DE3 'DFIN' DE3 ;
C4R3 = 'DROI' C4 R3 'DINI' DE1 'DFIN' DE5 ;
R3R2 = 'DROI' R3 R2 'DINI' DE3 'DFIN' DE3 ;

R3C4 = 'DROI' R3 C4 'DINI' DE5 'DFIN' DE1 ;
C4C5 = 'DROI' C4 C5 'DINI' DE3 'DFIN' DE3 ;
C5R4 = 'DROI' C5 R4 'DINI' DE1 'DFIN' DE5 ;
R4R3 = 'DROI' R4 R3 'DINI' DE3 'DFIN' DE3 ;

R4C5 = 'DROI' R4 C5 'DINI' DE5 'DFIN' DE1 ;
C5A3 = 'DROI' C5 A3 'DINI' DE3 'DFIN' DE1 ;
A3C6 = 'DROI' A3 C6 'DINI' DE1 'DFIN' DE5 ;
C6R4 = 'DROI' C6 R4 'DINI' DE1 'DFIN' DE3 ;

R5R4 = 'DROI' R5 R4 'DINI' DE3 'DFIN' DE3 ;
R4C6 = 'DROI' R4 C6 'DINI' DE3 'DFIN' DE1 ;
C6C7 = 'DROI' C6 C7 'DINI' DE3 'DFIN' DE3 ;
C7R5 = 'DROI' C7 R5 'DINI' DE1 'DFIN' DE3 ;

C8R5 = 'DROI' C8 R5 'DINI' DE1 'DFIN' DE4 ;
R5C7 = 'DROI' R5 C7 'DINI' DE3 'DFIN' DE1 ;
C7A4 = 'DROI' C7 A4 'DINI' DE4 'DFIN' DE1 ;
A4C8 = 'DROI' A4 C8 'DINI' DE1 'DFIN' DE3 ; 

C9R6 = 'DROI' C9 R6 'DINI' DE1 'DFIN' DE4 ;
R6R5 = 'DROI' R6 R5 'DINI' DE3 'DFIN' DE3 ;
R5C8 = 'DROI' R5 C8 'DINI' DE4 'DFIN' DE1 ;
C8C9 = 'DROI' C8 C9 'DINI' DE3 'DFIN' DE3 ;

C10R1 = 'DROI' C10 R1 'DINI' DE1 'DFIN' DE4 ;
R1R6 = 'DROI' R1 R6 'DINI' DE3 'DFIN' DE3 ;
R6C9 = 'DROI' R6 C9 'DINI' DE4 'DFIN' DE1 ;
C9C10 = 'DROI' C9 C10 'DINI' DE3 'DFIN' DE3 ;


COMPA = 'DALLER' A1C1 C1R1 R1C10 C10A1 'PLAN' ;
COMPB = 'DALLER' C1C2 C2R2 R2R1 R1C1 'PLAN' ;
COMPC = 'DALLER' C2A2 A2C3 C3R2 R2C2 'PLAN' ;
COMPD = 'DALLER' R2C3 C3C4 C4R3 R3R2 'PLAN' ;
COMPE = 'DALLER' R3C4 C4C5 C5R4 R4R3 'PLAN' ;
COMPF = 'DALLER' R4C5 C5A3 A3C6 C6R4 'PLAN' ;
COMPG = 'DALLER' R5R4 R4C6 C6C7 C7R5 'PLAN' ;
COMPH = 'DALLER' C8R5 R5C7 C7A4 A4C8 'PLAN' ;
COMPI = 'DALLER' C9R6 R6R5 R5C8 C8C9 'PLAN' ;
COMPJ = 'DALLER' C10R1 R1R6 R6C9 C9C10 'PLAN' ;

COMPA = 'COUL' COMPA 'ROUG' ;
COMPB = 'COUL' COMPB 'BLEU' ;
COMPC = 'COUL' COMPC 'JAUN' ;
COMPD = 'COUL' COMPD 'TURQ' ;
COMPE = 'COUL' COMPE 'BLAN' ;
COMPF = 'COUL' COMPF 'ROUG' ;
COMPG = 'COUL' COMPG 'VERT' ;
COMPH = 'COUL' COMPH 'BLAN' ;
COMPI = 'COUL' COMPI 'ROUG' ;
COMPJ = 'COUL' COMPJ 'BLEU' ;

MT = COMPA ET COMPB ET COMPC ET 
     COMPD ET COMPE ET COMPF ET COMPG ET COMPH ET COMPI ET COMPJ ;
MMT = Chan MT 'QUAF';
'ELIM' epsi MT ;
WALL = A1C1 ET C1C2 ET C2A2 ET A2C3 ET C3C4 ET C4C5 ET C5A3 ET
       A3C6 ET C6C7 ET C7A4 ET A4C8 ET C8C9 ET C9C10 ET C10A1 ;

$MT = 'MODE' MMT 'NAVIER_STOKES' DISCR ;
MT  = 'DOMA' $MT 'MAILLAGE' ;
CNT =  CONT MT ;


ENTREE = R2R1 ;
Mentree= 'CHAN' entree 'QUAF';
$ENTREE = 'MODE' MENTREE 'NAVIER_STOKES' DISCR ;
entree  = 'DOMA' $entree 'MAILLAGE';
SI (EGA SREC 'G') ;
SORTIE = R6R5 ;
BOUCHEE = R4R3 ;
SIGN = -1.0 ;
SINON ;
SORTIE = R4R3 ;
BOUCHEE = R6R5 ;
SIGN = 1.0 ;
FINSI ;
Msortie = 'CHAN' sortie 'QUAF';
$SORTIE = 'MODE' MSORTIE 'NAVIER_STOKES' DISCR ;
sortie  = 'DOMA' $sortie 'MAILLAGE';

Mwall = 'CHAN' wall 'QUAF';
$WALL = 'MODE' MWALL 'NAVIER_STOKES' DISCR ;
wall  = 'DOMA' $wall 'MAILLAGE';

VOL = 'DOMA' $MT 'VOLUME' ;
VOLT = 'SOMT' VOL ;
VOL = KCHT $MT 'SCAL' 'CENTRE' VOL ;

ELIM (Mmt et Mentree et Msortie et Mwall) epsi;

'MESSAGE' 'VOLUME TOTAL = ' VOLT ;

******************************************
* FRONTIERES POUR CONDITIONS AUX LIMITES *
******************************************
 
MUR_HOR = A1C1 ET C1C2 ET C2A2 ET A3C6 ET C6C7 ET C7A4 ET R5R4 ;
MUR_VER = A2C3 ET C3C4 ET C4C5 ET C5A3 ET A4C8 ET C8C9 
                ET C9C10 ET C10A1 ET R3R2 ET R1R6 ET BOUCHEE ;

*========================================================*
*               DONNEES DU CALCUL DE DISTRIBUTION        *
*========================================================*

*--------------------------------------
* temperature et pression de saturation
*--------------------------------------
T = 373.0 ;
ZZ = 7000.*(-1.)/T ;
PSAT = 1.055E21/(T ** 5.)* (EXP ZZ) ;
PSAT = PSAT*101325.0 ;

*--------------------------------------
* masses molaires et constantes des gaz
*--------------------------------------
MH2 = 2.0E-3 ;
MO2 = 32.0E-3 ;
MN2 = 28.0E-3 ;
MH2O = 18.0E-3 ;

Rg = 8.313 ;

RH2 = Rg/MH2 ;
RO2 = Rg/MO2 ;
RN2 = Rg/MN2 ;
RH2O = Rg/MH2O ;

*-------------------
* conditions air sec  
*-------------------
P1 = 1.E5 ;
T1 = 273.0 ;
XH2_SEC = 0.1 ;
XAIR_SEC = 1.0 - XH2_SEC ;

*-----------------------------
* pressions partielles des gaz
* fractions volumiques
*-----------------------------

PVAP = PSAT ;
PAIR = (XAIR_SEC*P1)*(T/T1) ;
PH2 = (XH2_SEC*P1)*(T/T1) ;
P = PVAP + PAIR + PH2 ;

XH2 = PH2/P ;
XAIR = PAIR/P ;
XH2O = PVAP/P ;
XO2 = 0.21*XAIR ;
XN2 = 1.0 - XH2 -XH2O - XO2 ;

*---------------------------------------------
* masse molaire du melange et constante du gaz
*---------------------------------------------
M = (XH2*MH2) + (XO2*MO2) + (XN2*MN2) + (XH2O*MH2O) ;
R = Rg/M ;

*------------------------------------------------------------
* densite du melange (loi d'etat GP)
*------------------------------------------------------------
PSI = R*T ;
RHO = P/PSI ;

*------------------------
* Caractéristiques du jet
*------------------------

XH2Jet = 0.0 ;
XH2OJet = 0.0 ;
XO2Jet = 0.21*( 1.0 - XH2Jet - XH2OJet ) ;
XN2Jet = 1.0 - XH2Jet - XH2OJet - XO2Jet ;
MJet = (XH2Jet*MH2) + (XH2OJet*MH2O) + (XO2Jet*MO2) +
       (XN2Jet*MN2) ;
RJet = Rg/MJet ;
TJet = 650.0 ;
PSIJet = Rjet*Tjet ;
RhoJet = P/PSIJet ;
GJet = 2.5 ;
GJetabs= GJet ;
UJet = GJet/RhoJet ;
KJet = 0.05*UJet*UJet ;
EJet = 0.02*(UJet**3.)/DINJ ;
Gjet = Gjet*SIGN ;

*------------------------------------------------
* fractions massiques des differents constituants
*------------------------------------------------
YH2 = XH2*MH2/M ;
YO2 = XO2*MO2/M ;
YN2 = XN2*MN2/M ;
YH2O = XH2O*MH2O/M ;
YH2Jet = XH2Jet*MH2/MJet ;
YH2OJet =XH2OJet*MH2O/MJet ;
YO2Jet = XO2Jet*MO2/MJet ;
YN2Jet =XN2Jet*MN2/MJet ;

*------------------------------------------------------------------
* Cp des différents constituants et du mélange (supposés constants)
*------------------------------------------------------------------
CPH2 = 13514.04 + (1.684537*T) ;
CPH2O = 1715.632 + (0.552805*T) ;
CPN2 = 1006.15 + (0.1387166*T) ;
CPO2 = 907.580 + (0.1420522*T) ;

CP = (YH2*CPH2) + (YH2O*CPH2O) + (YO2*CPO2) + (YN2*CPN2) ;

gamma = CP/(CP-R) ;

*----------------------------------------------------
* calcul d'une vitesse de référence (pour les tracés)
*----------------------------------------------------

uref = gamma*R*T ;
uref = uref**0.5 ;
uref = 0.01*uref ;
ampl = 1./uref ;

*------------------------------------------------------------------------
* Prandtl, Schmidt, viscosité (loi de Sutherland), conductivité thermique
* coefficients de diffusion moléculaire  
*------------------------------------------------------------------------
Pr = 0.7 ;
Prt = 1.0 ;
Sct = 1.0 ;

Mu = 1.716E-5*((T/273.0)**1.5)*(273.0+110.5)/(T+110.5) ;
Nu = Mu/RHO ;
Lambda = Mu*CP/Pr ;

alpha = Lambda/RHO/CP ;

DH2O = 2.55E-5 ;
DH2 = 7.12E-5 ;
DO2 = 2.06E-5 ;

PSIWALL = T*R ;
h = 10. ;

U0 = 1.0D0 ;
L0 = 1.0D0 ;

*==============================================================
* TABLE EQEX CONTENANT LES EQUATIONS A RESOUDRE :           
*  0) FILTRE K-E
*  1) QUANTITE DE MOUVEMENT     --- OPERATEUR NSKE
*  2) EQUATION DE L'ENERGIE     --- OPERATEUR TSCAL
*  3) TRANSPORT DES ESPECES     --- OPERATEUR TSCAL
* AINSI QUE 2 PROCEDURES (CALCUL1 ET CALCUL2) POUR L'ALGORITHME
* 'FAIBLE MACH'
*  4) CONDITIONS AUX LIMITES
*==============================================================

RV = 'EQEX' $MT 'ITMA' 100 'TFINAL' 10. 'ALFA' 0.8
 'ZONE' $MT 'OPER'      'CALCUL1'
 'ZONE' $MT 'OPER'      'FILTREKE' 'UN' L0 NU 'INCO' 'KN' 'EN'
 'OPTI' 'CONS'
 'ZONE' $MT 'OPER'      'NSKE' 'RHOC' MU 'MUT' 'UN' 'ROG'
                                'INCO' 'GN' 'KN' 'EN'
 'OPTI' 'NOCONS'
 'ZONE' $MT 'OPER'      'TSCAL' alpha 'UN' 'S' 'NUT' Prt 
                                'INCO' 'PSIN'
 'OPTI' 'EFM1' 'EXPL'
 'ZONE' $WALL 'OPER'    'ECHI' h PSIWALL 'INCO' 'PSIN'  ;

RV = 'EQEX' RV
 'ZONE' $MT 'OPER'      'TSCAL' DH2 'UN' 0.0 'NUT' SCT 
                                'INCO' 'YH2'
 'ZONE' $MT 'OPER'      'TSCAL' DO2 'UN' 0.0 'NUT' SCT 
                                'INCO' 'YO2'
 'ZONE' $MT 'OPER'      'TSCAL' DH2O 'UN' 0.0 'NUT' SCT 
                                'INCO' 'YH2O' 
 'ZONE' $MT 'OPER'      'CALCUL2' 
 'ZONE' $MT 'OPER'      'RECOMB' ;

  RV=EQEX RV OPTI EFM1 'CENTREE'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'GN' 'DELTAT' 'INCO' 'GN'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'PSIN' 'DELTAT' 'INCO' 'PSIN'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'YH2' 'DELTAT' 'INCO' 'YH2'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'YO2' 'DELTAT' 'INCO' 'YO2'
  ;
  RV=EQEX RV OPTI EFM1 'CENTREE'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'YH2O' 'DELTAT' 'INCO' 'YH2O'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
  ;
 RV = 'EQEX' RV
 'CLIM' 'GN' 'UIMP' MUR_VER 0.0
        'GN' 'VIMP' MUR_HOR 0.0
        'GN' 'UIMP' SORTIE  GJet
        'GN' 'VIMP' SORTIE  0.0
        'GN' 'UIMP' ENTREE  0.0
        'GN' 'VIMP' ENTREE  GJetabs ;
RV = 'EQEX' RV
 'CLIM' 'PSIN' 'TIMP' SORTIE  PSIJet 
        'YH2' 'TIMP' SORTIE YH2Jet
        'YO2' 'TIMP' SORTIE YO2Jet 
        'YH2O' 'TIMP' SORTIE YH2OJet 
        'KN' 'TIMP' SORTIE  KJet
        'EN' 'TIMP' SORTIE  EJet ;

*===========================================================
* TABLE EQEX POUR L'EQUATION DE POISSON (SOLVEUR ELLIPTIQUE)

RVP = EQEX 'OPTI' 'EF' KPRESS
 'ZONE' $MT    OPER  KBBT  -1. betap  INCO 'GN' 'PRES'
           'OPTI' 'EF' 'CENTREE' 'INCOD' KPRESS
 'ZONE' $MT    OPER  FIMP  'ZS'       INCO 'PRES'
;

    rvp.'METHINV'.TYPINV=1 ;
    rvp.'METHINV'.IMPINV=0 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;

  RV.'PROJ'     =RVP ;

*===========================================================

*==========================================================
* TABLE INCO CONTENANT LES INCONNUES ET DONNEES DU PROBLEME 
*==========================================================

  RV.'NOMVI' = 'GN' ;

  RV.'INCO' = 'TABLE' 'INCO' ;
  RVP.'INCO' = RV.'INCO' ;

* vitesse et quantité de mouvement
  RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.);
  RV.INCO.'GN' = kcht $MT VECT SOMMET (0. 0.);
* pression             
  RV.INCO.'PRES' = kcht $MT SCAL KPRESS 0.;
* température
  RV.INCO.'PSIN' = kcht $MT SCAL SOMMET   PSI ;
  RV.INCO.'PSI1' = kcht $MT SCAL SOMMET   PSI ;
  RV.INCO.'TN'   = kcht $MT SCAL SOMMET   T ;
  RV.INCO.'TPLAQUE' = T ;
* variables modèle de turbulence
  RV.INCO.'EN' = kcht $MT SCAL SOMMET   1.E-5 ;
  RV.INCO.'KN' = kcht $MT SCAL SOMMET   1.E-5 ;
  RV.INCO.'MUT' = kcht $MT SCAL CENTRE  Mu ;
  RV.INCO.'NUT' = kcht $MT SCAL CENTRE  Nu ;
* masse volumique
  RV.INCO.'RHO' = kcht $MT SCAL SOMMET  RHO    ;
  RV.INCO.'RHOC'= kcht $MT SCAL CENTRE  RHO    ;
  RV.INCO.'RHOM'= RHO ;
* fractions massiques
  RV.INCO.'YH2' = kcht $MT SCAL SOMMET  YH2   ;
  RV.INCO.'YO2' = kcht $MT SCAL SOMMET  YO2   ;
  RV.INCO.'YN2' = kcht $MT SCAL SOMMET  YN2   ;
  RV.INCO.'YH2O'= kcht $MT SCAL SOMMET  YH2O   ;
* termes sources (Poisson/énergie/QDM)
  RV.INCO.'S'=kcht $MT SCAL CENTRE   0.    ;
  RV.INCO.'Q'=kcht $MT SCAL SOMMET   0.    ;
  RV.INCO.'ROG'=kcht $MT VECT CENTRE   (0. 0.)     ;
  RV.INCO.'ZS' = KCHT $MT SCAL CENTRE 0.0 ;
* pression
  RV.INCO.'PM' = P ;

* progressions
  RV.INCO.'TEMPS' = PROG 0. ;
  RV.INCO.'XH2MAX' = PROG XH2 ;
  RV.INCO.'XH2MIN' = PROG XH2 ;
  RV.INCO.'TP'    = PROG RV.INCO.'TPLAQUE' ;
  RV.INCO.'DEBIT' = PROG 0. ;
  RV.INCO.'EFFIC' = PROG 0. ;
  RV.INCO.'TS'    = PROG T ;

MASINI = KOPS RV.INCO.'RHOC' '*' VOL ;
MASINI = SOMT MASINI ;

  RV.INCO.'PRESSIO' = PROG P ;
  RV.INCO.'MASSETOT' = PROG MASINI ;

RV . 'INCO' . 'A'      = 0.48D-8 ;
RV . 'INCO' . 'B'      = 0.58D-3 ;
RV . 'INCO' . 'XH2ON'  = 0.005 ;
RV . 'INCO' . 'XH2OFF' = 0.005 ;
RV . 'INCO' . 'BETA'   = 0.2 ;

exec rv ;

SI (RV.PASDETPS.'TPS' < 8.0) ;
        ERREUR 5 ;
FINSI ;
SI (MAXI(RV.INCO.'DEBIT') < 4.9D-2);
        ERREUR 5 ;
FINSI ;
SI (MAXI(RV.INCO.'TP') < 690.D0);
        ERREUR 5 ;
FINSI ;
SI (MAXI(RV.INCO.'TN') < 550.D0);
        ERREUR 5 ;
FINSI ;

FIN ;