Esope source files

   1 : C FROEMO    SOURCE    CHAT      05/01/13    00:09:08     5004           
   2 :       SUBROUTINE FROEMO(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               :  FROEMO
  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 ROE (avec gamma = const)
  20 : C                      Voir: M. Delanaye, These de doctorat, Université
  21 : C                            de Liege
  22 : C                      N.B.: SPLIT DE FRACTIONS MASSIQUES A LA 
  23 : C                            LARROUTUROU
  24 : C
  25 : C LANGUAGE          :  FORTRAN 77 
  26 : C      
  27 : C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/LTMF  
  28 : C
  29 : C************************************************************************
  30 : C
  31 : C APPELES 
  32 : C
  33 : C************************************************************************
  34 : C
  35 : C**** Entrées:
  36 : C
  37 : C     NESP            =  nombre d'especes considérées dans les Equations
  38 : C                        d'Euler
  39 : C
  40 : C     GAMG, GAMD      =  les "gamma" du gaz (gauche et droite)
  41 : C
  42 : C     ROG, ROD        =  les densités 
  43 : C
  44 : C     PG, PD          =  les pressions
  45 : C
  46 : C     UNG, UND        =  vitesses normales
  47 : C
  48 : C     UTG, UTD        =  vitesses tangentielles
  49 : C
  50 : C     YG, YD          =  tables des fractiones massiques
  51 : C
  52 : C**** Sorties:
  53 : C
  54 : C     FLU1            =  table du flux a l'interface dans le repaire
  55 : C                        (n,t), i.e.
  56 : C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  57 : C                         rho*un*y1, ...)
  58 : C
  59 : C     CELLT           =  condition de stabilité, i.e.
  60 : C
  61 : C                        dT/diamax < cellt
  62 : C      
  63 : C************************************************************************
  64 : C
  65 : C HISTORIQUE (Anomalies et modifications éventuelles)
  66 : C
  67 : C HISTORIQUE :  Créé le 11.09.2000
  68 : C
  69 : C************************************************************************
  70 : C
  71 : C N.B.: Toutes les variables sont DECLAREES
  72 : C                          
  73 :       IMPLICIT INTEGER(I-N)
  74 :       INTEGER NESP, I1
  75 :       REAL*8  GAMG,ROG,PG,UNG,UTG
  76 :      &     ,GAMD,ROD,PD,UND,UTD
  77 :      &     ,YG(*),YD(*),FLU1(*),CELLT
  78 :      &     ,ONE,ZERO,HALF,EPSI,TWO
  79 :      &     ,GAMMA
  80 :      &     ,UNG2,UTG2,AG,AG2,HG,SQROG
  81 :      &     ,UND2,UTD2,AD,AD2,HD,SQROD
  82 :      &     ,DUN,DUT,DRO,DP,ROTILD,OMEGA,UNTILD,UTTILD,HTILD,ECTILD
  83 :      &     ,ATILD2,ATILD,LAM1,LAM2,LAM3
  84 :      &     ,LAMD1,LAMD2,LAMD3,LAMG1,LAMG2,LAMG3
  85 :      &     ,DELTA1,DELTA2,DELTA3,DD1,DD2,DD3
  86 :      &     ,DW1,DW2,DW3,DW4,ZZ1,ZZ2,ZZ3,ZZ4
  87 :      &     ,DF(4),FG(4),FD(4),CELL,FLU1G,FLU1D
  88 : C
  89 :       PARAMETER(ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,EPSI=1.0D-14,
  90 :      *     TWO=2.0D0)
  91 : C
  92 : C**** YG, YD, FLU1, déjà definis avant
  93 : C
  94 : C     Dans le cas Euler monoespece, on doit
  95 : C     avoir  :
  96 : C     YG(1) = YD(1) = 0.0D0
  97 : C
  98 :       GAMMA = HALF*(GAMG+GAMD)
  99 :       UNG2 = UNG*UNG
 100 :       UTG2 = UTG*UTG
 101 :       AG2 = GAMMA*PG/ROG
 102 :       AG = SQRT(AG2)
 103 :       HG = AG2/(GAMMA-ONE)+HALF*(UNG2+UTG2)
 104 :       SQROG = SQRT(ROG)
 105 : C
 106 :       UND2 = UND*UND
 107 :       UTD2 = UTD*UTD
 108 :       AD2 = GAMMA*PD/ROD
 109 :       AD = SQRT(AD2)
 110 :       HD = AD2/(GAMMA-ONE)+HALF*(UND2+UTD2)
 111 :       SQROD = SQRT(ROD)
 112 : C
 113 :       DUN = UND - UNG
 114 :       DUT = UTD - UTG
 115 :       DRO = ROD - ROG
 116 :       DP = PD - PG
 117 : C
 118 :       ROTILD = SQROG*SQROD
 119 :       OMEGA = SQROG/( SQROG + SQROD )
 120 :       UNTILD = UNG*OMEGA + UND*(ONE-OMEGA)
 121 :       UTTILD = UTG*OMEGA + UTD*(ONE-OMEGA)
 122 :       HTILD = HG*OMEGA + HD*(ONE-OMEGA)
 123 :       ECTILD = HALF* (UNTILD*UNTILD + UTTILD * UTTILD)
 124 :       ATILD2 = (GAMMA-ONE)*(HTILD-ECTILD)
 125 :       ATILD = SQRT(ATILD2)
 126 : C
 127 :       LAM1 = UNTILD + ATILD
 128 :       LAM2 = UNTILD - ATILD
 129 :       LAM3 = UNTILD
 130 : C
 131 : C**** entropy fix
 132 : C
 133 :       LAMD1 = UND + AD
 134 :       LAMD2 = UND - AD
 135 :       LAMD3 = UND
 136 :       LAMG1 = UNG + AG
 137 :       LAMG2 = UNG - AG
 138 :       LAMG3 = UNG
 139 : C
 140 :       DELTA1 = LAMD1 - LAMG1
 141 :       DELTA2 = LAMD2 - LAMG2
 142 :       DELTA3 = LAMD3 - LAMG3
 143 : C
 144 :       DD1 = 4.0D0*MAX(DELTA1,ZERO)+EPSI
 145 :       DD2 = 4.0D0*MAX(DELTA2,ZERO)+EPSI
 146 :       DD3 = 4.0D0*MAX(DELTA3,ZERO)+EPSI
 147 : C
 148 :       LAM1 = ABS(LAM1)
 149 :       LAM2 = ABS(LAM2)
 150 :       LAM3 = ABS(LAM3)
 151 : C
 152 :       IF(LAM1.LE.(HALF*DD1)) LAM1=LAM1*LAM1/DD1+DD1/4.0D0
 153 :       IF(LAM2.LE.(HALF*DD2)) LAM2=LAM2*LAM2/DD2+DD2/4.0D0
 154 :       IF(LAM3.LE.(HALF*DD3)) LAM3=LAM3*LAM3/DD3+DD3/4.0D0
 155 : C
 156 :       DW1 = (DRO - (DP / ATILD2))
 157 :       DW2 = -1.0D0 * DUT * ROTILD
 158 :       CELL =  DP / (ROTILD * ATILD)
 159 :       DW3 = (DUN + CELL) * ROTILD / (2 * ATILD)
 160 :       DW4 = (CELL - DUN) * ROTILD / (2 * ATILD)
 161 : C
 162 :       ZZ1 = HALF*DW1*LAM3
 163 :       ZZ2 = HALF*DW2*LAM3
 164 :       ZZ3 = HALF*DW3*LAM1
 165 :       ZZ4 = HALF*DW4*LAM2
 166 : C
 167 : C     [ZZ1] = [RO]
 168 : C     [ZZ2] = [UN]
 169 : C     [ZZ3] = [ZZ4] = [RO] 
 170 : C
 171 :       DF(1) = ZZ1 + ZZ3 + ZZ4
 172 :       DF(2) = ZZ1*UNTILD + ((ZZ3+ZZ4)*UNTILD) +((ZZ3-ZZ4)*ATILD)
 173 :       DF(3) = DF(1) * UTTILD - ZZ2
 174 :       DF(4) = (ZZ1 * ECTILD) - (ZZ2 * UTTILD) + ((ZZ3 + ZZ4) * HTILD) 
 175 :      &      + ((ZZ3 - ZZ4) * ATILD * UNTILD)
 176 : C
 177 :       FD(1)=ROD*UND
 178 :       FD(2)=ROD*UND2+PD
 179 :       FD(3)=ROD*UND*UTD
 180 :       FD(4)=ROD*UND*(AD2/(GAMMA-ONE)+(HALF*(UND2+UTD2)))
 181 : C
 182 :       FG(1)=ROG*UNG
 183 :       FG(2)=ROG*UNG2+PG
 184 :       FG(3)=ROG*UNG*UTG
 185 :       FG(4)=ROG*UNG*(AG2/(GAMMA-ONE)+(HALF*(UNG2+UTG2)))
 186 : C
 187 :       FLU1(1) = HALF*(FG(1)+FD(1)) -  DF(1)
 188 :       FLU1(2) = HALF*(FG(2)+FD(2)) -  DF(2)
 189 :       FLU1(3) = HALF*(FG(3)+FD(3)) -  DF(3)
 190 :       FLU1(4) = HALF*(FG(4)+FD(4)) -  DF(4)
 191 : C
 192 :       FLU1G = HALF * (FLU1(1) + ABS(FLU1(1)))
 193 :       FLU1D = HALF * (FLU1(1) - ABS(FLU1(1)))
 194 :       DO I1 = 1, NESP, 1
 195 :          FLU1(I1+4) = FLU1G * YG(I1) + FLU1D * YD(I1)
 196 :       ENDDO
 197 : C
 198 : C**** cellt
 199 : C
 200 :       CELLT=TWO*MAX(LAM1,LAM2,LAM3,
 201 :      &     ABS(LAMD1),ABS(LAMD2),ABS(LAMD3),
 202 :      &     ABS(LAMG1),ABS(LAMG2),ABS(LAMG3))
 203 :       CELLT=1.0D0/CELLT
 204 : C
 205 :       RETURN
 206 :       END
 207 : 
 208 : 
 209 :  
 210 :  
 211 :  
 212 :