Sources Esope | Cast3M

flurie
C FLURIE    SOURCE    CHAT      05/01/13    00:04:53     5004
      SUBROUTINE FLURIE(NESP,
     &                  GAMG,ROG,PG,UNG,UTG,
     &                  GAMD,ROD,PD,UND,UTD,
     &                  YG,YD,FLUX,
     &                  CELLT,
     &                  LOGNC,LOGAN,MESERR)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FLURIE
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                      de Godunov (Solveur de Riemann exact).
C
C                      Parametrisation de Smoller
C
C                      (voir:
C                       1) J. SMOLLER, "Shock Waves and Reaction
C                          Diffusion Equations", Springer Verlag, 1983;
C
C                       2) Beccantini, Paillere,
C                           "Upwind Flux Splitting Schemes..."
C                           RAPPORT DMT 97//268
C                       )
C
C LANGAGE           :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI  DRN/DMT/SEMT/TTMF
C
C************************************************************************
C
C APPELES
C
C FLURIE
C   |
C   |
C    ---- RIECOM ---- RACC ---- WNVXC ---- VLH1
C           |
C           |
C            -------- VLH1
C           |
C           |
C            -------- VLF1
C           |
C           |
C            -------- VLF3
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     FLUX            =  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     CELLT           =  condition de stabilité, i.e.
C
C                        dT/diamax &lt; cellt
C
C     LOGNC           =  si .TRUE. -> no convergence
C
C     LOGAN           =  si .TRUE. -> une anomalie a été detectée
C
C     MESERR          =  message d'erreur
C
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créé le 7.1.98
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
      IMPLICIT INTEGER(I-N)
      INTEGER NESP, NMAX, NMAX0,  I1
      REAL*8  ZERO, ZERO0, EPSI, EPSI0
     &      ,GAMG,ROG,PG,UNG,UTG
     &      ,GAMD,ROD,PD,UND,UTD
     &      ,VGRIL
     &      ,ROR,UNR,UTR,ROETR,PR
     &      ,VREL, FLUM, RHTR
     &      ,CELLT, AR, GAMR
     &      ,YG(*),YD(*),FLUX(*)
C
      CHARACTER*(40) MESERR
      LOGICAL LOGNC, LOGAN, LOGETD, LOGLAR
      PARAMETER(NMAX=1000,ZERO = 1.0D-14,EPSI=1D-12)
C
C**** ZERO = tolérance d'egalite pour REAL*8
C
C     NMAX = nombre max d'itérations in racc.eso
C            pour calculer l'etats d'intersection
C            invariants de Riemann - conditions de
C            Rankine-Hugoniot.
C
C     EPSI = erreur tolérée in  racc.eso
C
C
C**** Initialisation de LOGNC, LOGAN,MESERR ne doit pas etre faite ici,
C     mais avant, i.e.
C
C      LOGNC = .FALSE.
C      LOGAN = .FALSE.
C      MESERR(1:40) = '                                        '
C
C
C**** Les fractions massiques
C
C     YG, YD déjà definis
C
C**** Le flux à l'interfaces sont dans le repaire (n,t)
C
C
C
C**** VGRIL = pour ALE
C
      VGRIL = 0.0D0
C
C**** On protege le "PARAMETER"
C
      NMAX0 = NMAX
      ZERO0 = ZERO
      EPSI0 = EPSI
C
C**** Calcul de l'etat a l'interface
C
      CALL RIECOM(NMAX0,EPSI0,ZERO0,
     &            VGRIL,
     &            GAMG,ROG,PG,UNG,UTG,
     &            GAMD,ROD,PD,UND,UTD,
     &            ROR,UNR,UTR,ROETR,PR,GAMR,LOGETD,
     &            LOGAN,LOGNC,MESERR)
C
C**** Valeurs à l' interface (x/t=VGRIL)
C
C     ROR    = rho
C     UNR =  = un
C     UTR =  = ut
C     ROETR  = rho*et
C     PR     = p
C     LOGETD = .FALSE. -> on est à gauche de la discontinuité de contact.
C              .TRUE.  -> on est à droite de la discontinuité de contact.
C
      IF(LOGNC .OR. LOGAN) GOTO 9999
C
C**** Les flux
C
      VREL = UNR - VGRIL
      FLUM = ROR * VREL
      RHTR = ROETR + PR
      FLUX(1) = FLUM
      FLUX(2) = FLUM * UNR + PR
      FLUX(3) = FLUM * UTR
      FLUX(4) = VREL * RHTR
C
C**** TEST
C
C     On peut demontrer que le sign de la vitesse
C     sur la discontinuité de contact est le meme
C     que le sign de ROR*(UNR-VGRIL)
C     (voir Larrouturou,
C     J. Comp. Phys., 95, 1991)
C
C     Par contre, si ROR = 0, ceci n'est plus vrai,
C     mais dans ce cas RHO*U*Y = 0
C
C
C     FLUM &lt; 0 -> On est à droite
C
      LOGLAR = FLUM .LT. 0.0D0
      IF((ABS(FLUM) .GT. 1.0D-8) .AND. (LOGLAR .NEQV. LOGETD))  THEN
         MESERR = 'RIEMANN, flurie.eso (Larrouturou)       '
         LOGAN = .TRUE.
         GOTO 9999
      ENDIF
C
C**** Fin test
C
C
      IF(NESP .GT. 0)THEN
         IF(LOGETD)THEN
            DO I1 = 1, NESP
               FLUX(I1+4) = FLUM * YD(I1)
            ENDDO
         ELSE
            DO I1 = 1, NESP
               FLUX(I1+4) = FLUM * YG(I1)
            ENDDO
         ENDIF
      ENDIF
C
C**** CFL
C
C PPw AR = SQRT (GAMR * PR / ROR)
C PPw CELLT = 1.0D0 / (ABS(VREL)+AR)
C+PPw Dans l'espace, on n'entend pas crier...
      IF(PR.EQ.0.D0.AND.ROR.EQ.0.D0)THEN
        AR=0.D0
      ELSE
        AR = SQRT (GAMR * PR / ROR)
      ENDIF
      IF(AR.EQ.0.D0.AND.VREL.EQ.0.D0)THEN
        CELLT=1.D0/ZERO
      ELSE
        CELLT = 1.0D0 / (ABS(VREL)+AR)
      ENDIF
C+PPw
C
 9999 CONTINUE
C
      RETURN
      END