* fichier : flamhms.dgibi
GRAPH = FAUX ;
***********************************************************
* COMBUSTION EN REGIME LAMINAIRE - MODELE TRAVIS *
* ON TESTE LES VALEURS DE PRESSION ET DE TEMPERATURE AICC *
***********************************************************
opti dime 2 ;
opti elem qua8 ;
opti isov suli ;
nx = 8 ;
ny = 8 ;
p1 = -1. 0 ;
p2 = 1. 0 ;
p3 = 1. 2.0 ;
p4 = -1. 2.0 ;
bas = p1 d nx p2 ;
dro = p2 d ny p3 ;
hau = p3 d nx p4 ;
gau = p4 d ny p1 ;
MT = 'DALLER' BAS DRO HAU GAU 'PLAN' ;
***************************
* MODELES 'NAVIER_STOKES' *
***************************
MMT = CHAN MT QUAF ;
$MT = 'MODE' MMT 'NAVIER_STOKES' 'MACRO' ;
DOMA $MT 'IMPR' ;
mt = doma $mt maillage ;
CNT = CONT MT ;
WALL = (BAS ET DRO ET HAU ET GAU) CHAN QUAF ;
ELIM (WALL ET MMT) 1.E-4;
$WALL = 'MODE' WALL 'NAVIER_STOKES' 'MACRO' ;
VOL = 'DOMA' $MT 'VOLUME' ;
VOLT = 'SOMT' VOL ;
**********************
* LIEU DE L'ALLUMAGE *
**********************
PA = (DOMA $MT SOMMET) POIN 'PROC' (0. 0.2) ;
*********************************************************************
* PROCEDURE CALCUL1
*********************************************************************
DEBPROC CALCUL1 ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
RVP = RV.'PRESSION' ;
SI ( RV.PASDETPS.'NUPASDT' < 30 ) ;
RV.PASDETPS.'DELTAT' = 0.001 ;
FINSI ;
* 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 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 ;
as2 ama1 = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1 ;
*********************************************************************
* PROCEDURE CALCUL2
*********************************************************************
DEBPROC CALCUL2 ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
RVP = RV.'PRESSION' ;
* on récupère le pas de temps
Dt = RV.PASDETPS.'DELTAT' ;
Dt = Dt * RV.'ALFA' ;
* combustion et mise a jour du cp et (gamma-1)/gamma
YH2 YO2 YN2 YH2O Q = 'FLAM' 'LAMINAIR'
(RV.INCO.'RHO') (RV.INCO.'CV') 'LINEAIRE' (RV.INCO.'R')
(RV.INCO.'YH2') (RV.INCO.'YO2') (RV.INCO.'YN2')
(RV.INCO.'YH2O') (RV.INCO.'TN') Dt ;
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 ;
RV.INCO.'Q' = KCHT $MT SCAL SOMMET Q ;
RC = NOEL $MT RV.INCO.'R' ;
CVC = NOEL $MT RV.INCO.'CV' ;
CPC = KOPS RC '+' CVC ;
RV.INCO.'alpha' = KOPS lambda '/' RV.INCO.'RHOC' ;
RV.INCO.'alpha' = KOPS RV.INCO.'alpha' '/' CPC ;
RV.INCO.'zuzu' = KOPS RC '/' CPC ;
* 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') ;
* Evolution de la pression moyenne
dPdt = ZZT ;
dPdt = dPdt*(RV.INCO.'PM')/VOLT/RV.INCO.'RHOM' ;
MESSAGE 'dPdt = ' dPdt ;
* terme source de l'équation d'énergie
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 '*' RV.INCO.'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 5 iterations
DD = RV.PASDETPS.'NUPASDT' ;
NN = DD/5 ;
LO = (DD-(5*NN)) EGA 0 ;
SI ( LO ) ;
MESSAGE '==========================================================' ;
MESSAGE 'TEMPS = ' RV.'PASDETPS'.'TPS' 'PRESSION = ' RV.INCO.'PM'
'MIN TEMP = ' (MINI RV.INCO.'TN') 'MAX TEMP = ' (MAXI RV.INCO.'TN') ;
MESSAGE 'MASSE INITIALE = ' MASINI ' MASSE ACTUELLE = '
(RV.INCO.'RHOM'*VOLT) ;
PP = PROG RV.INCO.'PM' ;
PT = PROG RV.PASDETPS.'TPS' ;
RV.INCO.'TEMPS' = RV.INCO.'TEMPS' ET PT ;
RV.INCO.'PRESS' = RV.INCO.'PRESS' ET PP ;
MESSAGE '==========================================================' ;
FINSI ;
as2 ama1 = 'KOPS' 'MATRIK' ;
FINPROC as2 ama1 ;
*========================================================*
* DONNEES DU CALCUL DE DISTRIBUTION *
*========================================================*
*--------------------------------------
* 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 ;
*----------------------------------------------------------
* fractions molaires des differents constituants du melange
*----------------------------------------------------------
XH2 = 0.075 ;
XH2O = 0. ;
XO2 = 0.21*(1.0 - XH2 - XH2O) ;
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 ;
*------------------------------------------------------------
* pression, temperature et densite du melange (loi d'etat GP)
*------------------------------------------------------------
P = 0.9E5 ;
T = 400.0 ;
PSI = R*T ;
RHO = P/PSI ;
PSIWALL = PSI ;
*------------------------------------------------
* fractions massiques des differents constituants
*------------------------------------------------
YH2 = XH2*MH2/M ;
YO2 = XO2*MO2/M ;
YN2 = XN2*MN2/M ;
YH2O = XH2O*MH2O/M ;
*------------------------------------------------------------------
* 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 ;
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 ;
DIF = 1.E-9 ;
*==============================================================
* 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' 5000 'ALFA' 0.6 'TFINAL' 10.0
'ZONE' $MT 'OPER' 'CALCUL1'
'OPTI' 'CONS'
'ZONE' $MT 'OPER' 'NS' MU 'UN' 'ROG' 'INCO' 'GN'
'OPTI' 'NOCONS'
'ZONE' $MT 'OPER' 'TSCAL' 'alpha' 'UN' 'S' 'INCO' 'PSIN'
;
RV = 'EQEX' RV
'OPTI' 'NOCONS'
'ZONE' $MT 'OPER' 'TSCAL' DIF 'UN' 0.0 'INCO' 'YH2'
'OPTI' 'NOCONS'
'ZONE' $MT 'OPER' 'TSCAL' DIF 'UN' 0.0 'INCO' 'YO2'
'OPTI' 'NOCONS'
'ZONE' $MT 'OPER' 'TSCAL' DIF 'UN' 0.0 'INCO' 'YH2O'
;
RV = 'EQEX' RV
'ZONE' $MT 'OPER' 'CALCUL2'
'OPTI' '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'
'ZONE' $MT 'OPER' 'DFDT' 1. 'YH2O' 'DELTAT' 'INCO' 'YH2O'
;
RV = 'EQEX' RV
'CLIM' 'GN' 'UIMP' (DRO ET GAU) 0.0
'GN' 'VIMP' (HAU ET BAS) 0.0 ;
*===========================================================
* TABLE EQPR POUR L'EQUATION DE POISSON (SOLVEUR ELLIPTIQUE)
*===========================================================
RVP = 'EQPR' $MT KTYPI 1
'ZONE' $MT 'OPER' 'PRESSION' 'ZS'
'PIMP' 0. ;
*==========================================================
* TABLE INCO CONTENANT LES INCONNUES ET DONNEES DU PROBLEME
*==========================================================
RV.'PRESSION' = RVP ;
RV.'NOMVI' = 'GN' ;
RV.'INCO' = 'TABLE' 'INCO' ;
RVP.'INCO' = RV.'INCO' ;
* vitesse et quantité de mouvement
RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 1.E-5);
RV.INCO.'GN' = kcht $MT VECT SOMMET (0. 1.E-5);
* 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.'R' = kcht $MT SCAL SOMMET R ;
RV.INCO.'CV' = kcht $MT SCAL SOMMET (CP/gamma) ;
* 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 ;
* diffusivite thermique
RV.INCO.'alpha' = KCHT $MT SCAL CENTRE alpha ;
RV.INCO.'zuzu' = KCHT $MT SCAL CENTRE ((gamma-1)/gamma) ;
* pression
RV.INCO.'PM' = P ;
* progressions contenant certaines données
RV.INCO.'TEMPS' = PROG 0. ;
RV.INCO.'PRESS' = PROG P ;
RV.INCO.'MTOTA' = PROG (RHO*VOLT) ;
RV.INCO.'YH2u' = RV.INCO.'YH2' ;
M = KOPS (RV.INCO.'YH2'/MH2) '+' (RV.INCO.'YO2'/MO2) ;
M = KOPS M '+' (RV.INCO.'YN2'/MN2) ;
M = KOPS M '+' (RV.INCO.'YH2O'/MH2O) ;
M = KOPS 1. '/' M ;
XH2 = KOPS RV.INCO.'YH2' '*' M ;
XH2 = XH2/MH2 ;
XH2 = XH2*100. ;
CA = MANU 'CHPO' PA 1 SCAL 3100. ;
TA = KCHT $MT 'SCAL' 'SOMMET' CA ;
RV.INCO.'TN' = (RV.INCO.'TN' ET TA) ;
SI ( GRAPH ) ;
trace rv.inco.'TN' MT 14 'TITR' 'temperature initiale [K]' ;
FINSI ;
RV.INCO.'PSIN' = KOPS RV.INCO.'R' '*' RV.INCO.'TN' ;
RV.INCO.'RHO' = KOPS RV.INCO.'PM' '/' RV.INCO.'PSIN' ;
RV.INCO.'RHOC' = NOEL $MT RV.INCO.'RHO' ;
VOL = DOMA $MT 'VOLUME' ;
VOLT = SOMT VOL ;
MESSAGE 'VOLUME TOTAL = ' VOLT ;
RHOM = SOMT (KOPS RV.INCO.'RHOC' '*' VOL) ;
RHOM = RHOM/VOLT ;
RHO = RHOM ;
MASINI = KOPS RV.INCO.'RHOC' '*' VOL ;
MASINI = SOMT MASINI ;
MESSAGE 'MASSE INITIALE = ' MASINI ;
RV.'TFINAL' = 1.0 ;
exec rv ;
SI ( GRAPH ) ;
M = KOPS (RV.INCO.'YH2'/MH2) '+' (RV.INCO.'YO2'/MO2) ;
M = KOPS M '+' (RV.INCO.'YN2'/MN2) ;
M = KOPS M '+' (RV.INCO.'YH2O'/MH2O) ;
M = KOPS 1. '/' M ;
XH2 = KOPS RV.INCO.'YH2' '*' M ;
XH2 = XH2/MH2 ;
XH2 = XH2*100. ;
TRACE XH2 MT CNT 14 'TITR' 'FRACTION VOL. HYDROGENE' ;
FINSI ;
RV.'TFINAL' = 3.0 ;
exec rv ;
SI ( GRAPH ) ;
TAB1 = TABLE ;
TAB1.1 = 'MARQ LOSA REGU' ;
TAB1.'TITRE' = TABLE ;
TAB1.'TITRE' . 1 = 'MOT' 'PRESSURE' ;
EVOL1 = 'EVOL' 'MANU' 'TIME [s]' RV.INCO.'TEMPS'
'Press. [Pa]' RV.INCO.'PRESS' ;
'DESS' EVOL1 'LEGE' TAB1 'MIMA' ;
FINSI ;
PM = RV.INCO.'PM' ;
TM = NOEL $MT RV.INCO.'TN' ;
TM = SOMT( KOPS TM '*' VOL) / VOLT ;
ERRP = (PM-2.574E5)/2.574E5 ;
ERRP = 100.*(ABS ERRP) ;
ERRT = (TM-1188.)/1188. ;
ERRT = 100.*(ABS ERRT) ;
MESSAGE 'PAICC CALCULEE = ' PM 'ERREUR = ' (ERRP) '%' ;
MESSAGE 'TAICC CALCULEE = ' TM 'ERREUR = ' (ERRT) '%' ;
SI ( ERRP > 2. ) ;
ERREUR 5 ;
FINSI ;
SI ( ERRT > 2. ) ;
ERREUR 5 ;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer