1 : C FAUSM2 SOURCE CHAT 05/01/12 23:56:44 5004
2 : SUBROUTINE FAUSM2(NESP,
3 : & GAMG,ROG,PG,UNG,UTG,
4 : & GAMD,ROD,PD,UND,UTD,
5 : & YG,YD,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 KONJA2
14 : C
15 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
16 : C
17 : C AUTEUR : S. KUDRIAKOV, 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 All Speeds"
36 : c J.R.Edwards and M.S.Liou
37 : c----------------------------------------------------------------------
38 : c INPUT:
39 : c
40 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
41 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
42 : c
43 : c beta -- parameter of the AUSM+ scheme in the Mach function;
44 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
45 : c
46 : c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
47 : c left cell;
48 : c
49 : c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
50 : c right cell;
51 : c
52 : c nvect -- normal vector to the interface (2 components in 2D);
53 : c
54 : c tvect -- tangential vector to the interface;
55 : c
56 : c ga -- ratio of the specific heats (assumed constant)
57 : c----------------------------------------------------------------------
58 : c
59 : c OUTPUT:
60 : c f -- numerical flux-function in the NORMAL DIRECTION
61 : c----------------------------------------------------------------------
62 : c Variable 'cans' (set after the descriptions of all the variables)
63 : c can be used for switching off the low-Mach number additions
64 : c by simply assigning 'cans=0.0d0'
65 : c----------------------------------------------------------------------
66 : IMPLICIT INTEGER(I-N)
67 : integer nesp,i
68 : real*8 gamg,rog,pg,ung,utg
69 : real*8 gamd,rod,pd,und,utd
70 : real*8 alpha,beta
71 : real*8 gm1l,gm1r,f(*)
72 : real*8 un_l, un_r, ut_l, ut_r
73 : real*8 ml,mr,Mplus,Mmin,mmid
74 : real*8 mpl_m, mmin_m,am
75 : real*8 rold_l,pold_l,eold_l
76 : real*8 rold_r,pold_r,eold_r
77 : real*8 Pplus,Pmin,pmid
78 : real*8 hr_l,hr_r,top,bot
79 : real*8 br1,br2,rum
80 : real*8 aleft, arigh,v_inf
81 : real*8 epsil,qq,amw,Mmin1,Mplus1
82 : real*8 fmid,mlw,mrw,termp
83 : real*8 ur_r,ur_l,urm,mhalf,mhalfr
84 : real*8 canc,cellt
85 : real*8 yg(*),yd(*)
86 : parameter(alpha = 0.1875D0, beta = 0.125D0)
87 : parameter(epsil = 1.0d0)
88 : canc=1.0d0
89 : c-------------------------------------------------------------
90 : gm1l=gamg-1.0d0
91 : gm1r=gamd-1.0d0
92 : c----------------------------
93 : rold_l=rog
94 : pold_l=pg
95 : c----------------------------
96 : rold_r=rod
97 : pold_r=pd
98 : c------------------------------------------------------------------
99 : c Computation of the specific total energy on the left and right.
100 : c------------------------------------------------------------------
101 : eold_l=(ung*ung+utg*utg)/2.0d0
102 : eold_l=eold_l+pold_l/(gm1l*rold_l)
103 : eold_r=(und*und+utd*utd)/2.0d0
104 : eold_r=eold_r+pold_r/(gm1r*rold_r)
105 : c------------------------------------------------------------------
106 : c Computing reference velocity and its derivatives
107 : c see Eq.(2) of the Ref.2).
108 : c------------------------------------------------------------------
109 : aleft=sqrt(gamg*pold_l/rold_l)
110 : arigh=sqrt(gamd*pold_r/rold_r)
111 : qq=sqrt(ung*ung+utg*utg)
112 : if(qq .lt. (epsil*v_inf)) then
113 : ur_l = epsil*v_inf
114 : else
115 : ur_l=qq
116 : endif
117 : c-----------------------------
118 : if(ur_l .ge. aleft) then
119 : ur_l=aleft
120 : endif
121 : c-------------------------------------------------------------------
122 : qq=sqrt(und*und+utd*utd)
123 : if(qq .lt. (epsil*v_inf)) then
124 : ur_r = epsil*v_inf
125 : else
126 : ur_r=qq
127 : endif
128 : c-------------------------
129 : if(ur_r .ge. arigh) then
130 : ur_r=arigh
131 : endif
132 : c-------------------------------------------------------------------
133 : c Reference velocity at the interface is taken as an average
134 : c of the reference velocities of the neighbouring cells
135 : c-------------------------------------------------------------------
136 : urm=0.5d0*(ur_l+ur_r)
137 : c-------------------------------------------------------------------
138 : c Computation of the speed of sound;
139 : c numerical speed of sound at the interface is taken as an average
140 : c of the speeds of sounds of the neighbouring cells
141 : c-------------------------------------------------------------------
142 : am=0.5d0*(aleft+arigh)
143 : c-------------------------------------------------------------------
144 : c Computing numerical Mach number; see p.370, under (A1).
145 : c-------------------------------------------------------------------
146 : un_l=ung
147 : un_r=und
148 : ut_l=utg
149 : ut_r=utd
150 : c-------------------------------------------------------------------
151 : ml=un_l/am
152 : mr=un_r/am
153 : mhalf=0.5d0*(un_l+un_r)/am
154 : c-------------------------------
155 : mhalfr=urm/am
156 : c-------------------------------------------------------------------
157 : c Scaling function for the speed of sound;see Eq.(32) of the Ref. 2)
158 : c-------------------------------------------------------------------
159 : top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
160 : top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
161 : bot=1.0d0+mhalfr*mhalfr
162 : if(abs(canc-0.0d0).lt.0.000001d0) then
163 : fmid=1.0d0
164 : else
165 : fmid=sqrt(top)/bot
166 : endif
167 : c------------------------------------------------------------------
168 : c 'New' speed of sound 'amw' defined as a product of the scaling
169 : c function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
170 : c------------------------------------------------------------------
171 : amw=fmid*am
172 : mlw=un_l/amw
173 : mrw=un_r/amw
174 : c--------------------------
175 : am=amw
176 : c-------------------------------------------------------------------
177 : c Redefinition of the numerical mach numbers
178 : c See Eqs.(33) and (34) of the Ref. 2)
179 : c-------------------------------------------------------------------
180 : if(abs(canc-0.0d0).lt.0.000001d0) then
181 : top=2.0d0
182 : bot=0.0d0
183 : else
184 : top=1.0d0+mhalfr*mhalfr
185 : bot=1.0d0-mhalfr*mhalfr
186 : endif
187 : ml=0.5d0*(top*mlw+bot*mrw)
188 : mr=0.5d0*(top*mrw+bot*mlw)
189 : c-------------------------------------------------------------------
189 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
190 : c see (19a) and (19b), p.367. of the Ref.1)
191 : c-------------------------------------------------------------------
192 : if(abs(ml) .ge. 1.0d0) then
193 : Mplus=(ml+abs(ml))/2.0d0
194 : else
195 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
196 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
197 : endif
198 : Mplus1=(ml+abs(ml))/2.0d0
199 : c-------------------------------------------------------------------
200 : if(abs(mr) .ge. 1.0d0) then
201 : Mmin=(mr-abs(mr))/2.0d0
202 : else
203 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
204 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
205 : endif
206 : Mmin1=(mr-abs(mr))/2.0d0
207 : c-----------------------------------------------------------
208 : c mmid is m_{1/2} (notation as in the paper),
209 : c see Eq.(13), p.366 of the Ref.1)
210 : c-----------------------------------------------------------
211 : mmid=Mplus+Mmin
212 : c----------------------------------------------------------------
213 : c computing the main convective variables
214 : c mpl_m is m^{+}_{1/2} (paper's notation) and
215 : c mmin_m is m^{-}_{1/2} (paper's notation),
216 : c see Eq. (A2) on p.370 of the Ref.1)
217 : c----------------------------------------------------------------
218 : termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
219 : termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
220 : c-------------------------------------------------------------
221 : if(mmid .ge. 0.0d0) then
222 : mpl_m = mmid
223 : else
224 : mpl_m = 0.0d0
225 : endif
226 : c------------------------------------------------------------------
227 : if(mmid .le. 0.0d0) then
228 : mmin_m = mmid
229 : else
230 : mmin_m = 0.0d0
231 : endif
232 : c---------------------------------------------------------------
232 : c Computing the calligraphic P+ and P- with their derivatives,
233 : c see (21a) & (21b) on p.368 of the Ref.1)
234 : c---------------------------------------------------------------
235 : if(ml .ge. 1.0d0) then
236 : Pplus = 1.0d0
237 : else
238 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
239 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
240 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
241 : else
242 : Pplus = 0.0d0
243 : endif
244 : endif
245 : c---------------------------------------------------------------
246 : if(mr .ge. 1.0d0) then
247 : Pmin = 0.0d0
248 : else
249 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
250 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
251 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
252 : else
253 : Pmin = 1.0d0
254 : endif
255 : endif
256 : c-------------------------------------------------------------------
257 : c computing pmid - p_{1/2}, see Eq.(20b), p.367 of the Ref.1)
258 : c-------------------------------------------------------------------
259 : pmid=Pplus*pold_l + Pmin*pold_r
260 : c---------------------------------------------------------------------
261 : rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
262 : c---------------------------------------------------------------------
263 : c Computing numerical fluxes
264 : c---------------------------------------------------------------------
265 : f(1)=rum
266 : c---------------------------------------------------------------------
267 : if(rum .ge. 0.0d0) then
268 : br1=rum*un_l
269 : else
270 : br1=rum*un_r
271 : endif
272 : f(2)=br1+pmid
273 : c-------------------------------------------------------------
274 : if(rum .ge. 0.0d0) then
275 : br2=rum*ut_l
276 : else
277 : br2=rum*ut_r
278 : endif
279 : f(3)=br2
280 : c-------------------------------------------------------------
281 : hr_l=(ung*ung+utg*utg)/2.0d0+gamg*pold_l/gm1l/rold_l
282 : hr_r=(und*und+utd*utd)/2.0d0+gamd*pold_r/gm1r/rold_r
283 : if(rum .ge. 0.0d0) then
284 : f(4)=rum*hr_l
285 : else
286 : f(4)=rum*hr_r
287 : endif
288 : c---------------------------------------------------------------------
289 : do 777 i=1,nesp
290 : if(rum .ge. 0.0d0) then
291 : br1=rum*yg(i)
292 : else
293 : br1=rum*yd(i)
294 : endif
295 : f(4+i)=br1
296 : 777 continue
297 : c----------------------------------------------------------------------
298 : cellt=1.0d0/(0.5d0*abs(un_l+un_r)+am)
299 : c----------------------------------------------------------------------
300 : return
301 : end
302 :
303 :
304 :
305 :
306 :
307 :
308 :
309 :
310 :
311 :
312 :
313 :
314 :
315 :
316 :
317 :
318 :
319 :
320 :
321 :
322 :
323 :
324 :
325 :
326 :
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 :
© Cast3M 2003 - All rights reserved.
Disclaimer