* fichier : injection.dgibi
************************************************************************
* NOM : Injection Air/Air (Air chaud dans Air froid sur 6
* secondes) dans une cavité carrée 2D plan
*
* DESCRIPTION : Cas test du modèle asymptotique à bas nombre de Mach
* pour un écoulement laminaire, confiné, monoespèce mais
* faisant intervenir la thermique
*
*
*
* LANGAGE : GIBIANE-CAST3M
* AUTEUR : Etienne STUDER (CEA/DEN/DM2S/SFME/LTMF)
* mél : mp1@semt2.smts.cea.fr
**********************************************************************
* VERSION : v1, 15/12/2005, version initiale
* HISTORIQUE : v1, ??/??/2003, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Toutes les procédures sont incluses dans le fichier
*
*
*
************************************************************************
interact= FAUX ;
'OPTION' 'DIME' 2 'ELEM' QUA8 'ISOV' 'SULI' ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PS' ;
'SINON' ;
'OPTION' 'TRAC' 'X' ;
'FINSI' ;
************************************************************************
* PROCEDURES LOCALES
************************************************************************
* Calcul d'un terme source VECTEUR Explicite
'DEBPROC' TSOUR ;
'ARGUMENT' RVX * 'TABLE' ;
RV = RVX . 'EQEX' ;
LSTINC = RVX.'LISTINCO' ;
NOMV = 'EXTRAIRE' LSTINC 1 ;
ROG = RVX.'ARG1' ;
RGX = 'EXCO' (RV.'INCO'. ROG) 'UX' ;
RGY = 'EXCO' (RV.'INCO'. ROG) 'UY' ;
RGX = RGX '*' ('DOMA' (RVX.'DOMZ') 'XXDIAGSI') ;
RGY = RGY '*' ('DOMA' (RVX.'DOMZ') 'XXDIAGSI') ;
RGX = 'NOMC' RGX '1UN' 'NATURE' 'DISCRET' ;
RGY = 'NOMC' RGY '2UN' 'NATURE' 'DISCRET' ;
ROG = RGX 'ET' RGY ;
AS AF = 'KOPS' 'MATRIK' ;
'FINPROC' ROG AF ;
'DEBP' EXEC ;
*
************************ K E X IMPL ************************************
*
* Résolution implicite d'un problème de mécanique des fluides décrit
* par l'intermédiaire d'une table structurée par l'opérateur EQEX.
*
* E/ rv : Table décrivant le problème à résoudre
* /S rv : A l'indice INCO, table contenant la solution au temps final
*
* L'utilisateur peut par l'intermédiaire de procédures personnelles
* stocker les résultats aux temps intermédiaires, faire des tracés, ...
*
************************ K E X E C ************************************
*
'ARGUMENT' rv*'TABLE ' ;
*
*========================
* TESTS et INITIALISATION
*========================
*
*
* On n'accepte que les formulations en modèle NAVIER_STOKES
'SI' ('EGA' (rv . 'NAVISTOK') 0) ;
'MESS' 'Pas de modèle NAVIER_STOKES' ;
'ERRE' 5 ;
'FINSI' ;
*
dimope = 'DIME' (rv . 'LISTOPER') ;
*
* Coefficient de relaxation mis à 1 si non initialisé
'SI' ('NON' ('EXISTE' rv 'OMEGA')) ;
'MESS' '** WARNING ** OMEGA NON défini --> OMEGA=1' ;
omeg = 1.D0 ;
'SINON' ;
omeg = rv . 'OMEGA' ;
'FINSI' ;
*
* Création de la table pour historique
'SI' ('NON' ('EXISTE' rv 'HIST')) ;
rv . 'HIST' = 'TABLE' ;
'FINSI' ;
*
* IMPKRES : Niveau d'impression pour KRES
* IMPTCRR : Fréquence d'impression pour TCRR (affichage résidu)
IMPKRES = 0 ;
IMPTCRR = RV . 'IMPR' ;
*
*========================
* RESOLUTION suivant EQEX
*========================
*
*
* 1
* Boucle en temps
*----------------
ITMA = rv . 'ITMA' ;
'SI' ('<EG' ITMA 1) ; ITMA = 1 ; 'FINSI' ;
'REPETER' bloc1 ITMA ;
*
* 2
* Boucle de point fixe interne à un pas de temps
*-----------------------------------------------
'REPETER' bloci (rv . 'NITER') ;
* st mat : Second membre et matrice des opérateurs DFDT
* sf mau : Second membre er matrice des opérateurs sauf DFDT
st mat = 'KOPS' 'MATRIK' ;
sf mau = 'KOPS' 'MATRIK' ;
*
* 3
* Boucle sur les opérateurs de EQEX
*----------------------------------
'REPETER' bloc2 dimope ;
nomper = 'EXTRAIRE' &bloc2 (rv . 'LISTOPER') ;
notable = 'MOT' ('TEXTE' ('CHAINE' &bloc2 nomper)) ;
'SI' ('EGA' nomper 'DFDT ') ;
msi mai = ('TEXTE' nomper) (rv . notable) ;
mat = mat 'ET' mai ;
st = st 'ET' msi ;
'SINO' ;
msi mai = ('TEXTE' nomper) (rv . notable) ;
mau = mau 'ET' mai ;
sf = sf 'ET' msi ;
'FINSI' ;
'FIN' bloc2 ;
*--------------------------------------
* Fin Boucle sur les opérateurs de EQEX
* 3
*
s2 = sf 'ET' st ;
ma1 = mau 'ET' mat ;
* ES preconditionnement MATRICE
'SI' ('EXISTE' rv 'ESPR') ;
ma1 = ma1 'ET' (rv.'ESPR') ;
'FINSI' ;
*
'SI' ('EXISTE' rv 'CLIM') ;
s1 = rv . 'CLIM' ;
'SINON' ;
s1 = ' ' ;
'FINSI' ;
rv . 'S2' = s2 ;
rv . 'METHINV' . 'MATASS' = ma1 ;
rv . 'METHINV' . 'MAPREC' = ma1 ;
res = 'KRES' ma1 'TYPI' (rv . 'METHINV')
'CLIM' s1
'SMBR' s2
'IMPR' IMPKRES ;
eps = 'TCRR' res omeg (rv . 'INCO') 'IMPR' IMPTCRR ;
'MENAGE' ;
'FIN' bloci ;
*---------------------------------------------------
* Fin Boucle de point fixe interne à un pas de temps
* 2
*
irt = 0 ;
'SI' ('EGA' (rv . 'ITMA') 0) ;
irt = 'TCNM' rv 'NOUP';
'SINON' ;
irt = 'TCNM' rv ;
'FINSI' ;
'MENAGE' ;
'SI' ('EGA' irt 1) ;
'MESSAGE' ' Temps final atteint : ' (rv . 'PASDETPS' . 'TPS') ;
'QUITTER' bloc1 ;
'FINSI' ;
'FIN' bloc1 ;
*--------------------
* Fin Boucle en temps
* 1
*
************************ K E X IMPL ************************************
'FINPROC' ;
*************************************************************************
* DEFINITION DU PROBLEME
*
*************************************************************************
* le raffinement: nombre de points dans une dimension
* en 1D la verticale et en 2D on prend le même nombre
* dans les deux directions
RAF = 20 ;
* Choix de la méthode d'inversion
* TYPI = 1 : méthode directe
* TYPI = 3 : méthode itérative
TYPI = 1 ;
* Le cas test peut être traité en 1D ou 2D plan
UND = VRAI ;
DEUXD = FAUX ;
* Dimensions géométriques
* L : largeur de l'injection au bas de la cavité (centrée)
* H : hauteur de la cavité
* LL : largeur de la cavité (=L en 1D)
L = 0.2 ;
H = 7.0 ;
'SI' (UND) ;
LL = L ;
'SINON' ;
'SI' (DEUXD) ;
LL = 3.0 ;
'SINON' ;
'MESSAGE' '==> Cas imprévu !!!!!!!!' ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
* Caractéristiques du gaz
* GAMA : coeffcient isentropique
* RAIR : constante du gaz (R/M) avec R constante des gaz et M la masse
* molaire
* PRANDTL : nombre de Prandtl
GAMA = 1.4 ;
RAIR = 287.0 ;
PRANDTL = 0.71 ;
* Conditions initiales
* P0 : pression (Pa)
* T0 : température (K)
* R0 : masse volumique (kg/m3)
P0 = 1.E5 ;
T0 = 300.0 ;
R0 = P0 '/' (RAIR '*' T0) ;
* conditions à l'injection
* MPH : débit surfacique kg/m²/s
* TH : température (K)
* V0 : vitesse à l'injection (m/s)
MPH = 1.0 ;
TH = 600.0 ;
V0 = MPH '*' RAIR '*' TH '/' P0 ;
* Choix du cas de calcul par l'intermédiaire de booléens
* C1 : Reynolds = 40 sans prise en compte de la gravité
* C2 : idem C1 mais avec ajout de la flottabilité
* C3 : idem C2 mais avec Reynolds = 800
C1 = VRAI ;
C2 = FAUX ;
C3 = FAUX ;
'SI' (C1) ;
REYN = 40.0 ;
GRAV = 0.0 ;
'FINSI' ;
'SI' (C2) ;
REYN = 40.0 ;
GRAV = 9.81 ;
'FINSI' ;
'SI' (C3) ;
REYN = 800.0 ;
GRAV = 9.81 ;
'FINSI' ;
* Calcul des propriétés physiques
* MU : viscosité kg/m/s
* LAMBDA : conductivité thermique W/m/K
* CPAIR et CVAIR : chaleur spécifique J/kg/K
MU = MPH '*' L '/' REYN ;
LAMBDA = MU '*' GAMA '*' RAIR '/' PRANDTL '/' (GAMA '-' 1.0) ;
CPAIR = GAMA '*' RAIR '/' (GAMA '-' 1.0) ;
CVAIR = CPAIR '-' RAIR ;
********************************************************************
* LE MAILLAGE
********************************************************************
LS2 = L '/' 2.0 ;
LLS2 = LL '/' 2.0 ;
LLS4 = LL '/' 4.0 ;
HS4 = H '/' 4.0 ;
HS2 = H '/' 2.0 ;
3HS4 = H '*' 3.0 '/' 4.0 ;
'DENS' (H '/' ('FLOTTANT' RAF)) ;
'SI' (UND) ;
P1 = (LS2 '*' (-1.0)) 0.0 ;
P2 = 0.0 0.0 ;
P3 = (LS2) 0.0 ;
P1P2 = 'DROIT' 1 P1 P2 ;
P2P3 = 'DROIT' 1 P2 P3 ;
P4 = 0.0 H ;
P5 = (LS2 '*' (-1.0)) H ;
P6 = LS2 H ;
P2P4 = 'DROIT' P2 P4 ;
P1P5 = 'DROIT' P1 P5 ;
P3P6 = 'DROIT' P3 P6 ;
P6P4 = 'DROIT' 1 P6 P4 ;
P4P5 = 'DROIT' 1 P4 P5 ;
BAS = P1P2 'ET' P2P3 ;
HAUT = P6P4 'ET' P4P5 ;
LMIL = P2P4 ;
CDRO = P3P6 ;
CGAU = 'INVERSE' P1P5 ;
VTF = BAS P3P6 HAUT ('INVERSE' P1P5) 'DALLER' 'PLAN' ;
'FINSI' ;
'SI' (DEUXD) ;
P1 = (LLS2 '*' (-1.0)) 0.0 ;
P2 = 0. 0. ;
P1A = (LLS4 '*' (-1.0)) 0.0 ;
P1B = (LS2 '*' (-1.0)) 0.0 ;
N1 = 'ENTIER' (LS2 '/' LL '*' RAF) ;
N2 = 'ENTIER' ((LLS4 '-' LS2) '/' LL '*' RAF) ;
N3 = ('ENTIER' (RAF '/' 2.0)) '-' N1 '-' N2 ;
P1P1A = 'DROIT' N3 P1 P1A ;
P1AP1B = 'DROIT' N2 P1A P1B ;
P1BP2 = 'DROIT' N1 P1B P2 ;
P3 = LLS2 0.0 ;
P3B = LLS4 0.0 ;
P3A = LS2 0.0 ;
P2P3A = 'DROIT' N1 P2 P3A ;
P3AP3B = 'DROIT' N2 P3A P3B ;
P3BP3 = 'DROIT' N3 P3B P3 ;
INJE = P1BP2 'ET' P2P3A ;
BAS = P1P1A 'ET' P1AP1B 'ET' P3AP3B 'ET' P3BP3 ;
NY = 'ENTIER' (RAF '/' 4.0) ;
P4 = 0.0 H ;
HAUT = 'INVERSE' ('PLUS' (BAS 'ET' INJE) P4) ;
P5 = 'POIN' HAUT 'PROC' ('PLUS' P1 P4 ) ;
P6 = 'POIN' HAUT 'PROC' ('PLUS' P3 P4) ;
CDRO = 'DROIT' ('ENTIER' RAF) P3 P6 ;
CGAU = 'DROIT' ('ENTIER' RAF) P5 P1 ;
VTF = (BAS 'ET' INJE) CDRO HAUT CGAU 'DALLER' 'PLAN' ;
'FINSI' ;
* Choix des éléments Vitesse/Pression
DISCR = 'QUAF' ;
KPRES = 'CENTREP1' ;
* Choix du décentrement des termes convectifs
KSUPG = 'CENTREE' ;
* Choix de la condensation ou non de la matrice masse
*KMASS = 'EF' ;
KMASS = 'EFM1' ;
* Choix de l'ordre du schéma en temps
* BDF1 : Euler
* BDF2 : schéma à 3 pas (ordre 2)
KBDF = 'BDF1' ;
*KBDF = 'BDF2' ;
* positionnement des constantes pour le schéma en temps BDF1 (EULER) ou
* BDF2
'SI' ('EGA' KBDF 'BDF1') ;
AN = 1.0 ;
ANM1 = (-1.0) ;
ANM2 = 0.0 ;
ADT = 1.0 ;
'SINON' ;
'SI' ('EGA' KBDF 'BDF2') ;
AN = 3.0 ;
ANM1 = (-4.0) ;
ANM2 = 1.0 ;
ADT = 2.0 ;
'SINON' ;
'MESSAGE' '==> ERREUR indice KBDF' ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
* Precision maillage
EPS0 = 1.D-7 ;
* Definition des objets modèles
MVTF = 'CHANGER' VTF 'QUAF' ;
$VTF = 'MODELISER' MVTF 'NAVIER_STOKES' DISCR ;
* Traitement de la geometrie
*-----------------------------------------------------------
* Creation de l'enveloppe du maillage FLUIDE
GEO = TABLE ;
GEO.'MVTF' = MVTF ;
GEO.'$VTF' = $VTF ;
GEO.'EPSI' = EPS0 ;
GEO.'MHAUT' = HAUT ;
GEO.'$HAUT' = 'MODE' GEO.'MHAUT' 'NAVIER_STOKES' DISCR ;
TOTO = 'DOMA' GEO.'$HAUT' 'CENTRE' ;
GEO.'MBAS' = BAS ;
GEO.'$BAS' = 'MODE' GEO.'MBAS' 'NAVIER_STOKES' DISCR ;
TOTO1 = 'DOMA' GEO.'$BAS' 'CENTRE' ;
'SI' (DEUXD) ;
GEO.'MINJE' = INJE ;
GEO.'$INJE' = 'MODE' GEO.'MINJE' 'NAVIER_STOKES' DISCR ;
TOTO2 = 'DOMA' GEO.'$INJE' 'CENTRE' ;
'FINSI' ;
GEO.'MCDRO' = CDRO ;
GEO.'$CDRO' = 'MODE' GEO.'MCDRO' 'NAVIER_STOKES' DISCR ;
TUTU1 = 'DOMA' GEO.'$CDRO' 'CENTRE' ;
GEO.'MCGAU' = CGAU ;
GEO.'$CGAU' = 'MODE' GEO.'MCGAU' 'NAVIER_STOKES' DISCR ;
TUTU2 = 'DOMA' GEO.'$CGAU' 'CENTRE' ;
'SI' (UND) ;
GEO.'MMENVF' = CGAU 'ET' BAS 'ET' CDRO 'ET' HAUT ;
'FINSI' ;
'SI' (DEUXD) ;
GEO.'MMENVF' = CGAU 'ET' BAS 'ET' CDRO 'ET' HAUT 'ET' INJE ;
'FINSI' ;
GEO.'$MENVF' = 'MODE' GEO.'MMENVF' 'NAVIER_STOKES' DISCR ;
TUTU3 = 'DOMA' GEO.'$MENVF' 'CENTRE' ;
DOMF = 'DOMA' GEO.'$VTF' 'FACE' ;
'SI' (UND) ;
'ELIMINATION' DOMF (TOTO 'ET' TOTO1
'ET' TUTU1 'ET' TUTU2 'ET' TUTU3) EPS0 ;
'FINSI' ;
'SI' (DEUXD) ;
'ELIMINATION' DOMF (TOTO 'ET' TOTO1 'ET' TOTO2
'ET' TUTU1 'ET' TUTU2 'ET' TUTU3) EPS0 ;
'FINSI' ;
*
* Récupération des maillages a partir des MODELES
GEO.'BAS' = 'DOMA' GEO.'$BAS' 'MAILLAGE' ;
GEO.'HAUT' = 'DOMA' GEO.'$HAUT' 'MAILLAGE' ;
GEO.'VTF' = 'DOMA' GEO.'$VTF' 'MAILLAGE' ;
GEO.'CGAU' = 'DOMA' GEO.'$CGAU' 'MAILLAGE' ;
GEO.'CDRO' = 'DOMA' GEO.'$CDRO' 'MAILLAGE' ;
'SI' (DEUXD) ;
GEO.'INJE' = 'DOMA' GEO.'$INJE' 'MAILLAGE' ;
'FINSI' ;
GEO.'MENVF' = 'DOMA' GEO.'$MENVF' 'MAILLAGE' ;
* Maillage pour la condition aux limites en pression
* Il faut un point à l'intérieur
LVTF = 'CHAN' 'LINEAIRE' GEO.'VTF' ;
PI0 = LVTF 'POIN' 'PROC' ( 0. H ) ;
MELTI = (GEO.'MVTF') 'ELEM' 'CONTENANT' PI0 ;
$ELTI = 'MODE' MELTI 'NAVIER_STOKES' DISCR ;
MTPI = 'DOMA' $ELTI KPRES ;
MTPI = 'ELEM' ('LECT' 1) MTPI ;
'ELIM' (MTPI 'ET' ('DOMA' GEO.'$VTF' KPRES )) EPS0 ;
* Maillage pour ouvrir en haut pour faire fuir et equilibrer
* (test de la pertinence de la contrainte de divergence
GEO.'HAUT1' = 'CHAN' GEO.'HAUT' 'POI1' ;
'SI' (UND) ;
GEO.'HAUT1' = 'DIFF' GEO.'HAUT1' ('ELEM' GEO.'HAUT1' 2) ;
'FINSI' ;
'SI' (DEUXD) ;
GEO.'HAUT1' = 'DIFF' GEO.'HAUT1' ('ELEM' GEO.'HAUT1' ('ENTIER' (RAF
'/' 2.0))) ;
'FINSI' ;
* Fonction de forme en 2D pour l'injection: profil parabolique
*
'SI' (DEUXD) ;
XX = 'COORDONNEE' 1 GEO.'INJE' ;
FX = XX '*' XX '*' (-1.0) '+' (L '**' 2.0 '/' 4.0) '*' (6.0 '/' (L
'*' L)) ;
'FINSI' ;
* MTPI = 'ELEM' ('LECT' 1) ('DOMA' GEO.'$VTF' KPRES) ;
* Definition du modèle physique table RV
*
* - Definition de constantes intervenant dans les equations
* FCPRECI = Frequence dans le methode iterative
* FCPRECT= idem
FCPRECI = 1 ;
FCPRECT = 1 ;
*---------------------------------------------------------------
*- EQUATIONS DE LA QDM
* Formulation non conservative
* Decentrement KSUPG
*
RV = 'EQEX' 'NITER' 1 'OMEGA' 1.0 'ITMA' 0
'OPTI' 'EF' 'IMPL' KSUPG
'ZONE' (GEO.'$VTF') 'OPER' 'LAPN' 'MU' 'INCO' 'UN'
'OPTI' 'EF' 'IMPL' KSUPG
'ZONE' (GEO.'$VTF') 'OPER' 'KONV' 'RHO' 'UN' 'MU' 'INCO' 'UN'
;
RV = 'EQEX' RV
'OPTI' 'EF' 'CENTREE' 'IMPL'
'ZONE' (GEO.'$VTF') 'OPER' TSOUR 'ROG' 'INCO' 'UN'
;
'SI' ('EGA' KBDF 'BDF2') ;
RV = 'EQEX' RV
'OPTI' KMASS 'CENTREE' 'BDF2'
'ZONE' (GEO.'$VTF') 'OPER' 'DFDT' 'RHO' 'UNM' 'UNMM' 'DT'
'INCO' 'UN'
;
'SINON' ;
RV = 'EQEX' RV
'OPTI' KMASS 'CENTREE'
'ZONE' (GEO.'$VTF') 'OPER' 'DFDT' 'RHO' 'UNM' 'DT' 'INCO' 'UN'
;
'FINSI' ;
* Conditions aux limites QDM
'SI' (UND) ;
RV = 'EQEX' RV 'CLIM'
'UN' 'UIMP' GEO.'CDRO' 0.
'UN' 'UIMP' GEO.'CGAU' 0.
'UN' 'UIMP' GEO.'HAUT' 0.
'UN' 'VIMP' GEO.'HAUT1' 0.
'UN' 'UIMP' GEO.'BAS' 0.
'UN' 'VIMP' GEO.'BAS' V0
;
'FINSI' ;
'SI' (DEUXD) ;
RV = 'EQEX' RV 'CLIM'
'UN' 'UIMP' GEO.'CDRO' 0.
'UN' 'VIMP' GEO.'CDRO' 0.
'UN' 'UIMP' GEO.'CGAU' 0.
'UN' 'VIMP' GEO.'CGAU' 0.
'UN' 'UIMP' GEO.'BAS' 0.
'UN' 'VIMP' GEO.'BAS' 0.
;
RV = 'EQEX' RV 'CLIM'
'UN' 'UIMP' GEO.'HAUT' 0.
'UN' 'VIMP' GEO.'HAUT1' 0.
'UN' 'UIMP' GEO.'INJE' 0.
'UN' 'VIMP' GEO.'INJE' (FX '*' V0)
;
'FINSI' ;
* Inconnues
RV.'INCO' = TABLE 'INCO' ;
RV.'INCO'.'UN'= 'KCHT' GEO.'$VTF' 'VECT' 'SOMMET' (0.0 0.0) ;
RV.'INCO'.'ROG'= 'KCHT' GEO.'$VTF' 'VECT' 'SOMMET' (0.0 0.0) ;
RV.'INCO'.'MU'= MU ;
RV.'INCO'.'UNM' = 'COPIER' RV.'INCO'.'UN' ;
RV.'INCO'.'UNMM' = 'COPIER' RV.'INCO'.'UN' ;
*
* Gestion du transitoire
* TFIN : temps final en secondes
* CCFL : condition CFL utilisée
*
TFIN = 6.0 ;
CCFL = 2.0 ;
'SI' (DEUXD) ;
DT0 = CCFL '*' R0 '*' (LL '/' ('FLOTTANT' RAF)) '/' MPH ;
'FINSI' ;
'SI' (UND) ;
DT0 = CCFL '*' R0 '*' (H '/' ('FLOTTANT' RAF)) '/' MPH ;
'FINSI' ;
NPAS = 'ENTIER' (TFIN '/' DT0) ;
RV.'INCO'.'DT' = ('FLOTTANT' (TFIN '/' NPAS)) ;
'MESSAGE' 'Pas de temps choisi DT = ' RV.'INCO'.'DT' ;
RV.'METHINV'.TYPINV = TYPI ;
RV.'METHINV'.TYPRENU = 'SLOANE' ;
RV.'METHINV'.IMPINV = 1 ;
RV.'METHINV'.NITMAX = 3000 ;
RV.'METHINV'.PRECOND = 5 ;
RV.'METHINV'.SCALING = 1 ;
RV.'METHINV'.OUBMAT = 1 ;
RV.'METHINV'.ILUTLFIL = 2.0 ;
RV.'METHINV'.'FCPRECI' = FCPRECI ;
RV.'METHINV'.'FCPRECT' = FCPRECT ;
RV.'METHINV'.RESID = 1.e-15 ;
RV = 'EQEX' RV
'OPTI' 'EF' KSUPG 'IMPL' KPRES
'ZONE' (GEO.'$VTF') 'OPER' 'KBBT' 'CF' 'INCO' 'UN' 'PRES'
;
RV.'INCO'.'CF' = 1.0 ;
* Terme source pour la divergence: DSRC
* volume fermé: on impose la pression en un point MTPI
* ou on ouvre le contour en 1 point
RV = 'EQEX' RV
'OPTI' 'EF' 'IMPL' 'INCOD' KPRES
'ZONE' (GEO.'$VTF') 'OPER' 'FIMP' 'DSRC' 'INCO' 'PRES'
;
RV.'INCO'.'DSRC' = 0.0 ;
*RV = 'EQEX' RV
* 'CLIM' 'PRES' 'TIMP' MTPI 0.
*;
*-------------------------------------------------------------
* Definition de l'inconnue de PRESSION
*
RV.'INCO'.'PRES' = 'KCHT' GEO.'$VTF' 'SCAL' KPRES 0.0 ;
*------------------------------------------------------------
* Equation en temperature: TABLE RTF
*------------------------------------------------------------
RT = 'EQEX' 'ITMA' 0
'OPTI' 'EF' 'IMPL' KSUPG 'NOCONS'
'ZONE' (GEO.'$VTF') 'OPER' 'LAPN' 'ALPHA' 'INCO' 'TF'
'OPTI' 'EF' 'IMPL' KSUPG 'NOCONS'
'ZONE' (GEO.'$VTF') 'OPER' 'KONV' 'RHO' 'UN' 'ALPHA' 'INCO' 'TF'
'OPTI' 'EF' 'CENTREE' 'IMPL'
'ZONE' (GEO.'$VTF') 'OPER' 'FIMP' 'STF' 'INCO' 'TF'
;
'SI' ('EGA' KBDF 'BDF2') ;
RT = 'EQEX' RT
'OPTI' KMASS 'CENTREE' 'BDF2'
'ZONE' (GEO.'$VTF') 'OPER' 'DFDT' 'RHO' 'TFNM' 'TFNMM' 'DT'
'INCO' 'TF' ;
'SINON' ;
RT = 'EQEX' RT
'OPTI' KMASS 'CENTREE'
'ZONE' (GEO.'$VTF') 'OPER' 'DFDT' 'RHO' 'TFNM' 'DT'
'INCO' 'TF' ;
'FINSI' ;
'SI' (UND) ;
RT = 'EQEX' RT
'CLIM' 'TF' 'TIMP' (GEO.'BAS') TH
;
'FINSI' ;
'SI' (DEUXD) ;
RT = 'EQEX' RT
'CLIM' 'TF' 'TIMP' (GEO.'INJE') TH
;
'FINSI' ;
* Table des inconnues
RT.'INCO' = RV.'INCO' ;
'SI' (DEUXD) ;
TIN = 'MANUEL' 'CHPO' GEO.'INJE' 'SCAL' TH ;
'FINSI' ;
'SI' (UND) ;
TIN = 'MANUEL' 'CHPO' GEO.'BAS' 'SCAL' TH ;
'FINSI' ;
RV.'INCO'.'TF' = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' T0 TIN ;
RV.'INCO'.'TFNM' = 'COPIER' RV.'INCO'.'TF' ;
RV.'INCO'.'TFNMM' = 'COPIER' RV.'INCO'.'TF' ;
RV.'INCO'.'STF' = 0.0 ;
RV.'INCO'.'ALPHA' = MU '/' PRANDTL ;
RT.'METHINV'.TYPINV = TYPI ;
RT.'METHINV'.TYPRENU = 'SLOANE' ;
RT.'METHINV'.IMPINV = 1 ;
RT.'METHINV'.NITMAX = 400 ;
RT.'METHINV'.PRECOND = 5 ;
RT.'METHINV'.SCALING = 1 ;
RT.'METHINV'.OUBMAT = 1 ;
RT.'METHINV'.ILUTLFIL = 1.5 ;
RT.'METHINV'.'FCPRECI' = FCPRECI ;
RT.'METHINV'.'FCPRECT' = FCPRECT ;
RT.'METHINV'.RESID = 1.e-15 ;
* Execution - Algorithme BAS-MACH
* Conditions initiales
RV.'INCO'.'PM' = 'PROG' P0 ;
RV.'INCO'.'PE' = 'PROG' P0 ;
*
RV.'INCO'.'RHOM' = 'PROG' R0 ;
RV.'INCO'.'REM' = 'PROG' (R0 '*' CVAIR '*' T0) ;
RHO = ('INVERSE' (RV.'INCO'.'TF' '*' RAIR ) ) '*' P0 ;
RV.'INCO'.'RHO' = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' RHO ;
RV.'INCO'.'RHONM' = 'COPIER' RV.'INCO'.'RHO' ;
RV.'INCO'.'RHONMM' = 'COPIER' RV.'INCO'.'RHO' ;
'MESSAGE' '==> Masse initiale (kg) : ' ('SOMT' (RV.'INCO'.'RHO'
'*' ('DOMA' (GEO.'$VTF') 'XXDIAGSI'))) ;
DIAG = 'DOMA' (GEO.'$VTF') 'XXDIAGSI' ;
RV.'INCO'.'LTPS' = 'PROG' 0.0 ;
RV.'INCO'.'ERRORM' = 'PROG' 0.0 ;
RV.'INCO'.'ERRORE' = 'PROG' 0.0 ;
'OPTION' 'ECHO' 1 ;
VTOTAL = 'SOMT' ('DOMA' GEO.'$VTF' 'XXDIAGSI') ;
'MESSAGE' '==> Volume Fluide (m3) : ' VTOTAL ;
'TEMPS' 'ZERO' ;
*
* Ajout d'une matrice nulle pour préconditionner le système lors de
* l'utilisation de méthodes itératives
*
'SI' ('EGA' TYPI 3) ;
RR = 'EQEX'
'OPTI' 'EF' KSUPG 'IMPL' KPRES
'ZONE' (GEO.'$VTF') 'OPER' 'KMAB' 1.0 'INCO' 'UN' 'PRES' ;
RR.'INCO' = RV.'INCO' ;
A1 B1 = 'KMAB' RR.'1KMAB' ;
M1 = 'KCHT' GEO.'$VTF' 'VECT' 'SOMMET' 'COMP' '1UN' '2UN' (0. 0.) ;
RV.'ESPR' = 'KOPS' 'CMCT' B1 B1 M1 ;
'FINSI' ;
'SI' (interact) ;
'TRACER' (GEO.'VTF') 'TITR' 'Maillage...' ;
'FINSI' ;
*
* Boucle en temps
*
'REPETER' BCLTPS NPAS ;
NBT = 'DIME' (RV.'INCO'.'PM') ;
* temps courant
TCOUR = ('EXTRAIRE' (RV.'INCO'.'LTPS') NBT)
'+' (RV.'INCO'.'DT') ;
'MESSAGE' '==> Instant courant : ' TCOUR ;
RV.'PASDETPS'.'TPS' = TCOUR ;
RV.'PASDETPS'.'NUPASDT' = (RV.'PASDETPS'.'NUPASDT') '+' 1 ;
RT.'PASDETPS'.'NUPASDT' = (RT.'PASDETPS'.'NUPASDT') '+' 1 ;
RV.'INCO'.'LTPS' = (RV.'INCO'.'LTPS') 'ET' ('PROG' TCOUR) ;
* grandeurs 0D à l'instant précédent
PMN = 'EXTRAIRE' (RV.'INCO'.'PM') NBT ;
RHOMN = 'EXTRAIRE' (RV.'INCO'.'RHOM') NBT ;
REMN = 'EXTRAIRE' (RV.'INCO'.'REM') NBT ;
* Calcul de la masse 0D à l'instant courant
RHOM = 'MAXIMUM' ('PROG' (RHOMN '+' ((RV.'INCO'.'DT') '*'L '*'
MPH '/' VTOTAL)) 0.0 ) ;
PEXACT = P0 '+' (GAMA '*' MPH '*' RAIR '*' TH '*' L '*' TCOUR
'/' VTOTAL) ;
* Energie 0D
REM = REMN '+' ((RV.'INCO'.'DT') '*'
(MPH '*' L '*' CPAIR '*' TH)
'/' VTOTAL) ;
* boucle interne
ERIM = 'PROG' ;
NBCLI = 100 ;
'REPETER' BCLI NBCLI ;
* Calcul d'une pression garantissant la conservation de la masse pour
* cela on derive l'integrale de la loi d'état par rapport au temps
TFOLD = 'COPIER' (RV.'INCO'.'TF') ;
IRTF = 'INVERSE' (RAIR '*' (RV.'INCO'.'TF')) ;
IRTF = 'SOMT' (DIAG '*' IRTF) ;
PM = RHOM '*' VTOTAL '/' IRTF ;
'MESSAGE' '==> Pression 0D (Pa) : ' PM PEXACT ;
* Nouvelle condition aux limites pour la vitesse à l'injection
UINJ = MPH '*' RAIR '*' TH '/' PM ;
'MESSAGE' '==> Vitesse injection : ' UINJ ;
RV.'CLIM' = (RV.'CLIM') '/' V0 '*' UINJ ;
V0 = UINJ ;
* Correction du champ de vitesse transportant compte tenu de
* cette nouvelle vitesse: RV.'INCO'.'UN'
UNX = 'EXCO' (RV.'CLIM') '1UN' ;
UNX = 'NOMC' 'UX' UNX 'NATURE' 'DISCRET' ;
UNY = 'EXCO' (RV.'CLIM') '2UN' ;
UNY = 'NOMC' 'UY' UNY 'NATURE' 'DISCRET' ;
'MESSAGE' ('MAXIMUM' ('EXCO' (RV.'INCO'.'UN') 'UY')) ;
RV.'INCO'.'UN' = 'KCHT' (GEO.'$VTF') 'VECT' 'SOMMET'
(RV.'INCO'.'UN') (UNX 'ET' UNY) ;
'MESSAGE' ('MAXIMUM' ('EXCO' (RV.'INCO'.'UN') 'UY')) ;
* Calcul du terme source de l'équation d'énergie
'SI' ('EGA' RV.'PASDETPS'.'NUPASDT' 2) ;
DPDT = (PM '-' PMN) '/' (RV.'INCO'.'DT') ;
'SINON' ;
PMNN = 'EXTRAIRE' (RV.'INCO'.'PM') (NBT '-' 1) ;
DPDT = ((AN '*' PM) '+' (ANM1 '*' PMN) '+' (ANM2
'*' PMNN)) '/' (ADT '*' (RV.'INCO'.'DT')) ;
'FINSI' ;
'MESSAGE' '==> DPDT C.M. : ' DPDT ;
RV.'INCO'.'STF' = DPDT '/' CPAIR ;
* Calcul de la densité à partir de la loi d'état et du terme source
* de la divergence
RHO = ('INVERSE' ((RV.'INCO'.'TF') '*' RAIR ) ) '*' PM ;
'SI' ('EGA' RV.'PASDETPS'.'NUPASDT' 2) ;
DRDT = RHO '-' (RV.'INCO'.'RHONM')
'/' RV.'INCO'.'DT' ;
'SINON' ;
DRDT = ((RHO '*' AN) '+' ((RV.'INCO'.'RHONM') '*'
ANM1) '+' ((RV.'INCO'.'RHONMM') '*' ANM2))
'/' (ADT '*' (RV.'INCO'.'DT')) ;
'FINSI' ;
'MESSAGE' ('MAXIMUM' DRDT) ('MINIMUM' DRDT) ;
GRHO = 'KOPS' RHO (GEO.'$VTF') 'GRADS' ;
DRDT = DRDT '+' ( 'PSCA' (RV.'INCO'.'UN') ('MOTS' 'UX' 'UY')
GRHO ('MOTS' 'UX' 'UY')) ;
'MESSAGE' ('MAXIMUM' DRDT) ('MINIMUM' DRDT) ;
DRDT = DRDT '*' ('INVERSE' RHO) ;
'MESSAGE' ('MAXIMUM' DRDT) ('MINIMUM' DRDT) ;
* Compatibilité (intégrale sur le volume de Div u = 0
INTD = 'SOMT' (DRDT '*' DIAG) ;
* calcul de l'integrale de u.n ds sur le contour
TATA = 'DBIT' (RV.'INCO'.'UN') GEO.'$MENVF' ;
'MESSAGE' '==> Compatibilité Div u : ' INTD TATA;
RV.'INCO'.'DSRC' = DRDT '-' (INTD '/' VTOTAL) '+' (TATA '/'
VTOTAL) ;
'MESSAGE' ('SOMT' ((RV.'INCO'.'DSRC') '*' DIAG)) ;
*
RV.'INCO'.'RHO' = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' RHO ;
* résolution de la QDM et de la temperature
* Pour le cas BDF2, le premier pas doit être traité en BDF1
PAS1 = 'EGA' &BCLI 1 ;
'SI' (('EGA' KBDF 'BDF2') 'ET' PAS1) ;
'SI' ('EGA' RV.'PASDETPS'.'NUPASDT' 2) ;
RV.'4DFDT'.'KOPT'.'ISCHT' = 0 ;
RT.'4DFDT'.'KOPT'.'ISCHT' = 0 ;
'OUBLIER' (RT.'4DFDT') 'ARG3' ;
'OUBLIER' (RT.'4DFDT') 'ARG4' ;
'OUBLIER' (RV.'4DFDT') 'ARG3' ;
'OUBLIER' (RV.'4DFDT') 'ARG4' ;
RT.'4DFDT'.'ARG3' = 'DT' ;
RV.'4DFDT'.'ARG3' = 'DT' ;
RT.'4DFDT'.'IARG' = 3 ;
RV.'4DFDT'.'IARG' = 3 ;
'FINSI' ;
'SI' ('EGA' RV.'PASDETPS'.'NUPASDT' 3) ;
RV.'4DFDT'.'KOPT'.'ISCHT' = 1 ;
RT.'4DFDT'.'KOPT'.'ISCHT' = 1 ;
RT.'4DFDT'.'ARG4' = 'DT' ;
RV.'4DFDT'.'ARG4' = 'DT' ;
RT.'4DFDT'.'ARG3' = 'TFNMM' ;
RV.'4DFDT'.'ARG3' = 'UNMM' ;
RT.'4DFDT'.'IARG' = 4 ;
RV.'4DFDT'.'IARG' = 4 ;
'FINSI' ;
'FINSI' ;
RV.'METHINV'.'XINIT' = ('NOMC' '1UN' ('EXCO' (RV.'INCO'.'UN')
'UX') 'NATURE' 'DISCRET') '+'
('NOMC' '2UN' ('EXCO' (RV.'INCO'.'UN')
'UY') 'NATURE' 'DISCRET') '+'('NOMC' 'PRES' (RV.'INCO'.'PRES')
'NATURE' 'DISCRET') ;
'TEMPS' 'ZERO' ;
EXEC RV ;
'MESSAGE' 'Temps CPU = ' (('TEMPS' 'NOEC') '/' 100.0) ' s' ;
RT.'METHINV'.'XINIT' = 'NOMC' 'TF' (RV.'INCO'.'TF') 'NATURE'
'DISCRET' ;
'TEMPS' 'ZERO' ;
EXEC RT ;
'MESSAGE' 'Temps CPU = ' (('TEMPS' 'NOEC') '/' 100.0) ' s' ;
RES1 = 'MAXIMUM' ('ABS' ((RV.'INCO'.'TF') '-' TFOLD)) ;
'MESSAGE' '==> Residu RES : ' RES1 ;
'SI' ( interact) ;
VUN1 = RV.'INCO'.'UN' ;
VVN1 = 'VECT' VUN1 0.05 'UX' 'UY' 'JAUN' ;
'TRACE' RV.'INCO'.'TF' (GEO.'VTF') ('CONT' (GEO.'VTF'))
VVN1 'TITR' 'Champ de vitesse' ;
'FINSI' ;
* Calcul du terme source de la QDM
RHOP = (-1.) '*' (R0 '-' RHO) ;
* Ajout du terme de compressibilité dans Div Tau
TAU = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' ((RV.'INCO'.'DSRC')
'*' (RV.'INCO'.'MU') '*' (-1.0) '/' 3.0 ) ;
'MESSAGE' ('MAXIMUM' TAU) ('MINIMUM' TAU) ;
DIVT = 'KOPS' TAU GEO.'$VTF' 'GRADS' ;
DIVTX = ('EXCO' DIVT 'UX') ;
DIVTY = ('EXCO' DIVT 'UY') ;
'MESSAGE' ('MAXIMUM' DIVTX) ('MINIMUM' DIVTX) ;
'MESSAGE' ('MAXIMUM' DIVTY) ('MINIMUM' DIVTY) ;
RHOP = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' RHOP ;
ROGX = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET' DIVTX ;
ROGY = 'KCHT' (GEO.'$VTF') 'SCAL' 'SOMMET'
(DIVTY '-' (GRAV '*' RHOP));
ROGX = 'NOMC' 'UX' ROGX 'NATU' 'DISCRET' ;
ROGY = 'NOMC' 'UY' ROGY 'NATU' 'DISCRET' ;
RV.'INCO'.'ROG'= 'KCHT' (GEO.'$VTF') 'VECT' 'SOMMET'
(ROGX 'ET' ROGY) ;
* Affichage des bilans masse et énergie 0D/MultiD
MAS0D = RHOM '*' VTOTAL ;
ENR0D = REM '*' VTOTAL ;
MMD = 'SOMT' (DIAG '*' RHO) ;
EMD = 'SOMT' (DIAG '*' RHO '*' CVAIR '*' (RV.'INCO'.'TF')) ;
'MESSAGE' 'Bilan Masse 0D: ' MAS0D ' MultiD: ' MMD ;
'MESSAGE' 'Bilan Energie 0D: ' ENR0D ' MultiD: ' EMD ;
ERRORM = (MMD '-' MAS0D) '/' MAS0D '*' 100.0 ;
ERRORE = (EMD '-' ENR0D) '/' ENR0D '*' 100.0 ;
ERIM = ERIM 'ET' ('PROG' ERRORE) ;
'SI' ('<' RES1 1.E-4) ;
'QUITTER' BCLI ;
'FINSI' ;
'FIN' BCLI ;
'LISTE' ERIM ;
* sauvegarde et avancement en temps
RV.'INCO'.'RHOM' = RV.'INCO'.'RHOM' 'ET' ('PROG' RHOM) ;
RV.'INCO'.'REM' = RV.'INCO'.'REM' 'ET' ('PROG' REM) ;
RV.'INCO'.'PM' = RV.'INCO'.'PM' 'ET' ('PROG' PM) ;
RV.'INCO'.'PE' = RV.'INCO'.'PE' 'ET' ('PROG' PEXACT) ;
RV.'INCO'.'ERRORM' = RV.'INCO'.'ERRORM' 'ET'
('PROG' ERRORM) ;
RV.'INCO'.'ERRORE' = RV.'INCO'.'ERRORE' 'ET'
('PROG' ERRORE) ;
* mise a jour des champ multiD pour DFDT
RV.'INCO'.'UNMM' = 'COPIER' RV.'INCO'.'UNM' ;
RV.'INCO'.'TFNMM' = 'COPIER' RV.'INCO'.'TFNM' ;
RV.'INCO'.'RHONMM' = 'COPIER' RV.'INCO'.'RHONM' ;
RV.'INCO'.'UNM' = 'COPIER' RV.'INCO'.'UN' ;
RV.'INCO'.'RHONM' = 'COPIER' RV.'INCO'.'RHO' ;
RV.'INCO'.'TFNM' = 'COPIER' RV.'INCO'.'TF' ;
'SI' ('>EG' TCOUR TFIN) ;
'QUITTER' BCLTPS ;
'FINSI' ;
'FIN' BCLTPS ;
'TEMPS' ;
'TEMPS' 'PLACE' ;
*
* Creation d'une droite pour profil
*
mt = GEO.'VTF' ;
m1 = 'CHANGER' mt 'LIGNE' ;
pt1 = 'POIN' mt 'PROC' P2 ;
pt2 = 'POIN' mt 'PROC' P4 ;
p1 = 'POIN' mt 'DROIT' pt1 pt2 0.0001 ;
LMIL = 'ELEM' m1 'APPUYE' p1 ;
*
* Profils sur LMIL
* Température (K)
*==========================================================
* Test de non régression
*==========================================================
EVVT = 'EVOL' 'VERT' 'CHPO' (RV.'INCO'.'TF') LMIL ;
TREF = 'PROG' 600.00 606.09 612.13 618.67
625.08 632.00 638.60 645.44 651.59
657.31 661.76 664.78 665.90 664.56 660.85
653.95 644.61 631.93 617.26 599.82
581.27 561.06 540.84 520.36 500.90
482.52 465.88 451.28 438.73 428.64
420.45 414.55 410.11 407.34 405.46
404.50 403.94 403.74 403.65 403.63 403.62 ;
TCOUR = 'EXTRAIRE' EVVT 'ORDO' ;
DT = 'MAXIMUM' (TCOUR '-' TREF) 'ABS' ;
'MESSAGE' 'Erreur max sur le profil de température = ' DT ;
DTMAX = 0.01 ;
TEST = DT < DTMAX ;
'SI' TEST ;
'ERREUR' 0 ;
'SINON' ;
'ERREUR' 5 ;
'FINSI' ;
*
* End of dgibi file INJECTION AIR/AIR
*
*
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer