1 : C FAU_UP SOURCE BECC 10/05/05 21:15:07 6674
2 : SUBROUTINE FAU_UP(nesp,
3 : & gaml,rhol,pl,unl,utl,
4 : & gamr,rhor,pr,unr,utr,
5 : & yl,yr,v_inf,f,
6 : & cellt)
7 : C************************************************************************
8 : C
9 : C PROJET : CASTEM 2000
10 : C
11 : C NOM : FAUSM2 ('FLUX for Bas Mach')
12 : C
13 : C DESCRIPTION : Voir
14 : C
15 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
16 : C
17 : C AUTEUR : A. Beccantini, DM2S/SFME/LTMF
18 : C
19 : C************************************************************************
20 : C
21 : c----------------------------------------------------------------------
22 : c GENERAL DESCRIPTION:
23 : c This subroutine provides the numerical flux function
24 : c defined at the cell interface; this flux is given in
25 : c the NORMAL DIRECTION (nvect)
26 : c The low-mach number corrections are made for the flux functions
27 : c
28 : c EQUATIONS: 2D Euler equations of gas dynamics
29 : c
30 : c
31 : c REFERENCE: 1) JCP, 129, 364-382 (1996)
32 : c " A Sequel to AUSM: AUSM+ ";
33 : c M.S.Liou
34 : c 2) AIAA Journal, Sept. 1998
35 : c "Low-Diffusion Flux-Splitting Methods for Flows at
36 : C All Speeds"
37 : c J.R.Edwards and M.S.Liou
38 : C 3) JCP, 214, 137-170 (2006)
39 : C A Sequel to AUSM, Part 2 AUSM+-up for all speeds
40 : c M.S.Liou
41 : C
42 : c----------------------------------------------------------------------
43 : c INPUT:
44 : c
45 : implicit integer(i-n)
46 : integer nesp
47 : real*8 gaml, rhol, pl, unl, utl
48 : & , gamr, rhor, pr, unr, utr
49 : & , yl(*), yr(*), v_inf, f(*)
50 : & , cellt
51 : C Some functions and parameters
52 : real*8 fau_m4, fau_p5, kp, ku, sigma
53 : C
54 : real*8 gm1l, gm1r, htotl, htotr, asonl, asonr, qq
55 : & , uco_l, uco_r, mco0, asonm, ml, mr, mbar2
56 : & , mzero, fson, alpha, rhom, mhalf, uhalf
57 : & , fp5p, fp5m, phalf
58 :
59 : PARAMETER(kp=0.25d0, ku=0.75d0, sigma=1.0D0)
60 : C PARAMETER(kp=0.25d0, ku=0.0d0, sigma=1.0D0)
61 : c-------------------------------------------------------------
62 : gm1l=gaml-1.0d0
63 : gm1r=gamr-1.0d0
64 : c------------------------------------------------------------------
65 : c Computation of the specific total enthalpies.
66 : c------------------------------------------------------------------
67 : htotl = ((unl*unl) + (utl*utl)) / 2.0d0
68 : htotl = htotl + (pl / (gm1l*rhol))
69 : htotl = htotl + (pl / rhol)
70 : htotr = ((unr*unr) + (utr*utr)) / 2.0d0
71 : htotr = htotr + (pr / (gm1r*rhor))
72 : htotr = htotr + (pr / rhor)
73 : * write(*,*) htotl, htotr
74 : c------------------------------------------------------------------
75 : c Computing reference velocity
76 : c see Eq.(70) of the Ref.3).
77 : c------------------------------------------------------------------
78 : asonl=sqrt(gaml*pl/rhol)
79 : asonr=sqrt(gamr*pr/rhor)
80 : qq=sqrt((unl*unl)+(utl*utl))
81 : c------------------------------------------------------------------
82 : if(qq .lt. v_inf) then
83 : uco_l = v_inf
84 : else
85 : uco_l = qq
86 : endif
87 : c-----------------------------
88 : if(uco_l .ge. asonl) then
89 : uco_l = asonl
90 : endif
91 : c-------------------------------------------------------------------
92 : qq=sqrt((unr*unr)+(utr*utr))
93 : if(qq .lt. v_inf) then
94 : uco_r = v_inf
95 : else
96 : uco_r = qq
97 : endif
98 : c-------------------------
99 : if(uco_r .ge. asonr) then
100 : uco_r = asonr
101 : endif
102 : c-------------------------------------------------------------------
103 : c mco0 is the equivalent of M_{\inf} in Ref 3.
104 : c-------------------------------------------------------------------
105 : mco0 = max((uco_r / asonr), (uco_l / asonl))
106 : c write(*,*) 'mco0 = ', mco0
107 : c-------------------------------------------------------------------
108 : c Computation of the speed of sound;
109 : c numerical speed of sound at the interface is taken as an average
110 : c of the speeds of sounds of the neighbouring cells
111 : c-------------------------------------------------------------------
112 : asonm=0.5d0*(asonl + asonr)
113 : c-------------------------------------------------------------------
114 : c Computing numerical Mach number
115 : c ml, mr, mbar2 (used in the term for Mu, which does not involves
116 : c the cut-off speed, see equation (70), ref 3.)
117 : c Then, even if mc0 = 1, the term involving kp acts!
118 : c-------------------------------------------------------------------
119 : ml = unl / asonm
120 : mr = unr / asonm
121 : mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm)
122 : * write(*,*) asonm
123 : * write(*,*) ml, mr, mbar2
124 : c-------------------------------
125 : c-------------------------------------------------------------------
126 : c Scaling function for the speed of sound
127 : c We redefine mzero (see (71), ref 3).
128 : c-------------------------------------------------------------------
129 : C
130 : mzero = max((mco0 * mco0), mbar2)
131 : mzero = min(1.0D0, mzero)
132 : mzero = mzero ** 0.5D0
133 : C
134 : fson = mzero * (2.0D0 - mzero)
135 : alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson))
136 : C write(*,*) 'alpha = ', alpha
137 : C-------------------------------------------------------------------
138 : C mhalf (see (73), ref 3)
139 : C First, contribution of Mu; then contribution of ml and mr.
140 : C-------------------------------------------------------------------
141 : C
142 : rhom = 0.5D0 * (rhol + rhor)
143 : mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0)
144 : mhalf = -1.0D0 * mhalf * kp * (pr - pl) /
145 : & (fson * rhom * asonm * asonm)
146 : C
147 : mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0)
148 : c write(*,*) 'mhalf = ', mhalf
149 : uhalf = mhalf * asonm
150 : c write(*,*) 'unl, unr ' , unl, unr
151 : c write(*,*) 'mhalf, uhalf ', mhalf, uhalf
152 : C-------------------------------------------------------------------
153 : C phalf
154 : C-------------------------------------------------------------------
155 : fp5p = fau_p5(ml,1.0D0,alpha)
156 : fp5m = fau_p5(mr,-1.0D0,alpha)
157 : phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm *
158 : & (unr - unl)
159 : phalf = phalf + (fp5p * pl) + (fp5m * pr)
160 : c write(*,*) 'phalf = ', phalf
161 : c---------------------------------------------------------------------
162 : c Computing numerical fluxes
163 : c---------------------------------------------------------------------
164 : if(uhalf .ge. 0.0d0) then
165 : f(1) = uhalf * rhol
166 : f(2) = (f(1) * unl) + phalf
167 : f(3) = f(1) * utl
168 : f(4) = f(1) * htotl
169 : do iesp = 1, nesp, 1
170 : f(4 + iesp) = f(1) * yl(iesp)
171 : enddo
172 : else
173 : f(1) = uhalf * rhor
174 : f(2) = (f(1) * unr) + phalf
175 : f(3) = f(1) * utr
176 : f(4) = f(1) * htotr
177 : do iesp = 1, nesp, 1
178 : f(4 + iesp) = f(1) * yr(iesp)
179 : enddo
180 : endif
181 : c----------------------------------------------------------------------
182 : cellt = 1.0d0 / ((0.5d0*abs(unl+unr))+asonm)
183 : c write(*,*) 'cellt = ', cellt
184 : c----------------------------------------------------------------------
185 : return
186 : end
187 :
© Cast3M 2003 - All rights reserved.
Disclaimer