C FAU_UP SOURCE BECC 10/05/05 21:15:07 6674 SUBROUTINE FAU_UP(nesp, & gaml,rhol,pl,unl,utl, & gamr,rhor,pr,unr,utr, & yl,yr,v_inf,f, & cellt) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : FAUSM2 ('FLUX for Bas Mach') C C DESCRIPTION : Voir C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. Beccantini, 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 C All Speeds" c J.R.Edwards and M.S.Liou C 3) JCP, 214, 137-170 (2006) C A Sequel to AUSM, Part 2 AUSM+-up for all speeds c M.S.Liou C c---------------------------------------------------------------------- c INPUT: c implicit integer(i-n) integer nesp real*8 gaml, rhol, pl, unl, utl & , gamr, rhor, pr, unr, utr & , yl(*), yr(*), v_inf, f(*) & , cellt C Some functions and parameters real*8 fau_m4, fau_p5, kp, ku, sigma C real*8 gm1l, gm1r, htotl, htotr, asonl, asonr, qq & , uco_l, uco_r, mco0, asonm, ml, mr, mbar2 & , mzero, fson, alpha, rhom, mhalf, uhalf & , fp5p, fp5m, phalf PARAMETER(kp=0.25d0, ku=0.75d0, sigma=1.0D0) C PARAMETER(kp=0.25d0, ku=0.0d0, sigma=1.0D0) c------------------------------------------------------------- gm1l=gaml-1.0d0 gm1r=gamr-1.0d0 c------------------------------------------------------------------ c Computation of the specific total enthalpies. c------------------------------------------------------------------ htotl = ((unl*unl) + (utl*utl)) / 2.0d0 htotl = htotl + (pl / (gm1l*rhol)) htotl = htotl + (pl / rhol) htotr = ((unr*unr) + (utr*utr)) / 2.0d0 htotr = htotr + (pr / (gm1r*rhor)) htotr = htotr + (pr / rhor) * write(*,*) htotl, htotr c------------------------------------------------------------------ c Computing reference velocity c see Eq.(70) of the Ref.3). c------------------------------------------------------------------ asonl=sqrt(gaml*pl/rhol) asonr=sqrt(gamr*pr/rhor) qq=sqrt((unl*unl)+(utl*utl)) c------------------------------------------------------------------ if(qq .lt. v_inf) then uco_l = v_inf else uco_l = qq endif c----------------------------- if(uco_l .ge. asonl) then uco_l = asonl endif c------------------------------------------------------------------- qq=sqrt((unr*unr)+(utr*utr)) if(qq .lt. v_inf) then uco_r = v_inf else uco_r = qq endif c------------------------- if(uco_r .ge. asonr) then uco_r = asonr endif c------------------------------------------------------------------- c mco0 is the equivalent of M_{\inf} in Ref 3. c------------------------------------------------------------------- mco0 = max((uco_r / asonr), (uco_l / asonl)) c write(*,*) 'mco0 = ', mco0 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------------------------------------------------------------------- asonm=0.5d0*(asonl + asonr) c------------------------------------------------------------------- c Computing numerical Mach number c ml, mr, mbar2 (used in the term for Mu, which does not involves c the cut-off speed, see equation (70), ref 3.) c Then, even if mc0 = 1, the term involving kp acts! c------------------------------------------------------------------- ml = unl / asonm mr = unr / asonm mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm) * write(*,*) asonm * write(*,*) ml, mr, mbar2 c------------------------------- c------------------------------------------------------------------- c Scaling function for the speed of sound c We redefine mzero (see (71), ref 3). c------------------------------------------------------------------- C mzero = max((mco0 * mco0), mbar2) mzero = min(1.0D0, mzero) mzero = mzero ** 0.5D0 C fson = mzero * (2.0D0 - mzero) alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson)) C write(*,*) 'alpha = ', alpha C------------------------------------------------------------------- C mhalf (see (73), ref 3) C First, contribution of Mu; then contribution of ml and mr. C------------------------------------------------------------------- C rhom = 0.5D0 * (rhol + rhor) mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0) mhalf = -1.0D0 * mhalf * kp * (pr - pl) / & (fson * rhom * asonm * asonm) C mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0) c write(*,*) 'mhalf = ', mhalf uhalf = mhalf * asonm c write(*,*) 'unl, unr ' , unl, unr c write(*,*) 'mhalf, uhalf ', mhalf, uhalf C------------------------------------------------------------------- C phalf C------------------------------------------------------------------- fp5p = fau_p5(ml,1.0D0,alpha) fp5m = fau_p5(mr,-1.0D0,alpha) phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm * & (unr - unl) phalf = phalf + (fp5p * pl) + (fp5m * pr) c write(*,*) 'phalf = ', phalf c--------------------------------------------------------------------- c Computing numerical fluxes c--------------------------------------------------------------------- if(uhalf .ge. 0.0d0) then f(1) = uhalf * rhol f(2) = (f(1) * unl) + phalf f(3) = f(1) * utl f(4) = f(1) * htotl do iesp = 1, nesp, 1 f(4 + iesp) = f(1) * yl(iesp) enddo else f(1) = uhalf * rhor f(2) = (f(1) * unr) + phalf f(3) = f(1) * utr f(4) = f(1) * htotr do iesp = 1, nesp, 1 f(4 + iesp) = f(1) * yr(iesp) enddo endif c---------------------------------------------------------------------- cellt = 1.0d0 / ((0.5d0*abs(unl+unr))+asonm) c write(*,*) 'cellt = ', cellt c---------------------------------------------------------------------- return end