Sources Esope | Cast3M

flhu21
C FLHU21    SOURCE    CHAT      05/01/13    00:03:38     5004
      SUBROUTINE FLHU21(NESP,
     &                  GAMG,ROG,PG,UNG,UTG,UVG,
     &                  GAMD,ROD,PD,UND,UTD,UVD,
     &                  YG,YD,FLU1,FLU2,
     &                  CELLT,
     &                  LOGNC,MESERR,LOGAN)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FLHUS1
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 mèthode
C                      HUS de Coquel - Liou.
C
C                      Hybridation de van Leer "FVS"
C                      avec Osher "FDS" in  "reverse order decomposition."
C
C                      (voir:
C                       1)  Beccantini, Paillere,
C                           "Upwind Flux Splitting Schemes..."
C                           RAPPORT DMT 97//268
C                       2)  Coquel, Liou
C                           "Hybrid Upwind Splitting (HUS) by a Field-by-
C                           Field Decomposition "
C                           1993. NASA TM 106843)
C
C LANGUAGE          :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/TTMF
C
C************************************************************************
C
C APPELES
C
C FLHUS1 ------ FLUPVL
C    |
C     --------- FLUMVL
C    |
C     --------- RACHUS
C
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, i.e.
C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
C                         rho*un*y1, ...)
C
C     IFLU2           =  table de travaille, utilizé ici mais definie
C                        avant (pour gagner du temps  CPU).
C
C     CELLT           =  condition de stabilité, i.e.
C
C                        dT/diamax &lt; cellt
C
C     LOGNC, MESSER   =  paramètres qui segnalent des eventuels problèmes
C                        de convergence
C
C     LOGAN           =  .TRUE. -> anomalie detectée
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créé le 6.1.98
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
 
      IMPLICIT INTEGER(I-N)
      INTEGER NESP, I1, NMAX, NMAX0
      REAL*8  GAMG,ROG,PG
     &       ,GAMD,ROD,PD
     &       ,GM1G,AG2,AG,UNG,UTG,UVG,MG,HTG
     &       ,GM1D,AD2,AD,UND,UTD,UVD,MD,HTD
     &       ,UTG2,UTD2,UVG2,UVD2
     &       ,ZERO, ZERO0
     &       ,DU
     &       ,COEFFG, COEFFD, UTILG, UTILD
     &       ,USTAR,RSTARG,RSTARD,PSTAR,CELLG,CELLD
     &       ,ASTARG,ASTARD,MSTARG,MSTARD,USTAR2,HTSG,HTSD
     &       ,CELL
     &       ,AMG,AMD,CELLT0,CELLT
     &       ,YG(*),YD(*),FLU1(*),FLU2(*)
      CHARACTER*(40) MESERR
      LOGICAL LOGNC, LARCO, LOGI, LOGAN
C
C**** ZERO = tolerance d'egalite pour REAL*8
C     NMAX = numero max d'iterations in RACHUS
C            pour calculer l'etats d'intersection
C            des invariants de Riemann
C
C     LARCO = Logical; si .TRUE. on calcule le flux
C             relative aux fractiones massiques avec
C             la correction de Larrouturou
C             (voir J. Comp. Phys., 95, 1991)
C
      PARAMETER(NMAX=1000,ZERO = 1.0D-12,LARCO = .TRUE.)
C
      LOGNC  = .FALSE.
      MESERR = '                                        '
C
C**** YG, YD, FLU1, FLU2 déjà defini
C
C
C**** Calcul du flux avec le FVS de van Leer
C
C
C**** Onde de "gauche" a "droite"
C
      GM1G = GAMG - 1.0D0
      AG2 = GAMG * PG/ ROG
      AG = SQRT(AG2)
      UTG2 = UTG*UTG
      UVG2= UVG*UVG
      MG = UNG / AG
      HTG = AG2 /GM1G + 0.5D0 * (UNG*UNG + UTG2 + UVG2)
      AMG = ABS(MG)
      CELLT = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0)
     &        /((AMG+1.0D0)*AG)
C
      CALL FLUPV2(NESP,GAMG,ROG,MG,AG,UTG,UVG,HTG,YG,FLU1)
C
C
C**** Onde de "droite" a "gauche".
C
      GM1D = GAMD - 1.0D0
      AD2 = GAMD * PD/ ROD
      AD = SQRT(AD2)
      UTD2 = UTD*UTD
      UVD2 = UVD*UVD
      MD = UND / AD
      HTD = AD2 / GM1D +0.5D0 * (UND*UND + UTD2 + UVD2)
      AMD = ABS(MD)
      CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0)
     &        /((AMD+1.0D0)*AD)
      CELLT = MIN(CELLT,CELLT0)
C
      CALL FLUMV2(NESP,GAMD,ROD,MD,AD,UTD,UVD,HTD,YD,FLU2)
C
      DO I1 = 1, NESP+5
        FLU1(I1) = FLU1(I1) + FLU2(I1)
      ENDDO
C
C**** CALCUL des etats gauche et droite sur le champ linearment
C     dégénéré en imposant la 'reverse order decomposition'
C     d'Osher, i.e. la decomposition d'Osher 'physique'
C
C
C**** Osher reversal composition:
C
C     Vide -> Les invariants de Riemann ne se coupent pas
C          -> RETURN
C
      COEFFG  = 2.0D0 / GM1G
      COEFFD  = 2.0D0 / GM1D
      UTILG = UNG + COEFFG * AG
      UTILD = UND - COEFFD * AD
      DU = UTILG - UTILD
      IF(DU .GT. ZERO)THEN
C
C
C******* Calcul des intersections des Invariants de Riemann
C
C******* Protection des PARAMETER NMAX, ZERO
C         N.B. Les subroutines pouvent les changer
C
         NMAX0 = NMAX
         ZERO0 = ZERO
C
         CALL  RACHUS(NMAX0,ZERO0,
     &             GAMG,COEFFG,AG,UTILG,PG,
     &             GAMD,COEFFD,AD,UTILD,PD,
     &             USTAR,PSTAR,LOGNC,MESERR,LOGAN)
C
C******* Pas de convergence de la mèthode de Newton-Rapson en RACHUS
C
         IF(LOGNC) GOTO 9999
C
C******* Calcul des etats a Gauche et a Droite du champ linearment dégénéré
C
         CELLG = (PSTAR/PG)**(1.0D0/GAMG)
         CELLD = (PSTAR/PD)**(1.0D0/GAMD)
         ASTARG=AG*CELLG**(1.0D0/COEFFG)
         ASTARD=AD*CELLD**(1.0D0/COEFFD)
         RSTARG=ROG*CELLG
         RSTARD=ROD*CELLD
         MSTARG = USTAR / ASTARG
         MSTARD = USTAR / ASTARD
         USTAR2 = USTAR*USTAR
         HTSG = (ASTARG * ASTARG) / GM1G + 0.5D0*(USTAR2+UTG2+UVG2)
         HTSD = (ASTARD * ASTARD) / GM1D + 0.5D0*(USTAR2+UTD2+UVD2)
C
C******* CFL
C
         AMG = ABS(MSTARG)
         CELLT0 = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0)
     &        /((AMG+1.0D0)*ASTARG)
         CELLT = MIN(CELLT,CELLT0)
         AMD = ABS(MSTARD)
         CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0)
     &        /((AMD+1.0D0)*ASTARD)
         CELLT = MIN(CELLT,CELLT0)
C
C******* Partie antidiffusive du flux
C
         IF(USTAR .GT. 0.0D0)THEN
            CALL FLUMV2(NESP,
     &               GAMD,RSTARD,MSTARD,ASTARD,UTD,UVD,HTSD,YD,FLU2)
C
            DO I1 = 1, NESP+5
               FLU1(I1) = FLU1(I1) - FLU2(I1)
            ENDDO
C
            CALL FLUMV2(NESP,
     &               GAMG,RSTARG,MSTARG,ASTARG,UTG,UVG,HTSG,YG,FLU2)
C
            DO I1 = 1, NESP+5
               FLU1(I1) = FLU1(I1) + FLU2(I1)
            ENDDO
         ELSE
            CALL FLUPV2(NESP,
     &              GAMD,RSTARD,MSTARD,ASTARD,UTD,UVD,HTSD,YD,FLU2)
C
            DO I1 = 1, NESP+5
               FLU1(I1) = FLU1(I1) + FLU2(I1)
            ENDDO
C
            CALL FLUPV2(NESP,
     &               GAMG,RSTARG,MSTARG,ASTARG,UTG,UVG,HTSG,YG,FLU2)
C
            DO I1 = 1, NESP+5
               FLU1(I1) = FLU1(I1) - FLU2(I1)
            ENDDO
         ENDIF
C
C******* Correction de Larrouturou pour la TVD des fractions massiques.
C                                  pour la vitesse tangentielle
C
C        N.B. Pour le solver de Riemann exact cette propriete est
C        valide
C        Dans le cas de vide ce n'est pas necessaire:
C        VanLeer a la  TVD des fractions massiques.
C
         LOGI = LARCO
         IF(LOGI)THEN
            CELL = FLU1(1)
            IF(CELL .GT. 0.0D0)THEN
               FLU1(3) = CELL * UTG
               FLU1(4) = CELL * UVG
               DO I1 = 1, NESP
                  FLU1(I1+5) = CELL * YG(I1)
               ENDDO
            ELSE
               FLU1(3) = CELL * UTD
               FLU1(4) = CELL * UVD
               DO I1 = 1, NESP
                  FLU1(I1+5) = CELL * YD(I1)
               ENDDO
            ENDIF
         ENDIF
       ENDIF
C
 9999 CONTINUE
      RETURN
      END