Sources Esope | Cast3M

fausm2
C FAUSM2    SOURCE    CHAT      05/01/12    23:56:44     5004
       SUBROUTINE FAUSM2(NESP,
     &                  GAMG,ROG,PG,UNG,UTG,
     &                  GAMD,ROD,PD,UND,UTD,
     &                  YG,YD,V_INF,F,
     &                  CELLT)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FAUSM2 ('FLUX for Bas Mach')
C
C DESCRIPTION       :  Voir KONJA2
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
C
C************************************************************************
C
c----------------------------------------------------------------------
c GENERAL DESCRIPTION:
c   This subroutine provides the numerical flux function
c   defined at the cell interface; this flux is given in
c   the NORMAL DIRECTION (nvect)
c   The low-mach number corrections are made for the flux functions
c
c EQUATIONS: 2D Euler equations of gas dynamics
c
c
c REFERENCE: 1) JCP, 129, 364-382 (1996)
c              " A Sequel to AUSM: AUSM+ ";
c              M.S.Liou
c            2) AIAA Journal, Sept. 1998
c          "Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
c            J.R.Edwards and M.S.Liou
c----------------------------------------------------------------------
c INPUT:
c
c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
c                 ( -3/4 &lt;= alpha &lt;= 3/16 ) (imposed as a parameter)
c
c      beta    -- parameter of the AUSM+ scheme in the Mach function;
c                 ( -1/16 &lt;= beta &lt;= 1/2 )  (imposed as a parameter)
c
c      wvec_l  -- vector of the primitive variables (rho,ux,uy,p) at the
c                 left cell;
c
c      wvec_r  -- vector of the primitive variables (rho,ux,uy,p) at the
c                 right cell;
c
c      nvect   -- normal vector to the interface (2 components in 2D);
c
c      tvect   -- tangential vector to the interface;
c
c      ga      -- ratio of the specific heats (assumed constant)
c----------------------------------------------------------------------
c
c OUTPUT:
c         f    -- numerical flux-function in the NORMAL DIRECTION
c----------------------------------------------------------------------
c Variable 'cans' (set after the descriptions of all the variables)
c  can be used for switching off the low-Mach number additions
c  by simply assigning 'cans=0.0d0'
c----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
       integer nesp,i
       real*8 gamg,rog,pg,ung,utg
       real*8 gamd,rod,pd,und,utd
       real*8 alpha,beta
       real*8 gm1l,gm1r,f(*)
       real*8 un_l, un_r, ut_l, ut_r
       real*8 ml,mr,Mplus,Mmin,mmid
       real*8 mpl_m, mmin_m,am
       real*8 rold_l,pold_l,eold_l
       real*8 rold_r,pold_r,eold_r
       real*8 Pplus,Pmin,pmid
       real*8 hr_l,hr_r,top,bot
       real*8 br1,br2,rum
       real*8 aleft, arigh,v_inf
       real*8 epsil,qq,amw,Mmin1,Mplus1
       real*8 fmid,mlw,mrw,termp
       real*8 ur_r,ur_l,urm,mhalf,mhalfr
       real*8 canc,cellt
       real*8 yg(*),yd(*)
       parameter(alpha = 0.1875D0, beta = 0.125D0)
       parameter(epsil = 1.0d0)
       canc=1.0d0
c-------------------------------------------------------------
       gm1l=gamg-1.0d0
       gm1r=gamd-1.0d0
c----------------------------
       rold_l=rog
       pold_l=pg
c----------------------------
       rold_r=rod
       pold_r=pd
c------------------------------------------------------------------
c Computation of the specific total energy on the left and right.
c------------------------------------------------------------------
       eold_l=(ung*ung+utg*utg)/2.0d0
       eold_l=eold_l+pold_l/(gm1l*rold_l)
       eold_r=(und*und+utd*utd)/2.0d0
       eold_r=eold_r+pold_r/(gm1r*rold_r)
c------------------------------------------------------------------
c   Computing reference velocity and its derivatives
c   see Eq.(2) of the Ref.2).
c------------------------------------------------------------------
       aleft=sqrt(gamg*pold_l/rold_l)
       arigh=sqrt(gamd*pold_r/rold_r)
       qq=sqrt(ung*ung+utg*utg)
       if(qq .lt. (epsil*v_inf)) then
          ur_l = epsil*v_inf
       else
          ur_l=qq
       endif
c-----------------------------
       if(ur_l .ge. aleft) then
          ur_l=aleft
       endif
c-------------------------------------------------------------------
       qq=sqrt(und*und+utd*utd)
       if(qq .lt. (epsil*v_inf)) then
          ur_r = epsil*v_inf
       else
          ur_r=qq
       endif
c-------------------------
       if(ur_r .ge. arigh) then
          ur_r=arigh
       endif
c-------------------------------------------------------------------
c  Reference velocity at the interface is taken as an average
c  of the reference velocities of the neighbouring cells
c-------------------------------------------------------------------
       urm=0.5d0*(ur_l+ur_r)
c-------------------------------------------------------------------
c Computation of the speed of sound;
c numerical speed of sound at the interface is taken as an average
c of the speeds of sounds of the neighbouring cells
c-------------------------------------------------------------------
       am=0.5d0*(aleft+arigh)
c-------------------------------------------------------------------
c Computing numerical Mach number; see p.370, under (A1).
c-------------------------------------------------------------------
       un_l=ung
       un_r=und
       ut_l=utg
       ut_r=utd
c-------------------------------------------------------------------
       ml=un_l/am
       mr=un_r/am
       mhalf=0.5d0*(un_l+un_r)/am
c-------------------------------
       mhalfr=urm/am
c-------------------------------------------------------------------
c Scaling function for the speed of sound;see Eq.(32) of the Ref. 2)
c-------------------------------------------------------------------
       top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
       top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
       bot=1.0d0+mhalfr*mhalfr
       if(abs(canc-0.0d0).lt.0.000001d0) then
         fmid=1.0d0
       else
         fmid=sqrt(top)/bot
       endif
c------------------------------------------------------------------
c 'New' speed of sound 'amw' defined as a product of the scaling
c  function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
c------------------------------------------------------------------
       amw=fmid*am
       mlw=un_l/amw
       mrw=un_r/amw
c--------------------------
       am=amw
c-------------------------------------------------------------------
c  Redefinition of the numerical mach numbers
c  See Eqs.(33) and (34) of the Ref. 2)
c-------------------------------------------------------------------
       if(abs(canc-0.0d0).lt.0.000001d0) then
         top=2.0d0
         bot=0.0d0
       else
         top=1.0d0+mhalfr*mhalfr
         bot=1.0d0-mhalfr*mhalfr
       endif
       ml=0.5d0*(top*mlw+bot*mrw)
       mr=0.5d0*(top*mrw+bot*mlw)
c-------------------------------------------------------------------
c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
c  see (19a) and (19b), p.367. of the Ref.1)
c-------------------------------------------------------------------
       if(abs(ml) .ge. 1.0d0) then
         Mplus=(ml+abs(ml))/2.0d0
       else
         Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
         Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
       endif
       Mplus1=(ml+abs(ml))/2.0d0
c-------------------------------------------------------------------
       if(abs(mr) .ge. 1.0d0) then
         Mmin=(mr-abs(mr))/2.0d0
       else
         Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
         Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
       endif
       Mmin1=(mr-abs(mr))/2.0d0
c-----------------------------------------------------------
c mmid is m_{1/2} (notation as in the paper),
c see Eq.(13), p.366 of the Ref.1)
c-----------------------------------------------------------
       mmid=Mplus+Mmin
c----------------------------------------------------------------
c computing the main convective variables
c mpl_m is m^{+}_{1/2} (paper's notation) and
c mmin_m  is m^{-}_{1/2} (paper's notation),
c see Eq. (A2) on p.370 of the Ref.1)
c----------------------------------------------------------------
       termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
       termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
c-------------------------------------------------------------
       if(mmid .ge. 0.0d0) then
         mpl_m = mmid
       else
         mpl_m = 0.0d0
       endif
c------------------------------------------------------------------
       if(mmid .le. 0.0d0) then
         mmin_m = mmid
       else
         mmin_m = 0.0d0
       endif
c---------------------------------------------------------------
c Computing the calligraphic P+ and P- with their derivatives,
c see (21a) & (21b) on p.368 of the Ref.1)
c---------------------------------------------------------------
       if(ml .ge. 1.0d0) then
         Pplus = 1.0d0
       else
         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
           Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
           Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
         else
           Pplus = 0.0d0
         endif
       endif
c---------------------------------------------------------------
       if(mr .ge. 1.0d0) then
         Pmin = 0.0d0
       else
         if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
           Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
           Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
         else
           Pmin = 1.0d0
         endif
       endif
c-------------------------------------------------------------------
c computing pmid - p_{1/2}, see Eq.(20b), p.367 of the Ref.1)
c-------------------------------------------------------------------
       pmid=Pplus*pold_l + Pmin*pold_r
c---------------------------------------------------------------------
       rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
c---------------------------------------------------------------------
c Computing numerical fluxes
c---------------------------------------------------------------------
       f(1)=rum
c---------------------------------------------------------------------
       if(rum .ge. 0.0d0) then
         br1=rum*un_l
       else
         br1=rum*un_r
       endif
       f(2)=br1+pmid
c-------------------------------------------------------------
       if(rum .ge. 0.0d0) then
         br2=rum*ut_l
       else
         br2=rum*ut_r
       endif
       f(3)=br2
c-------------------------------------------------------------
       hr_l=(ung*ung+utg*utg)/2.0d0+gamg*pold_l/gm1l/rold_l
       hr_r=(und*und+utd*utd)/2.0d0+gamd*pold_r/gm1r/rold_r
       if(rum .ge. 0.0d0) then
         f(4)=rum*hr_l
       else
         f(4)=rum*hr_r
       endif
c---------------------------------------------------------------------
       do 777 i=1,nesp
        if(rum .ge. 0.0d0) then
          br1=rum*yg(i)
        else
          br1=rum*yd(i)
        endif
        f(4+i)=br1
 777   continue
c----------------------------------------------------------------------
       cellt=1.0d0/(0.5d0*abs(un_l+un_r)+am)
c----------------------------------------------------------------------
       return
       end