C FROEBM SOURCE KLOCZKO 05/06/14 21:15:08 5111 C FROEBM SOURCE BECC 02/06/18 21:15:55 4365 SUBROUTINE FROEBM(NESP, & GAMG,ROG,PG,UNG,UTG, & GAMD,ROD,PD,UND,UTD, & YG,YD,VINF,FLU, & CELLT) C C C PROJET : CASTEM 2000 C C NOM : FROEBM C C DESCRIPTION : Formulation Volumes Finis pour les Equations C d'Euler Multi-Espèces relatives à un mélange C de gaz parfaits. C C Calcul du flux aux interfaces avec la méthode C de Roe-Turkel (cf. Thèse Cécile Viozat) adaptée C aux écoulements à faible nombre de Mach. C C LANGUAGE : FORTRAN 77 C C AUTEUR : T. KLOCZKO DM2S/SFME/LTMF C C************************************************************************ C C APPELES : C C************************************************************************ C C**** Entrées: C C NESP = nombre d'especes considérées dans les Equations C d'Euler C C GAMG, GAMD = les "gamma" du gaz (gauche et droite) C C ROG, ROD = les densités C C PG, PD = les pressions C C UNG, UND = vitesses normales C C UTG, UTD = vitesses tangentielles C C VINF = vitesse de cut-off C C**** Sorties: C C FLU = table du flux a l'interface dans le repaire C (n,t), i.e. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht, C rho*un*y1, ...) C C CELLT = condition de stabilité, i.e. C C dT/diamax < cellt C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : 29.04.05 C C C C************************************************************************ C C N.B.: Toutes les variables sont DECLAREES C C C IMPLICIT NONE C C*** Déclaration des indices de boucles et du nombre d'espèces INTEGER NESP,I1,I2,I3 C C*** Déclaration des variables primitives REAL*8 ROG,PG,UNG,UTG,RETG,GAMG, & ROD,PD,UND,UTD,RETD,GAMD C C*** Déclaration des variables de cacul des incréments conservatifs REAL*8 XP(4) C C*** Déclaration des variables de l'état moyen de Roe REAL*8 SQROG,SQROD,ROM,ROS, & UNM,UTM, & HG,HD,HM, & CM,GAMM,GAM1 C C*** Déclaration de la condition de stabilité REAL*8 CELLT C C*** Délaration des variables de calcul du Mach de référence REAL*8 CG,CD,PHIM,VINF C C*** Déclaration des variables de calcul des valeurs propres REAL*8 VP(4),XU C C*** Déclaration des variables de calcul des matrices de passage REAL*8 R,S,TU,Q2, & MPI(4,4),MP(4,4) C C*** Déclaration des variables de calcul de la dissipation REAL*8 VDIS(4),DS(4) C C*** Déclaration des variables de calcul du flux numérique REAL*8 FC(4),FLU(*) C C*** Déclaration de la partie multiespèce REAL*8 YG(*),YD(*),FLUM C C C C*** DEBUT DES CALCULS *** C C C**** We include in the cut-off UNG,UND,UTG,UTD in order to C avoid low diffusivity in stagnation regions C VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0)) VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0)) C C**** Calcul des Deltas XP C RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG) & + (0.5D0*ROG*((UNG*UNG)+(UTG*UTG))) RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD) & + (0.5D0*ROD*((UND*UND)+(UTD*UTD))) C XP(1) = ROG - ROD XP(2) = (ROG * UNG) - (ROD * UND) XP(3) = (ROG * UTG) - (ROD * UTD) XP(4) = RETG - RETD C C*** Calcul de l'état moyen de Roe C GAMM = SQRT(GAMG) * SQRT(GAMD) GAM1 = GAMM - 1.D0 C SQROG = SQRT(ROG) SQROD = SQRT(ROD) ROM = SQROG * SQROD ROS = SQROG + SQROD C UNM = ((SQROG * UNG) + (SQROD * UND)) / ROS UTM = ((SQROG * UTG) + (SQROD * UTD)) / ROS C HG = (1.D0 / ROG) * (RETG + PG) HD = (1.D0 / ROD) * (RETD + PD) HM = ((SQROG * HG) + (SQROD * HD)) / ROS C CM = SQRT(GAM1 * (HM - 0.5D0 * (UNM**2 + UTM**2))) C C*** Calcul de la condition de stabilité C CELLT = 1.D0 / (CM + ABS(UNM)) C C*** Calcul du nombre de Mach de référence PHIM C CG = SQRT((GAMG * PG) / ROG) CD = SQRT((GAMD * PD) / ROD) IF(VINF.GE.MAX(CD,CG))THEN PHIM = 1.D0 ELSE PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD)) ENDIF C C*** Calcul des valeurs propres C XU = ((1.D0 - PHIM**2) * UNM)**2 + 4.D0 * (PHIM * CM)**2 C VP(1) = UNM VP(2) = UNM VP(3) = 0.5D0 * ((1.D0 + PHIM**2) * UNM + SQRT(XU)) VP(4) = 0.5D0 * ((1.D0 + PHIM**2) * UNM - SQRT(XU)) C C*** Calcul des matrices de passages C R = VP(3) - UNM * PHIM**2 S = VP(4) - UNM * PHIM**2 TU = 0.5D0 * (VP(4) - VP(3)) Q2 = 0.5D0 * (UNM**2 + UTM**2) C MPI(1,1) = 1.D0 - GAM1 * Q2 / (CM**2) MPI(1,2) = GAM1 * UNM / (CM**2) MPI(1,3) = GAM1 * UTM / (CM**2) MPI(1,4) = -GAM1 / (CM**2) MPI(2,1) = UTM MPI(2,2) = 0.D0 MPI(2,3) = -1.D0 MPI(2,4) = 0.D0 MPI(3,1) = (S * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU MPI(3,2) = -(S * UNM * GAM1 + (PHIM*CM)**2) / TU MPI(3,3) = -(S * UTM * GAM1) / TU MPI(3,4) = S * GAM1 /TU MPI(4,1) = -(R * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU MPI(4,2) = (R * UNM * GAM1 + (PHIM*CM)**2) / TU MPI(4,3) = (R * UTM * GAM1) / TU MPI(4,4) = -R * GAM1 / TU C MP(1,1) = 1.D0 MP(1,2) = 0.D0 MP(1,3) = 0.5D0 / (PHIM*CM)**2 MP(1,4) = 0.5D0 / (PHIM*CM)**2 MP(2,1) = UNM MP(2,2) = 0.D0 MP(2,3) = 0.5D0 * (UNM + R) / (PHIM*CM)**2 MP(2,4) = 0.5D0 * (UNM + S) / (PHIM*CM)**2 MP(3,1) = UTM MP(3,2) = -1.D0 MP(3,3) = 0.5D0 * UTM / (PHIM*CM)**2 MP(3,4) = 0.5D0 * UTM / (PHIM*CM)**2 MP(4,1) = Q2 MP(4,2) = -UTM MP(4,3) = 0.5D0 * (HM + R * UNM) / (PHIM*CM)**2 MP(4,4) = 0.5D0 * (HM + S * UNM) / (PHIM*CM)**2 C C*** Calcul de la dissipation C DO I1 = 1,4 VDIS(I1) = 0.5D0 * ABS(VP(I1)) ENDDO C DO I1 = 1,4 DS(I1) = 0.D0 DO I2 = 1,4 DO I3 = 1,4 DS(I1) = DS(I1) + MP(I1,I2) * VDIS(I2) * MPI(I2,I3) * XP(I3) ENDDO ENDDO ENDDO C C*** Calcul du flux convectif C FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND)) FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD)) FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD)) FC(4) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD))) C C*** Calcul du flux numérique C DO I1 = 1,4 FLU(I1) = FC(I1) + DS(I1) ENDDO C C*** Partie multiespèce C FLUM = FLU(1) IF(FLUM .GT. 0)THEN DO I1 = 1, NESP, 1 FLU(4+I1)=FLUM * YG(I1) ENDDO ELSE DO I1 = 1, NESP, 1 FLU(4+I1)=FLUM * YD(I1) ENDDO ENDIF C RETURN END