Sources Esope | Cast3M

froemo
C FROEMO    SOURCE    CHAT      05/01/13    00:09:08     5004
      SUBROUTINE FROEMO(NESP,
     &     GAMG,ROG,PG,UNG,UTG,
     &     GAMD,ROD,PD,UND,UTD,
     &     YG,YD,FLU1,
     &     CELLT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FROEMO
C
C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations
C                      d'Euler Multi-Especes relatives à un melange
C                      de gaz ideals.
C
C                      Calcul du flux aux interfaces avec la methode
C                      de ROE (avec gamma = const)
C                      Voir: M. Delanaye, These de doctorat, Université
C                            de Liege
C                      N.B.: SPLIT DE FRACTIONS MASSIQUES A LA
C                            LARROUTUROU
C
C LANGUAGE          :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/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     YG, YD          =  tables des fractiones massiques
C
C**** Sorties:
C
C     FLU1            =  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 :  Créé le 11.09.2000
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
      IMPLICIT INTEGER(I-N)
      INTEGER NESP, I1
      REAL*8  GAMG,ROG,PG,UNG,UTG
     &     ,GAMD,ROD,PD,UND,UTD
     &     ,YG(*),YD(*),FLU1(*),CELLT
     &     ,ONE,ZERO,HALF,EPSI,TWO
     &     ,GAMMA
     &     ,UNG2,UTG2,AG,AG2,HG,SQROG
     &     ,UND2,UTD2,AD,AD2,HD,SQROD
     &     ,DUN,DUT,DRO,DP,ROTILD,OMEGA,UNTILD,UTTILD,HTILD,ECTILD
     &     ,ATILD2,ATILD,LAM1,LAM2,LAM3
     &     ,LAMD1,LAMD2,LAMD3,LAMG1,LAMG2,LAMG3
     &     ,DELTA1,DELTA2,DELTA3,DD1,DD2,DD3
     &     ,DW1,DW2,DW3,DW4,ZZ1,ZZ2,ZZ3,ZZ4
     &     ,DF(4),FG(4),FD(4),CELL,FLU1G,FLU1D
C
      PARAMETER(ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,EPSI=1.0D-14,
     *     TWO=2.0D0)
C
C**** YG, YD, FLU1, déjà definis avant
C
C     Dans le cas Euler monoespece, on doit
C     avoir  :
C     YG(1) = YD(1) = 0.0D0
C
      GAMMA = HALF*(GAMG+GAMD)
      UNG2 = UNG*UNG
      UTG2 = UTG*UTG
      AG2 = GAMMA*PG/ROG
      AG = SQRT(AG2)
      HG = AG2/(GAMMA-ONE)+HALF*(UNG2+UTG2)
      SQROG = SQRT(ROG)
C
      UND2 = UND*UND
      UTD2 = UTD*UTD
      AD2 = GAMMA*PD/ROD
      AD = SQRT(AD2)
      HD = AD2/(GAMMA-ONE)+HALF*(UND2+UTD2)
      SQROD = SQRT(ROD)
C
      DUN = UND - UNG
      DUT = UTD - UTG
      DRO = ROD - ROG
      DP = PD - PG
C
      ROTILD = SQROG*SQROD
      OMEGA = SQROG/( SQROG + SQROD )
      UNTILD = UNG*OMEGA + UND*(ONE-OMEGA)
      UTTILD = UTG*OMEGA + UTD*(ONE-OMEGA)
      HTILD = HG*OMEGA + HD*(ONE-OMEGA)
      ECTILD = HALF* (UNTILD*UNTILD + UTTILD * UTTILD)
      ATILD2 = (GAMMA-ONE)*(HTILD-ECTILD)
      ATILD = SQRT(ATILD2)
C
      LAM1 = UNTILD + ATILD
      LAM2 = UNTILD - ATILD
      LAM3 = UNTILD
C
C**** entropy fix
C
      LAMD1 = UND + AD
      LAMD2 = UND - AD
      LAMD3 = UND
      LAMG1 = UNG + AG
      LAMG2 = UNG - AG
      LAMG3 = UNG
C
      DELTA1 = LAMD1 - LAMG1
      DELTA2 = LAMD2 - LAMG2
      DELTA3 = LAMD3 - LAMG3
C
      DD1 = 4.0D0*MAX(DELTA1,ZERO)+EPSI
      DD2 = 4.0D0*MAX(DELTA2,ZERO)+EPSI
      DD3 = 4.0D0*MAX(DELTA3,ZERO)+EPSI
C
      LAM1 = ABS(LAM1)
      LAM2 = ABS(LAM2)
      LAM3 = ABS(LAM3)
C
      IF(LAM1.LE.(HALF*DD1)) LAM1=LAM1*LAM1/DD1+DD1/4.0D0
      IF(LAM2.LE.(HALF*DD2)) LAM2=LAM2*LAM2/DD2+DD2/4.0D0
      IF(LAM3.LE.(HALF*DD3)) LAM3=LAM3*LAM3/DD3+DD3/4.0D0
C
      DW1 = (DRO - (DP / ATILD2))
      DW2 = -1.0D0 * DUT * ROTILD
      CELL =  DP / (ROTILD * ATILD)
      DW3 = (DUN + CELL) * ROTILD / (2 * ATILD)
      DW4 = (CELL - DUN) * ROTILD / (2 * ATILD)
C
      ZZ1 = HALF*DW1*LAM3
      ZZ2 = HALF*DW2*LAM3
      ZZ3 = HALF*DW3*LAM1
      ZZ4 = HALF*DW4*LAM2
C
C     [ZZ1] = [RO]
C     [ZZ2] = [UN]
C     [ZZ3] = [ZZ4] = [RO]
C
      DF(1) = ZZ1 + ZZ3 + ZZ4
      DF(2) = ZZ1*UNTILD + ((ZZ3+ZZ4)*UNTILD) +((ZZ3-ZZ4)*ATILD)
      DF(3) = DF(1) * UTTILD - ZZ2
      DF(4) = (ZZ1 * ECTILD) - (ZZ2 * UTTILD) + ((ZZ3 + ZZ4) * HTILD)
     &      + ((ZZ3 - ZZ4) * ATILD * UNTILD)
C
      FD(1)=ROD*UND
      FD(2)=ROD*UND2+PD
      FD(3)=ROD*UND*UTD
      FD(4)=ROD*UND*(AD2/(GAMMA-ONE)+(HALF*(UND2+UTD2)))
C
      FG(1)=ROG*UNG
      FG(2)=ROG*UNG2+PG
      FG(3)=ROG*UNG*UTG
      FG(4)=ROG*UNG*(AG2/(GAMMA-ONE)+(HALF*(UNG2+UTG2)))
C
      FLU1(1) = HALF*(FG(1)+FD(1)) -  DF(1)
      FLU1(2) = HALF*(FG(2)+FD(2)) -  DF(2)
      FLU1(3) = HALF*(FG(3)+FD(3)) -  DF(3)
      FLU1(4) = HALF*(FG(4)+FD(4)) -  DF(4)
C
      FLU1G = HALF * (FLU1(1) + ABS(FLU1(1)))
      FLU1D = HALF * (FLU1(1) - ABS(FLU1(1)))
      DO I1 = 1, NESP, 1
         FLU1(I1+4) = FLU1G * YG(I1) + FLU1D * YD(I1)
      ENDDO
C
C**** cellt
C
      CELLT=TWO*MAX(LAM1,LAM2,LAM3,
     &     ABS(LAMD1),ABS(LAMD2),ABS(LAMD3),
     &     ABS(LAMG1),ABS(LAMG2),ABS(LAMG3))
      CELLT=1.0D0/CELLT
C
      RETURN
      END