Esope source files

   1 : C FDV22F    SOURCE    BECC      05/12/02    21:15:12     5258           
   2 :       SUBROUTINE fDV22F(INDMET,
   3 :      &                  alphaL, uvnL, uvtL, ulnL, ultL, pL, TvL, TlL, 
   4 :      &                  rvL, rlL,
   5 :      &                  alphaR, uvnR, uvtR, ulnR, ultR, pR, TvR, TlR, 
   6 :      &                  rvR, rlR,
   7 :      &                  F)
   8 : C***********************************************************************
   9 : C NOM         : fADV22F 
  10 : C DESCRIPTION : Calculation of the Two phase AUSMDV flux
  11 : C               (2D two fluid model)
  12 : C
  13 : C LANGAGE     : ESOPE
  14 : C AUTEUR      : José R. Garcia Cascales (CEA/DEN/DM2S/SFME/LTMF)
  15 : C               mél : fd1@semt2.smts.cea.fr
  16 : C***********************************************************************
  17 : C APPELES          : 
  18 : C APPELES (E/S)    : 
  19 : C     inputs       : INDMET, variable to say if we want 
  20 : C                    preconditioned version or not
  21 : C                    primitive variable at left side
  22 : C                    AL, UVNL, UVTL, ULNL, ULTL, PL, TVL, TLL, 
  23 : C                    RVG, RLG,
  24 : C                    primitive variables at right side
  25 : C                    AR, UVNR, UVTR, ULNR, ULTR, PR, TVR, TLR, 
  26 : C                    RVR, RLR,
  27 : C
  28 : C     output       : F, 2D AUSMDV or "preconditioned" AUSMDV flux
  29 : C***********************************************************************
  30 : C SYNTAXE GIBIANE    :
  31 : C ENTREES            :
  32 : C ENTREES/SORTIES    :
  33 : C SORTIES            :
  34 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  35 : C***********************************************************************
  36 : C VERSION    : v1, 13/03/2002, version initiale
  37 : C***********************************************************************
  38 : *
  39 :       real*8 F(8)
  40 :       real*8 alphaL, alphaR, uvnL, uvnR, uvtL, uvtR,
  41 :      &       ulnL, ulnR, ultL, ultR, pL, pR, evL, evR, elL,
  42 :      &       elR, hvL, hvR, hlL, hlR, TvL, TvR, TlL, TlR
  43 :       integer INDMET
  44 : C variables for AUSMDV scheme
  45 :       real*8 rvL, rvR, rlL, rlR, avL, avR, alL, alR
  46 :       real*8 alm, avm, MlL, MlR, MvL, MvR, avpm, alpm, 
  47 :      &       dotmv, dotml, Mvplus, Mvmin, Mlplus, Mlmin, 
  48 :      &       Pvplus, Pvmin, Plplus, Plmin,
  49 :      &       fvL, fvR, flL, flR, wvplus, wvmin, wlplus, wlmin
  50 : 
  51 :       real*8 gamv, cpv, gaml, cpl, pil
  52 :       PARAMETER(gamv = 1.4D0, cpv = 1008.D0, 
  53 :      &          gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)
  54 : 
  55 :       evL = cpv*TvL/gamv
  56 :       evR = cpv*TvR/gamv
  57 :       elL = cpl*TlL/gaml + pil/rlL
  58 :       elR = cpl*TlR/gaml + pil/rlR
  59 : 
  60 :       hvL = evL + pL/rvL
  61 :       hvR = evR + pR/rvR
  62 :       hlL = elL + pL/rlL
  63 :       hlR = elR + pR/rlR
  64 : 
  65 : C Speed of sound at left and rigth sides
  66 :       avL = dsqrt((gamv - 1.D0)*Cpv*TvL)
  67 :       avR = dsqrt((gamv - 1.D0)*Cpv*TvR)
  68 :       avm = dsqrt(avL*avR)
  69 :       alL = dsqrt((gaml - 1.D0)*Cpl*TlL)
  70 :       alR = dsqrt((gaml - 1.D0)*Cpl*TlR)
  71 :       alm = dsqrt(alL*alR)
  72 : 
  73 : C Match number at both sides
  74 :       MlL = ulnL/alm
  75 :       MlR = ulnR/alm
  76 :       MvL = uvnL/avm
  77 :       MvR = uvnR/avm
  78 : 
  79 : C AUSMDV scheme
  80 : C      fvL = pL/(alphaL*rvL)
  81 : C      fvR = pR/(alphaR*rvR)
  82 : C      flL = pL/((1.D0 - alphaL)*rlL)
  83 : C      flR = pR/((1.D0 - alphaR)*rlR)
  84 : 
  85 :       fvL = 1.D0/(alphaL*rvL)
  86 :       fvR = 1.D0/(alphaR*rvR)
  87 :       flL = 1.D0/((1.D0 - alphaL)*rlL)
  88 :       flR = 1.D0/((1.D0 - alphaR)*rlR)
  89 : 
  90 : C     fvL = pL/(rvL)
  91 : C     fvR = pR/(rvR)
  92 : C     flL = pL/(rlL)
  93 : C     flR = pR/(rlR)
  94 : 
  95 : C     fvL = 1.D0/(rvL)
  96 : C     fvR = 1.D0/(rvR)
  97 : C     flL = 1.D0/(rlL)
  98 : C     flR = 1.D0/(rlR)
  99 : 
 100 :       wvplus = 2.D0*fvL/(fvL + fvR)
 101 :       wvmin = 2.D0*fvR/(fvL + fvR)
 102 :       wlplus = 2.D0*flL/(flL + flR)
 103 :       wlmin = 2.D0*flR/(flL + flR)
 104 : 
 104 :       call mpADV(wvplus, wvmin, MvL, MvR, Mvplus, Mvmin, Pvplus, Pvmin)
 104 :       call mpADV(wlplus, wlmin, MlL, MlR, Mlplus, Mlmin, Plplus, Plmin)
 105 : 
 106 : C AUSM schemes mass fluxes
 107 :       dotmv = avm*(alphaL*rvL*Mvplus + alphaR*rvR*Mvmin)
 108 :       dotml = alm*((1.D0 - alphaL)*rlL*Mlplus + (1.D0 - alphaR)*rlR
 109 :      $     *Mlmin)
 110 :       if(INDMET .EQ. 4)then
 111 :          dotml = dotml + 0.01d0*(pL - pR)
 112 :       end if
 113 : C AUSM schemes pressure fluxes
 114 :       avpm = Pvplus*alphaL*pL + Pvmin*alphaR*pR
 115 :       alpm = Plplus*(1.D0 - alphaL)*pL + Plmin*(1.D0 - alphaR)
 116 :      $     *pR
 117 : 
 118 : C Definition of convective part of the flux
 119 :       if (dotmv .GT. 0.D0) then
 120 : 
 121 :          F(1) = dotmv*1.D0
 122 :          F(3) = dotmv*uvnL
 123 :          F(5) = dotmv*uvtL
 124 :          F(7) = dotmv*(hvL + (uvnL**2 + uvtL**2)/2.D0)
 125 :       else
 126 :          F(1) = dotmv*1.D0
 127 :          F(3) = dotmv*uvnR
 128 :          F(5) = dotmv*uvtR
 129 :          F(7) = dotmv*(hvR + (uvnR**2 + uvtR**2)/2.D0)
 130 :       end if
 131 :       
 132 :       if (dotml .GT. 0.D0) then
 133 :          F(2) = dotml*1.D0
 134 :          F(4) = dotml*ulnL
 135 :          F(6) = dotml*ultL
 136 :          F(8) = dotml*(hlL + (ulnL**2 + ultL**2)/2.D0)
 137 :       else
 138 :          F(2) = dotml*1.D0
 139 :          F(4) = dotml*ulnR
 140 :          F(6) = dotml*ultR
 141 :          F(8) = dotml*(hlR + (ulnR**2 + ultR**2)/2.D0)
 142 :       end if
 143 : 
 144 : C Addition of the pressure term
 145 : 
 146 :       F(3) = F(3) + avpm
 147 :       F(4) = F(4) + alpm
 148 : 
 149 :       END
 150 : 
 151 :