Esope source files

   1 : C FLURI2    SOURCE    CHAT      05/01/13    00:04:50     5004           
   2 : C FLURIE    SOURCE    BECC      97/12/23    22:14:53     3017           
   3 :       SUBROUTINE FLURI2(NESP,
   4 :      &                  GAMG,ROG,PG,UNG,UTG,UVG,
   5 :      &                  GAMD,ROD,PD,UND,UTD,UVD,
   6 :      &                  YG,YD,FLUX,
   7 :      &                  CELLT,
   8 :      &                  LOGNC,LOGAN,MESERR)
   9 : C
  10 : C************************************************************************
  11 : C
  12 : C PROJET            :  CASTEM 2000
  13 : C
  14 : C NOM               :  FLURIE
  15 : C
  16 : C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations 
  17 : C                      d'Euler Multi-Especes relatives à un melange 
  18 : C                      de gaz ideals.
  19 : C
  20 : C                      Calcul du flux aux interfaces avec la mèthode
  21 : C                      de Godunov (Solveur de Riemann exact).
  22 : C
  23 : C                      Parametrisation de Smoller
  24 : C      
  25 : C                      (voir:
  26 : C                       1) J. SMOLLER, "Shock Waves and Reaction 
  27 : C                          Diffusion Equations", Springer Verlag, 1983;
  28 : C
  29 : C                       2) Beccantini, Paillere,
  30 : C                           "Upwind Flux Splitting Schemes..."
  31 : C                           RAPPORT DMT 97//268
  32 : C                       )      
  33 : C
  34 : C LANGAGE           :  FORTRAN 77 
  35 : C      
  36 : C AUTEUR            :  A. BECCANTINI  DRN/DMT/SEMT/TTMF  
  37 : C
  38 : C************************************************************************
  39 : C
  40 : C APPELES
  41 : C
  42 : C FLURIE
  43 : C   |   
  44 : C   |
  45 : C    ---- RIECOM ---- RACC ---- WNVXC ---- VLH1
  46 : C           |
  47 : C           |
  48 : C            -------- VLH1
  49 : C           |
  50 : C           |
  51 : C            -------- VLF1
  52 : C           |
  53 : C           |      
  54 : C            -------- VLF3
  55 : C      
  56 : C************************************************************************
  57 : C
  58 : C**** Entrées:
  59 : C
  60 : C     NESP            =  nombre d'especes considérées dans les Equations
  61 : C                        d'Euler
  62 : C
  63 : C     GAMG, GAMD      =  les "gamma" du gaz (gauche et droite)
  64 : C
  65 : C     ROG, ROD        =  les densités 
  66 : C
  67 : C     PG, PD          =  les pressions
  68 : C
  69 : C     UNG, UND        =  vitesses normales
  70 : C
  71 : C     UTG, UTD        =  vitesses tangentielles
  72 : C
  73 : C     UVG, UVD        =  vitesses tangentielles
  74 : C
  75 : C     YG, YD          =  tables des fractiones massiques
  76 : C
  77 : C**** Sorties:
  78 : C
  79 : C     FLUX            =  table du flux a l'interface, i.e.
  80 : C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  81 : C                         rho*un*y1, ...)
  82 : C
  83 : C     CELLT           =  condition de stabilité, i.e.
  84 : C
  85 : C                        dT/diamax < cellt
  86 : C
  87 : C     LOGNC           =  si .TRUE. -> no convergence
  88 : C
  89 : C     LOGAN           =  si .TRUE. -> une anomalie a été detectée
  90 : C
  91 : C     MESERR          =  message d'erreur
  92 : C
  93 : C      
  94 : C************************************************************************
  95 : C
  96 : C HISTORIQUE (Anomalies et modifications éventuelles)
  97 : C
  98 : C HISTORIQUE :  Créé le 7.1.98 
  99 : C
 100 : C************************************************************************
 101 : C
 102 : C N.B.: Toutes les variables sont DECLAREES
 103 : C                          
 104 :       IMPLICIT INTEGER(I-N)
 105 :       INTEGER NESP, NMAX, NMAX0,  I1
 106 :       REAL*8  ZERO, ZERO0, EPSI, EPSI0
 107 :      &      ,GAMG,ROG,PG,UNG,UTG,UVG
 108 :      &      ,GAMD,ROD,PD,UND,UTD,UVD
 109 :      &      ,VGRIL
 110 :      &      ,ROR,UNR,UTR,UVR,ROETR,PR
 111 :      &      ,VREL, FLUM, RHTR
 112 :      &      ,CELLT, AR, GAMR
 113 :      &      ,YG(*),YD(*),FLUX(*)
 114 : C
 115 :       CHARACTER*(40) MESERR
 116 :       LOGICAL LOGNC, LOGAN, LOGETD, LOGLAR
 117 :       PARAMETER(NMAX=1000,ZERO = 1.0D-14,EPSI=1D-12)
 118 : C
 119 : C**** ZERO = tolérance d'egalite pour REAL*8
 120 : C      
 121 : C     NMAX = nombre max d'itérations in racc.eso
 122 : C            pour calculer l'etats d'intersection
 123 : C            invariants de Riemann - conditions de
 124 : C            Rankine-Hugoniot. 
 125 : C                     
 126 : C     EPSI = erreur tolérée in  racc.eso
 127 : C
 128 : C
 129 : C**** Initialisation de LOGNC, LOGAN,MESERR ne doit pas etre faite ici,
 130 : C     mais avant, i.e.
 131 : C
 132 : C      LOGNC = .FALSE.
 133 : C      LOGAN = .FALSE.
 134 : C      MESERR(1:40) = '                                        '
 135 : C      
 136 : C
 137 : C**** Les fractions massiques
 138 : C
 139 : C     YG, YD déjà definis 
 140 : C      
 141 : C**** Le flux à l'interfaces sont dans le repaire (n,t)
 142 : C
 143 : C      
 144 : C
 145 : C**** VGRIL = pour ALE
 146 : C
 147 :       VGRIL = 0.0D0
 148 : C
 149 : C**** On protege le "PARAMETER"
 150 : C      
 151 :       NMAX0 = NMAX
 152 :       ZERO0 = ZERO
 153 :       EPSI0 = EPSI
 154 : C
 155 : C**** Calcul de l'etat a l'interface
 156 : C
 156 :       CALL RIECO2(NMAX0,EPSI0,ZERO0,
 157 :      &            VGRIL,
 158 :      &            GAMG,ROG,PG,UNG,UTG,UVG,
 159 :      &            GAMD,ROD,PD,UND,UTD,UVD,
 160 :      &            ROR,UNR,UTR,UVR,ROETR,PR,GAMR,LOGETD,
 161 :      &            LOGAN,LOGNC,MESERR)
 162 : C
 163 : C**** Valeurs à l' interface (x/t=VGRIL)
 164 : C
 165 : C     ROR    = rho
 166 : C     UNR =  = un
 167 : C     UTR =  = ut      
 168 : C     ROETR  = rho*et
 169 : C     PR     = p      
 170 : C     LOGETD = .FALSE. -> on est à gauche de la discontinuité de contact.
 171 : C              .TRUE.  -> on est à droite de la discontinuité de contact. 
 172 : C
 173 :       IF(LOGNC .OR. LOGAN) GOTO 9999
 174 : C
 175 : C**** Les flux      
 176 : C
 177 :       VREL = UNR - VGRIL
 178 :       FLUM = ROR * VREL 
 179 :       RHTR = ROETR + PR
 180 :       FLUX(1) = FLUM
 181 :       FLUX(2) = FLUM * UNR + PR
 182 :       FLUX(3) = FLUM * UTR
 183 :       FLUX(4) = FLUM * UVR
 184 :       FLUX(5) = VREL * RHTR
 185 : C
 186 : C**** TEST
 187 : C
 188 : C     On peut demontrer que le sign de la vitesse
 189 : C     sur la discontinuité de contact est le meme
 190 : C     que le sign de ROR*(UNR-VGRIL) 
 191 : C     (voir Larrouturou,
 192 : C     J. Comp. Phys., 95, 1991)
 193 : C
 194 : C     Par contre, si ROR = 0, ceci n'est plus vrai,
 195 : C     mais dans ce cas RHO*U*Y = 0
 196 : C
 197 : C
 198 : C     FLUM < 0 -> On est à droite
 199 : C
 200 :       LOGLAR = FLUM .LT. 0.0D0      
 201 :       IF((ABS(FLUM) .GT. 1.0D-8) .AND. (LOGLAR .NEQV. LOGETD))  THEN
 202 :          MESERR = 'RIEMANN, flurie.eso (Larrouturou)       '
 203 :          LOGAN = .TRUE.
 204 :          GOTO 9999
 205 :       ENDIF
 206 : C
 207 : C**** Fin test
 208 : C
 209 : C           
 210 :       IF(NESP .GT. 0)THEN
 211 :          IF(LOGETD)THEN
 212 :             DO I1 = 1, NESP
 213 :                FLUX(I1+5) = FLUM * YD(I1)
 214 :             ENDDO
 215 :          ELSE
 216 :             DO I1 = 1, NESP
 217 :                FLUX(I1+5) = FLUM * YG(I1)
 218 :             ENDDO
 219 :          ENDIF
 220 :       ENDIF
 221 : C
 222 : C**** CFL
 223 : C
 224 :       AR = SQRT (GAMR * PR / ROR)
 225 :       CELLT = 1.0D0 / (ABS(VREL)+AR)
 226 : C     
 227 :  9999 CONTINUE
 228 : C      
 229 :       RETURN
 230 :       END
 231 : 
 232 : 
 233 :  
 234 :