Sources Esope | Cast3M

froebm
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 &lt; 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