Esope source files

   1 : C FLHLLC    SOURCE    BECC      10/05/05    21:15:10     6674           
   2 :       SUBROUTINE FLHLLC(NESP,
   3 :      &     GAMG,ROG,PG,UNG,UTG,
   4 :      &     GAMD,ROD,PD,UND,UTD,
   5 :      &     YG,YD,FLU1,
   6 :      &     CELLT)
   7 : C
   8 : C************************************************************************
   9 : C
  10 : C PROJET            :  CASTEM 2000
  11 : C
  12 : C NOM               :  FLHLLC
  13 : C
  14 : C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations 
  15 : C                      d'Euler Multi-Especes relatives à un melange 
  16 : C                      de gaz ideals.
  17 : C
  18 : C                      Calcul du flux aux interfaces avec la methode
  19 : C                      de HLLC (avec gamma = const)
  20 : C                      Voir: * Batten et al. SIAM J. Sci. Comp.
  21 : C                              Vol. 18 pp 1553-6, 1997
  22 : C                            * Batten et al. JCP, 1997
  23 : C                              Vol. 137 pp 38--
  24 : C
  25 : C                      N.B.: SPLIT DE FRACTIONS MASSIQUES ET UNT A LA 
  26 : C                            LARROUTUROU
  27 : C
  28 : C LANGUAGE          :  FORTRAN 77 
  29 : C      
  30 : C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/LTMF  
  31 : C
  32 : C************************************************************************
  33 : C
  34 : C APPELES 
  35 : C
  36 : C************************************************************************
  37 : C
  38 : C**** Entrées:
  39 : C
  40 : C     NESP            =  nombre d'especes considérées dans les Equations
  41 : C                        d'Euler
  42 : C
  43 : C     GAMG, GAMD      =  les "gamma" du gaz (gauche et droite)
  44 : C
  45 : C     ROG, ROD        =  les densités 
  46 : C
  47 : C     PG, PD          =  les pressions
  48 : C
  49 : C     UNG, UND        =  vitesses normales
  50 : C
  51 : C     UTG, UTD        =  vitesses tangentielles
  52 : C
  53 : C     YG, YD          =  tables des fractiones massiques
  54 : C
  55 : C**** Sorties:
  56 : C
  57 : C     FLU1            =  table du flux a l'interface dans le repaire
  58 : C                        (n,t), i.e.
  59 : C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  60 : C                         rho*un*y1, ...)
  61 : C
  62 : C     CELLT           =  condition de stabilité, i.e.
  63 : C
  64 : C                        dT/diamax < cellt (temps/longueur)
  65 : C      
  66 : C************************************************************************
  67 : C
  68 : C HISTORIQUE (Anomalies et modifications éventuelles)
  69 : C
  70 : C HISTORIQUE :  Créé le 11.09.2000
  71 : C
  72 : C************************************************************************
  73 : C
  74 : C N.B.: Toutes les variables sont DECLAREES
  75 : C                          
  76 : C     It exactly captures:
  77 : C
  78 : C     Stationary shocks (test for instance
  79 : C
  80 : C       GAMG=1.4D0
  81 : C       GAMD=1.4D0
  82 : C       ROG=1.0D0
  83 : C       ROD=4.5584471D0
  84 : C       PG=1.0D5
  85 : C       PD=0.18279373D7
  86 : C       UNG=0.14877919D4
  87 : C       UND=0.32638130D3
  88 : C       UTG=0.0D0
  89 : C       UTD=0.0D0
  90 : C     ) and stationary contacts
  91 : C
  92 :       IMPLICIT INTEGER(I-N)
  93 :       INTEGER NESP
  94 :      &     , IESP
  95 :       REAL*8  GAMG,ROG,PG,UNG,UTG
  96 :      &     ,GAMD,ROD,PD,UND,UTD
  97 :      &     ,YG(*),YD(*),FLU1(*),CELLT
  98 :      &     ,HALF
  99 :       PARAMETER (HALF=0.5D0)
 100 :       REAL*8 USGGM1,USGDM1,HTG,HTD,REING,REIND,SQRG,SQRD,USOSQR
 101 :      &     ,HTAVE,UNAVE,REIAVE,USGAVE,GAMAVE,ROAVE,PAVE,ASONAV
 102 :      &     ,ASONG,ASOND,SG,SD,SM,PSTAR,RETG,RETD,OMEGA
 103 :      &     ,ROGS,RUNGS,RETGS,RODS,RUNDS,RETDS
 104 : C
 105 :       USGGM1=1.0D0/(GAMG-1.0D0)
 106 :       USGDM1=1.0D0/(GAMD-1.0D0)
 107 : C
 108 : C     Evaluation of the Roe average state to compute the non-linear wave
 109 : C     speeds (According to Shuye, JCP 142, pag. 217--, 1998
 110 : C
 111 :       REING=PG*USGGM1
 112 :       REIND=PD*USGDM1
 113 :       HTG=(HALF*UNG*UNG)+(GAMG*REING/ROG)
 114 :       HTD=(HALF*UND*UND)+(GAMD*REIND/ROD)
 115 :       SQRG=ROG**0.5D0
 116 :       SQRD=ROD**0.5D0
 117 :       USOSQR=1.0D0/(SQRG+SQRD)
 118 :       HTAVE=((SQRG*HTG)+(SQRD*HTD))*USOSQR
 119 :       UNAVE=((SQRG*UNG)+(SQRD*UND))*USOSQR
 120 :       USGAVE=((SQRG*USGGM1)+(SQRD*USGDM1))*USOSQR
 121 :       GAMAVE=(1.0D0/USGAVE)+1.0D0
 122 :       REIAVE=((SQRG*REING)+(SQRD*REIND))*USOSQR
 123 : C    
 124 :       ROAVE=REIAVE*GAMAVE/(HTAVE-(HALF*UNAVE*UNAVE))
 125 :       PAVE=REIAVE*(GAMAVE-1.0D0)
 126 : C
 127 :       ASONAV=(GAMAVE*PAVE/ROAVE)**0.5D0
 128 : C
 129 : C**** Evaluation of the speed of the Genuinely-Non-Linear Waves
 130 : C
 131 :       ASONG=(GAMG*PG/ROG)**0.5D0
 132 :       ASOND=(GAMD*PD/ROD)**0.5D0
 133 : C
 134 :       SG=UNG-ASONG
 135 :       SG=MIN(SG,(UNAVE-ASONAV))
 136 :       SD=UND+ASOND
 137 :       SD=MAX(SD,(UNAVE+ASONAV))
 138 : C
 139 : C**** Speed of the contact discontinuity
 140 : C
 141 :       SM=ROD*UND*(SD-UND)-ROG*UNG*(SG-UNG)+PG-PD
 142 :       SM=SM/((ROD*(SD-UND))-(ROG*(SG-UNG)))
 143 : C
 144 : C**** Interfacial pressure
 145 : C
 146 :       PSTAR=ROG*(UNG-SG)*(UNG-SM)+PG
 147 :       RETG=(USGGM1*PG)+(HALF*ROG*((UNG*UNG)+(UTG*UTG)))
 148 :       RETD=(USGDM1*PD)+(HALF*ROD*((UND*UND)+(UTD*UTD)))
 149 : C
 150 :       IF(SG .GT. 0)THEN
 151 :          FLU1(1)=ROG*UNG
 152 :          FLU1(2)=ROG*UNG*UNG+PG
 153 :          FLU1(4)=UNG*(RETG+PG)
 154 :       ELSEIF(SM .GT. 0)THEN
 155 :          OMEGA=1.0D0/(SG-SM)
 156 :          ROGS=ROG*(SG-UNG)
 157 :          ROGS=ROGS*OMEGA
 158 :          RUNGS=((SG-UNG)*(ROG*UNG))+(PSTAR - PG)
 159 :          RUNGS=RUNGS*OMEGA
 160 :          RETGS=((SG-UNG)*RETG)+(PSTAR*SM)-(PG*UNG)
 161 :          RETGS=RETGS*OMEGA
 162 : C
 163 :          FLU1(1)=ROGS*SM
 164 :          FLU1(2)=RUNGS*SM+PSTAR
 165 :          FLU1(4)=(RETGS+PSTAR)*SM
 166 : C
 167 :       ELSEIF(SD .GT. 0)THEN
 168 :          OMEGA=1.0D0/(SD-SM)
 169 :          RODS=ROD*(SD-UND)
 170 :          RODS=RODS*OMEGA
 171 :          RUNDS=((SD-UND)*(ROD*UND))+(PSTAR - PD)
 172 :          RUNDS=RUNDS*OMEGA
 173 :          RETDS=((SD-UND)*RETD)+(PSTAR*SM)-(PD*UND)
 174 :          RETDS=RETDS*OMEGA
 175 : C
 176 :          FLU1(1)=RODS*SM
 177 :          FLU1(2)=RUNDS*SM+PSTAR
 178 :          FLU1(4)=(RETDS+PSTAR)*SM
 179 :       ELSE
 180 :          FLU1(1)=ROD*UND
 181 :          FLU1(2)=ROD*UND*UND+PD
 182 :          FLU1(4)=UND*(RETD+PD)
 183 :       ENDIF
 184 : C
 185 :       FLU1(3)=(FLU1(1)*(UTG+UTD))+(ABS(FLU1(1))*(UTG-UTD))
 186 :       FLU1(3)=HALF*FLU1(3)
 187 : C
 188 :       DO IESP=1,NESP,1
 189 :          FLU1(4+IESP)=(FLU1(1)*(YG(IESP)+YD(IESP)))+
 190 :      &        (ABS(FLU1(1))*(YG(IESP)-YD(IESP)))
 191 :          FLU1(4+IESP)=HALF*FLU1(4+IESP)
 192 :       ENDDO
 193 : C
 194 :       CELLT=1.0D0/(MAX(SG,SD))
 195 : C
 196 :       RETURN
 197 :       END
 198 : 
 199 : 
 200 :  
 201 :  
 202 :  
 203 :  
 204 :