* fichier : injection.dgibi ************************************************************************ * Section : Fluides Transitoire ************************************************************************ ************************************************************************ * 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' (' 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 ; TABTPS = TEMP 'NOEC'; tcpu = TABTPS.'TEMPS_CPU'.'INITIAL'; 'MESSAGE' 'Temps CPU = ' tcpu ' ms' ; RT.'METHINV'.'XINIT' = 'NOMC' 'TF' (RV.'INCO'.'TF') 'NATURE' 'DISCRET' ; 'TEMPS' 'ZERO' ; EXEC RT ; TABTPS = TEMP 'NOEC'; tcpu = TABTPS.'TEMPS_CPU'.'INITIAL'; 'MESSAGE' 'Temps CPU = ' tcpu ' ms' ; 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' ;