1 : C FAUSM4 SOURCE CHAT 05/01/12 23:56:50 5004
2 : C FAUSM3 SOURCE BECC 01/07/27 21:16:20 4170
3 : SUBROUTINE FAUSM4(NESP,
4 : & GAMG,ROG,PG,UNG,UTG,UVG,
5 : & GAMD,ROD,PD,UND,UTD,UVD,
6 : & YG,YD,V_INF,F,
7 : & CELLT)
8 : C************************************************************************
9 : C
10 : C PROJET : CASTEM 2000
11 : C
12 : C NOM : FAUSM4
13 : C
14 : C DESCRIPTION : Voir KONJA2
15 : C
16 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
17 : C
18 : C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
19 : C
20 : C************************************************************************
21 : C
22 : c----------------------------------------------------------------------
23 : c GENERAL DESCRIPTION:
24 : c This subroutine provides the numerical flux function
25 : c defined at the cell interface; this flux is given in
26 : c the NORMAL DIRECTION (nvect)
27 : c The low-mach number corrections are made for the flux functions
28 : c
29 : c EQUATIONS: 3D Euler equations of gas dynamics
30 : c
31 : c
32 : c REFERENCE: 1) JCP, 129, 364-382 (1996)
33 : c " A Sequel to AUSM: AUSM+ ";
34 : c M.S.Liou
35 : c 2) AIAA Journal, Sept. 1998
36 : c "Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
37 : c J.R.Edwards and M.S.Liou
38 : c----------------------------------------------------------------------
39 : c INPUT:
40 : c
41 : c NESP -- number of species in the multispecies gas
42 : c
43 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
44 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
45 : c
46 : c beta -- parameter of the AUSM+ scheme in the Mach function;
47 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
48 : c
49 : c (gamg,rhog,pg,ung,utg,uvg) -- vector of the primitive variables
50 : c at the left cell;
51 : c
52 : c (gamd,rhod,pd,und,utd,uvd) -- vector of the primitive variables
53 : c at the right cell;
54 : c
55 : c yg -- vector of the mass fractionc of the species
56 : c at the left cell;
57 : c
58 : c yd -- vector of the mass fractions of the species
59 : c at the right cell;
60 : c
61 : c v_inf -- parameter for choosing the reference velocity
62 : c when the magnitude of the physical velocity
63 : c is close to zero
64 : c----------------------------------------------------------------------
65 : c
66 : c OUTPUT:
67 : c f -- numerical flux-function in the NORMAL DIRECTION
68 : c----------------------------------------------------------------------
69 : c Variable 'cans' (set after the descriptions of all the variables)
70 : c can be used for switching off the low-Mach number additions
71 : c by simply assigning 'cans=0.0d0'
72 : c----------------------------------------------------------------------
73 : IMPLICIT INTEGER(I-N)
74 : integer nesp,i
75 : real*8 gamg,rog,pg,ung,utg,uvg
76 : real*8 gamd,rod,pd,und,utd,uvd
77 : real*8 alpha,beta,upr_l,upr_r
78 : real*8 gm1l,gm1r,f(*)
79 : real*8 un_l, un_r, ut_l, ut_r
80 : real*8 uv_l, uv_r
81 : real*8 ml,mr,Mplus,Mmin,mmid
82 : real*8 mpl_m, mmin_m,am
83 : real*8 rold_l,pold_l,eold_l
84 : real*8 rold_r,pold_r,eold_r
85 : real*8 Pplus,Pmin,pmid
86 : real*8 hr_l,hr_r,top,bot
87 : real*8 br1,br2
88 : real*8 aleft, arigh
89 : real*8 epsil,qq,amw,Mmin1,Mplus1
90 : real*8 fmid,mlw,mrw,termp
91 : real*8 ur_r,ur_l,urm,mhalf,mhalfr
92 : real*8 canc,cellt,v_inf,rum
93 : real*8 yg(*),yd(*)
94 : parameter(alpha = 0.1875D0, beta = 0.125D0)
95 : parameter(epsil = 1.0d0)
96 : canc=1.0d0
97 : upr_l=0.0d0
98 : upr_r=0.0d0
99 : c-------------------------------------------------------------
100 : gm1l=gamg-1.0d0
101 : gm1r=gamd-1.0d0
102 : c----------------------------
103 : rold_l=rog
104 : pold_l=pg
105 : c----------------------------
106 : rold_r=rod
107 : pold_r=pd
108 : c------------------------------------------------------------------
109 : c Computation of the specific total energy on the left and right.
110 : c------------------------------------------------------------------
111 : eold_l=(ung*ung+utg*utg+uvg*uvg)/2.0d0
112 : eold_l=eold_l+pold_l/(gm1l*rold_l)
113 : eold_r=(und*und+utd*utd+uvd*uvd)/2.0d0
114 : eold_r=eold_r+pold_r/(gm1r*rold_r)
115 : c------------------------------------------------------------------
116 : c Computing reference velocity and its derivatives
117 : c see Eq.(2) of the Ref.2).
118 : c------------------------------------------------------------------
119 : aleft=sqrt(gamg*pold_l/rold_l)
120 : arigh=sqrt(gamd*pold_r/rold_r)
121 : qq=sqrt(ung*ung+utg*utg+uvg*uvg)
122 : if(qq .lt. (epsil*v_inf)) then
123 : ur_l = epsil*v_inf
124 : else
125 : ur_l=qq
126 : endif
127 : c-------------------------------------------------------------------
128 : if(ur_l .ge. aleft) then
129 : ur_l=aleft
130 : endif
131 : c-------------------------------------------------------------------
132 : if(ur_l .lt. upr_l) ur_l=upr_l
133 : c-------------------------------------------------------------------
134 : qq=sqrt(und*und+utd*utd+uvd*uvd)
135 : if(qq .lt. (epsil*v_inf)) then
136 : ur_r = epsil*v_inf
137 : else
138 : ur_r=qq
139 : endif
140 : c------------------------------------------------------------------
141 : if(ur_r .ge. arigh) then
142 : ur_r=arigh
143 : endif
144 : c------------------------------------------------------------------
145 : if(ur_r .lt. upr_r) ur_r=upr_r
146 : c-------------------------------------------------------------------
147 : c Reference velocity at the interface is taken as an average
148 : c of the reference velocities of the neighbouring cells
149 : c-------------------------------------------------------------------
150 : urm=0.5d0*(ur_l+ur_r)
151 : c-------------------------------------------------------------------
152 : c Computation of the speed of sound;
153 : c numerical speed of sound at the interface is taken as an average
154 : c of the speeds of sounds of the neighbouring cells
155 : c-------------------------------------------------------------------
156 : am=0.5d0*(aleft+arigh)
157 : c-------------------------------------------------------------------
158 : c Computing numerical Mach number; see p.370, under (A1).
159 : c-------------------------------------------------------------------
160 : un_l=ung
161 : un_r=und
162 : c--------------
163 : ut_l=utg
164 : ut_r=utd
165 : c--------------
166 : uv_l=uvg
167 : uv_r=uvd
168 : c-------------------------------------------------------------------
169 : ml=un_l/am
170 : mr=un_r/am
171 : mhalf=0.5d0*(un_l+un_r)/am
172 : c-------------------------------
173 : mhalfr=urm/am
174 : c-------------------------------------------------------------------
175 : c Scaling function for the speed of sound;see Eq.(32) of the Ref. 2)
176 : c-------------------------------------------------------------------
177 : top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
178 : top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
179 : bot=1.0d0+mhalfr*mhalfr
180 : if(abs(canc-0.0d0).lt.0.000001d0) then
181 : fmid=1.0d0
182 : else
183 : fmid=sqrt(top)/bot
184 : endif
185 : c------------------------------------------------------------------
186 : c 'New' speed of sound 'amw' defined as a product of the scaling
187 : c function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
188 : c------------------------------------------------------------------
189 : amw=fmid*am
190 : mlw=un_l/amw
191 : mrw=un_r/amw
192 : c--------------------------
193 : am=amw
194 : c-------------------------------------------------------------------
195 : c Redefinition of the numerical mach numbers
196 : c See Eqs.(33) and (34) of the Ref. 2)
197 : c-------------------------------------------------------------------
198 : if(abs(canc-0.0d0).lt.0.000001d0) then
199 : top=2.0d0
200 : bot=0.0d0
201 : else
202 : top=1.0d0+mhalfr*mhalfr
203 : bot=1.0d0-mhalfr*mhalfr
204 : endif
205 : ml=0.5d0*(top*mlw+bot*mrw)
206 : mr=0.5d0*(top*mrw+bot*mlw)
207 : c-------------------------------------------------------------------
207 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
208 : c see (19a) and (19b), p.367. of the Ref.1)
209 : c-------------------------------------------------------------------
210 : if(abs(ml) .ge. 1.0d0) then
211 : Mplus=(ml+abs(ml))/2.0d0
212 : else
213 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
214 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
215 : endif
216 : Mplus1=(ml+abs(ml))/2.0d0
217 : c-------------------------------------------------------------------
218 : if(abs(mr) .ge. 1.0d0) then
219 : Mmin=(mr-abs(mr))/2.0d0
220 : else
221 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
222 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
223 : endif
224 : Mmin1=(mr-abs(mr))/2.0d0
225 : c-----------------------------------------------------------
226 : c mmid is m_{1/2} (notation as in the paper),
227 : c see Eq.(13), p.366 of the Ref.1)
228 : c-----------------------------------------------------------
229 : mmid=Mplus+Mmin
230 : c----------------------------------------------------------------
231 : c computing the main convective variables
232 : c mpl_m is m^{+}_{1/2} (paper's notation) and
233 : c mmin_m is m^{-}_{1/2} (paper's notation),
234 : c see Eq. (A2) on p.370 of the Ref.1)
235 : c----------------------------------------------------------------
236 : termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
237 : termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
238 : c-------------------------------------------------------------
239 : if(mmid .ge. 0.0d0) then
240 : mpl_m = mmid
241 : else
242 : mpl_m = 0.0d0
243 : endif
244 : c------------------------------------------------------------------
245 : if(mmid .le. 0.0d0) then
246 : mmin_m = mmid
247 : else
248 : mmin_m = 0.0d0
249 : endif
250 : c---------------------------------------------------------------
250 : c Computing the calligraphic P+ and P- with their derivatives,
251 : c see (21a) & (21b) on p.368 of the Ref.1)
252 : c---------------------------------------------------------------
253 : if(ml .ge. 1.0d0) then
254 : Pplus = 1.0d0
255 : else
256 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
257 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
258 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
259 : else
260 : Pplus = 0.0d0
261 : endif
262 : endif
263 : c---------------------------------------------------------------
264 : if(mr .ge. 1.0d0) then
265 : Pmin = 0.0d0
266 : else
267 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
268 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
269 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
270 : else
271 : Pmin = 1.0d0
272 : endif
273 : endif
274 : c-------------------------------------------------------------------
275 : c computing pmid - p_{1/2}, see Eq.(20b), p.367 of the Ref.1)
276 : c-------------------------------------------------------------------
277 : pmid=Pplus*pold_l + Pmin*pold_r
278 : c-------------------------------------------------------------------
279 : rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
280 : c-------------------------------------------------------------------
281 : c Computing numerical fluxes
282 : c-------------------------------------------------------------------
283 : f(1)=rum
284 : c-------------------------------------------------------------------
285 : if(rum .ge. 0.0d0) then
286 : br1=rum*un_l
287 : else
288 : br1=rum*un_r
289 : endif
290 : f(2)=br1+pmid
291 : c-------------------------------------------------------------
292 : if(rum .ge. 0.0d0) then
293 : br2=rum*ut_l
294 : else
295 : br2=rum*ut_r
296 : endif
297 : f(3)=br2
298 : c-------------------------------------------------------------
299 : if(rum .ge. 0.0d0) then
300 : br2=rum*uv_l
301 : else
302 : br2=rum*uv_r
303 : endif
304 : f(4)=br2
305 : c-------------------------------------------------------------
306 : hr_l=(ung*ung+utg*utg+uvg*uvg)/2.0d0+gamg*pold_l/gm1l/rold_l
307 : hr_r=(und*und+utd*utd+uvd*uvd)/2.0d0+gamd*pold_r/gm1r/rold_r
308 : if(rum .ge. 0.0d0) then
309 : f(5)=rum*hr_l
310 : else
311 : f(5)=rum*hr_r
312 : endif
313 : c---------------------------------------------------------------------
314 : do 777 i=1,nesp
315 : if(rum .ge. 0.0d0) then
316 : br1=rum*yg(i)
317 : else
318 : br1=rum*yd(i)
319 : endif
320 : f(5+i)=br1
321 : 777 continue
322 : c----------------------------------------------------------------------
323 : cellt=1.0d0/(0.5d0*abs(un_l+un_r)+am)
324 : c----------------------------------------------------------------------
325 : return
326 : end
327 :
328 :
329 :
330 :
331 :
332 :
333 :
334 :
335 :
336 :
337 :
338 :
339 :
340 :
341 :
342 :
343 :
344 :
345 :
346 :
347 :
348 :
349 :
350 :
351 :
352 :
353 :
354 :
355 :
356 :
357 :
358 :
359 :
360 :
361 :
362 :
363 :
364 :
365 :
366 :
367 :
368 :
369 :
370 :
371 :
372 :
373 :
374 :
375 :
376 :
377 :
378 :
© Cast3M 2003 - All rights reserved.
Disclaimer