* fichier : basmachQ.dgibi
*******************************************************************
* JEU DE DONNEES POUR TESTER LE BON FONCTIONNEMENT DES OPERATEURS *
* NSKE, PRESSION, TSCAL. INJECTION D'UN MELANGE DE GAZ DANS UN *
* VOLUME FERME ET PRESSURISATION. ALGORITHME 'FAIBLES MACH' POUR *
* LES EQUATIONS NAVIER-STOKES FAIBLEMENT COMPRESSIBLES + MODELE *
* DE TURBULENCE K-E. ON TESTE SUR LA CONSERVATION DE LA MASSE *
* H. PAILLERE SEMT/TTMF AVRIL 1997 *
* FORMULATION EF ELEMENTS QUA8 *
*******************************************************************
*************************************************************
* PROCEDURE PROCINJ *
*************************************************************
DEBPROC PROCINJ ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
$MT = RV.'DOMAINE' ;
* recuperation des anciennes conditions aux 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' ;
* modification des conditions aux limites :
RHOJet = RV.INCO.'PM'/PSIJET ;
UJet = GJet/RHOJet ;
KJet = 0.05*UJet*UJet ;
EJet = 0.02*(UJet**3.)/DINJ ;
cl_KN = 0.D0 * cl_KN + KJet ;
cl_EN = 0.D0 * cl_EN + EJet ;
* reimposition des conditions aux limites :
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 ;
*tc mise en commentaire du finsi ci-dessous
*FINSI ;
FINPROC ;
*********************************************************************
* PROCEDURE CALCUL1
*********************************************************************
DEBPROC CALCUL1 ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
RVP = RV.'PRESSION' ;
$MT = RV.'DOMAINE' ;
RV.PASDETPS.'DELTAT' = 0.01 ;
* 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 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 (RV.'DOMAINE') SCAL CENTRE 0. ;
rogy = KCHT (RV.'DOMAINE') SCAL CENTRE (-9.81*DRHO) ;
rogx = nomc 'UX' rogx 'NATU' 'DISCRET' ;
rogy = nomc 'UY' rogy 'NATU' 'DISCRET' ;
RV.INCO.'ROG' = KCHT (RV.'DOMAINE') VECT CENTRE (rogx ET rogy) ;
* calcul la vitesse à partir de la qdm
GN = RV.INCO.'GN' ;
gnx= kcht (rv.'DOMAINE') scal sommet (exco 'UX' gn) ;
gny= kcht (rv.'DOMAINE') 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 (RV.'DOMAINE') 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') ;
FINPROC ;
*********************************************************************
* PROCEDURE CALCUL2
*********************************************************************
DEBPROC CALCUL2 ;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
RVP = RV.'PRESSION' ;
$MT = RV.'DOMAINE' ;
* on récupère le pas de temps
Dt = RV.PASDETPS.'DELTAT' ;
Dt = Dt * RV.'ALFA' ;
* 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' $INJ ;
RV.INCO.'DEBIT' = q0 ;
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 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' ;
MESSAGE 'MASSE INI =' MASINI ' MASSE TOTALE = ' (RV.INCO.'RHOM'*VOLT) ;
PT = PROG RV.'PASDETPS'.'TPS' ;
PP = PROG RV.INCO.'PM' ;
PMTO = PROG (RV.INCO.'RHOM'*VOLT) ;
RV.INCO.'TEMPS' = RV.INCO.'TEMPS' ET PT ;
RV.INCO.'PRESS' = RV.INCO.'PRESS' ET PP ;
RV.INCO.'MTOTA' = RV.INCO.'MTOTA' ET PMTO ;
MESSAGE '==========================================================' ;
FINSI ;
FINPROC ;
***************************************************************
* JEU DE DONNEES NURETH-8 : CALCUL DE DISTRIBUTION/COMBUSTION *
* ALGORITHME 'COMPRESSIBLE FAIBLE MACH' SEMI-IMPLICITE *
* === DISTRIBUTION === *
***************************************************************
opti dime 2 ;
opti elem qua8 ;
opti isov suli ;
*========================================================*
* GENERATION DU MAILLAGE *
*========================================================*
L_A = 4.0 ;
H_A = 4.0 ;
L_B = 4.0 ;
H_B = H_A ;
HINJ = 1.0 ;
DINJ = 1.0 ;
HAB = 1. ;
DAB = 2.0 ;
LAB = 1. ;
DX = 1. ;
DY = 1. ;
epsi = (DX/10.) ;
******************
* COMPARTIMENT A *
******************
A7 = 0. 0. ;
A8 = HINJ 0. ;
A1 = (HINJ+DINJ) 0. ;
A2 = L_A 0.0 ;
A3 = L_A HAB ;
A4 = L_A (HAB+DAB) ;
A5 = L_A H_A ;
A6 = 0. H_A ;
N78 = 'ENTIER' (HINJ/DX) ;
N81 = 'ENTIER' (DINJ/DX) ;
N12 = 'ENTIER' ((L_A-HINJ-DINJ)/DX) ;
N23 = 'ENTIER' (HAB/DY) ;
N34 = 'ENTIER' (DAB/DY) ;
N45 = 'ENTIER' ((H_A-HAB-DAB)/DY) ;
N56 = N78 + N81 + N12 ;
N67 = N23 + N34 + N45 ;
A7A8 = A7 'DROI' N78 A8 ;
A8A1 = A8 'DROI' N81 A1 ;
A1A2 = A1 'DROI' N12 A2 ;
A2A3 = A2 'DROI' N23 A3 ;
A3A4 = A3 'DROI' N34 A4 ;
A4A5 = A4 'DROI' N45 A5 ;
A5A6 = A5 'DROI' N56 A6 ;
A6A7 = A6 'DROI' N67 A7 ;
BAS_A = A7A8 ET A8A1 ET A1A2 ;
DRO_A = A2A3 ET A3A4 ET A4A5 ;
HAU_A = A5A6 ;
GAU_A = A6A7 ;
COMP_A = 'DALLER' BAS_A DRO_A HAU_A GAU_A 'PLAN' ;
******************
* COMPARTIMENT B *
******************
B1 = (L_A+LAB) 0. ;
B2 = (L_A+LAB+L_B) 0. ;
B3 = (L_A+LAB+L_B) H_B ;
B6 = (L_A+LAB) H_B ;
B7 = (L_A+LAB) (HAB+DAB) ;
B8 = (L_A+LAB) HAB ;
N23 = 'ENTIER' (H_B/DY) ;
N36 = 'ENTIER' (L_B/DX) ;
N67 = 'ENTIER' ((H_B-HAB-DAB)/DY) ;
N78 = 'ENTIER' (DAB/DY) ;
N81 = N23 - N67 - N78 ;
N12 = N36 ;
B1B2 = B1 'DROI' N12 B2 ;
B2B3 = B2 'DROI' N23 B3 ;
B3B6 = B3 'DROI' N36 B6 ;
B6B7 = B6 'DROI' N67 B7 ;
B7B8 = B7 'DROI' N78 B8 ;
B8B1 = B8 'DROI' N81 B1 ;
BAS_B = B1B2 ;
DRO_B = B2B3 ;
HAU_B = B3B6 ;
GAU_B = B6B7 ET B7B8 ET B8B1 ;
COMP_B = 'DALLER' BAS_B DRO_B HAU_B GAU_B 'PLAN' ;
***************
* JONCTION AB *
***************
NX = 'ENTIER' (LAB/DX) ;
NY = 'ENTIER' (DAB/DY) ;
BAS_AB = A3 'DROI' NX B8 ;
DRO_AB = B8 'DROI' NY B7 ;
HAU_AB = B7 'DROI' NX A4 ;
GAU_AB = A4 'DROI' NY A3 ;
J_AB = 'DALLER' BAS_AB DRO_AB HAU_AB GAU_AB 'PLAN' ;
********************
* MAILLAGE COMPLET *
********************
COMP_A = COMP_A 'COUL' 'ROUG' ;
COMP_B = COMP_B 'COUL' 'VERT' ;
J_AB = J_AB 'COUL' 'BLAN' ;
MT = COMP_A ET COMP_B ET J_AB ;
'ELIM' epsi MT ;
'ORIE' MT ;
CNT = 'CONT' MT ;
******************
* TABLES DOMAINE *
******************
$MT = 'DOMA' MT epsi 'MACRO' ;
$INJ = 'DOMA' BAS_A 'MACRO' 'INCL' $MT epsi ;
$COMP_A = DOMA COMP_A 'MACRO' 'INCL' $MT epsi ;
$COMP_B = DOMA COMP_B 'MACRO' 'INCL' $MT epsi ;
$J_AB = DOMA J_AB 'MACRO' 'INCL' $MT epsi ;
MT = $MT.MAILLAGE ;
VOL = 'DOMA' $MT 'VOLUME' ;
VOLT = 'SOMT' VOL ;
VOL = KCHT $MT 'SCAL' 'CENTRE' VOL ;
'MESSAGE' 'VOLUME TOTAL = ' VOLT ;
******************************************
* FRONTIERES POUR CONDITIONS AUX LIMITES *
******************************************
MUR_HOR = A7A8 ET A1A2 ET BAS_AB ET B1B2 ET B3B6
ET HAU_AB ET A5A6 ;
MUR_VER = A2A3 ET B8B1 ET B2B3 ET B6B7 ET A4A5 ET A6A7 ;
BRECHE = A8A1 ;
MUR_HOR = MUR_HOR 'COUL' 'VERT' ;
MUR_VER = MUR_VER 'COUL' 'ROUG' ;
BRECHE = BRECHE 'COUL' 'JAUNE' ;
*========================================================*
* 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.0 ;
XH2O = 0.1 ;
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 = 1.0E5 ;
T = 400.0 ;
PSI = R*T ;
RHO = P/PSI ;
*------------------------
* Caractéristiques du jet
*------------------------
XH2Jet = 0.1 ;
XH2OJet = 0.3 ;
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 = 10.*DINJ/(DINJ+(DX*0.5)) ;
UJet = GJet/RhoJet ;
KJet = 0.05*UJet*UJet ;
EJet = 0.02*(UJet**3.)/DINJ ;
*------------------------------------------------
* 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 ;
*==============================================================
* 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' 30 'ALFA' 0.9
'ZONE' $MT 'OPER' 'PROCINJ'
'ZONE' $MT 'OPER' 'CALCUL1'
'ZONE' $MT 'OPER' 'FILTREKE' UJet DINJ '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' ;
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' ;
RV = 'EQEX' RV
'CLIM' 'GN' 'UIMP' MUR_VER 0.0
'GN' 'VIMP' MUR_HOR 0.0
'GN' 'UIMP' BRECHE 0.0
'GN' 'VIMP' BRECHE GJet
'KN' 'TIMP' BRECHE KJet
'EN' 'TIMP' BRECHE EJet ;
RV = 'EQEX' RV
'CLIM' 'PSIN' 'TIMP' BRECHE PSIJet
'YH2' 'TIMP' BRECHE YH2Jet
'YO2' 'TIMP' BRECHE YO2Jet
'YH2O' 'TIMP' BRECHE YH2OJet ;
*===========================================================
* 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. 0.);
RV.INCO.'GN' = kcht $MT VECT SOMMET (0. 0.);
* température
RV.INCO.'PSIN' = kcht $MT SCAL SOMMET PSI ;
RV.INCO.'PSI1' = kcht $MT SCAL SOMMET PSI ;
* 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 contenant certaines données
RV.INCO.'TEMPS' = PROG 0. ;
RV.INCO.'PRESS' = PROG P ;
RV.INCO.'MTOTA' = PROG (RHO*VOLT) ;
* debit de masse
RV.INCO.'DEBIT' = 0. ;
MASINI = KOPS RV.INCO.'RHOC' '*' VOL ;
MASINI = SOMT MASINI ;
exec rv ;
MTHEO = MASINI + (RV.INCO.'DEBIT' * RV.PASDETPS.'TPS') ;
MESSAGE MTHEO ;
SI ( (MAXI RV.INCO.'MTOTA') < (0.95*MTHEO) ) ;
ERREUR 5 ;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer