1 : C AUSMF SOURCE CHAT 05/01/12 21:29:07 5004
2 : c----------------------------------------------------------------------
3 : c GENERAL DESCRIPTION:
4 : c This subroutine provides the jacobian which is the derivatives
5 : c of the numerical flux function of the AUSM+ scheme
6 : c defined at the cell interface
7 : c with respect to the primitive variables of the left
8 : c cell (relative to the cell interface).
9 : c
10 : c EQUATIONS: 2D Euler equations of gas dynamics
11 : c
12 : c
13 : c REFERENCE: JCP, 129, 364-382 (1996)
14 : c " A Sequel to AUSM: AUSM+ ".
15 : c----------------------------------------------------------------------
16 : c INPUT:
17 : c
18 : c
19 : c (rog,uxg,uyg,pg) -- vector of the primitive variables at the
20 : c left cell;
21 : c
22 : c (rod,uxd,uyd,pd) -- vector of the primitive variables at the
23 : c right cell;
24 : c
25 : c (cnx,cny) -- normal vector to the interface (2 components in 2D);
26 : c
27 : c (ctx,cty) -- tangential vector to the interface;
28 : c
29 : c ga -- ratio of the specific heats (assumed constant)
30 : c----------------------------------------------------------------------
31 : c
32 : c OUTPUT:
33 : c
34 : c dfro -- vector of length 4 which represents derivatives of
35 : c flux-function for conservative variable (rho)
36 : c with respect to the primitive variables
37 : c
38 : c dfrun -- vector of length 4 which represents derivatives of
39 : c flux-function for conservative variable (rho u_n)
40 : c with respect to the primitive variables
41 : c
42 : c dfrut -- vector of length 4 which represents derivatives of
43 : c flux-function for conservative variable (rho u_t)
44 : c with respect to the primitive variables
45 : c
46 : c dfret -- vector of length 4 which represents derivatives of
47 : c flux-function for conservative variable (rho e_t)
48 : c with respect to the primitive variables
49 : c----------------------------------------------------------------------
50 : SUBROUTINE AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
51 : & ga,cnx,cny,ctx,cty,
52 : & dfro,dfrun,dfrut,dfret)
53 : IMPLICIT INTEGER(I-N)
54 : real*8 rog,uxg,uyg,pg
55 : real*8 rod,uxd,uyd,pd
56 : real*8 cnx,cny,ctx,cty
57 : real*8 jl(4,4),f(4)
58 : real*8 dfro(4),dfrun(4)
59 : real*8 dfrut(4),dfret(4)
60 : real*8 alpha,beta
61 : real*8 ga, gm1,temph
62 : real*8 n_x,n_y,t_x,t_y
63 : real*8 un_l, un_r, ut_l, ut_r
64 : real*8 ml,mr,Mplus,Mmin,mmid
65 : real*8 mpl_m, mmin_m,am
66 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
67 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
68 : real*8 Pplus,Pmin,pmid
69 : real*8 hr_l,hr_r
70 : real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
71 : real*8 aleft, arigh
72 : real*8 damr_l,damu_l
73 : real*8 damv_l,damp_l
74 : real*8 dmlr_l,dmlu_l
75 : real*8 dmlv_l,dmlp_l
76 : real*8 dmrr_l,dmru_l
77 : real*8 dmrv_l,dmrp_l
78 : real*8 dMpr_l,dMpu_l
79 : real*8 dMpv_l,dMpp_l
80 : real*8 dMmr_l,dMmu_l
81 : real*8 dMmv_l,dMmp_l
82 : real*8 dmir_l,dmiu_l
83 : real*8 dmiv_l,dmip_l
84 : real*8 d3mr_l,d3mu_l
85 : real*8 d3mv_l,d3mp_l
86 : real*8 d2mr_l,d2mu_l
87 : real*8 d2mv_l,d2mp_l
88 : real*8 dPpr_l,dPpu_l
89 : real*8 dPpv_l,dPpp_l
90 : real*8 dPmr_l,dPmu_l
91 : real*8 dPmv_l,dPmp_l
92 : real*8 dpir_l,dpiu_l
93 : real*8 dpiv_l,dpip_l
94 : integer i
95 : c
96 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
97 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
98 : c
99 : c beta -- parameter of the AUSM+ scheme in the Mach function;
100 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
101 : c
102 : parameter(alpha = 0.1875D0, beta = 0.125D0)
103 : c-----------------------
104 : gm1=ga-1.0d0
105 : c-----------------------
106 : n_x=cnx
107 : n_y=cny
108 : t_x=ctx
109 : t_y=cty
110 : c----------------------------
111 : c----------------------------
112 : rold_l=rog
113 : uold_l=uxg
114 : vold_l=uyg
115 : pold_l=pg
116 : c-----------------------
117 : rold_r=rod
118 : uold_r=uxd
119 : vold_r=uyd
120 : pold_r=pd
121 : c------------------------------------------------------------------
122 : c Computation of the specific total energy on the left and right.
123 : c------------------------------------------------------------------
124 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
125 : eold_l=eold_l+pold_l/(gm1*rold_l)
126 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
127 : eold_r=eold_r+pold_r/(gm1*rold_r)
128 : c-------------------------------------------------------------------
129 : c Computation of the speed of sound and its derivatives;
130 : c numerical speed of sound at the interface is taken as an average
131 : c of the speeds of sounds of the neighbouring cells
132 : c-------------------------------------------------------------------
133 : aleft=dsqrt(ga*pold_l/rold_l)
134 : arigh=dsqrt(ga*pold_r/rold_r)
135 : am=0.5d0*(aleft+arigh)
136 : c-------------------------------------------------------------------
137 : damr_l=-aleft/(4.0d0*rold_l)
138 : damu_l=0.0d0
139 : damv_l=0.0d0
140 : damp_l=ga/(4.0d0*aleft*rold_l)
141 : c-------------------------------------------------------------------
142 : c Computing numerical Mach number and its derivatives,
143 : c see p.370, under (A1).
144 : c-------------------------------------------------------------------
145 : un_l=uold_l*n_x+vold_l*n_y
146 : un_r=uold_r*n_x+vold_r*n_y
147 : ut_l=uold_l*t_x+vold_l*t_y
148 : ut_r=uold_r*t_x+vold_r*t_y
149 : c-------------------------------------------------------------------
150 : ml=un_l/am
151 : mr=un_r/am
152 : c-------------------------------------------------------------------
152 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
153 : c see (19a) and (19b), p.367.
154 : c-------------------------------------------------------------------
155 : if(dabs(ml) .ge. 1.0d0) then
156 : Mplus=(ml+dabs(ml))/2.0d0
157 : else
158 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
159 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
160 : endif
161 : c-------------------------------------------------------------------
162 : if(dabs(mr) .ge. 1.0d0) then
163 : Mmin=(mr-dabs(mr))/2.0d0
164 : else
165 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
166 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
167 : endif
168 : c-------------------------------------------------------------------
169 : c Derivatives of ml and mr with respect to both sets of primitive
170 : c variables: left and right.
171 : c-------------------------------------------------------------------
172 : temp_l=-un_l/(am*am)
173 : temp_r=-un_r/(am*am)
174 : c--------
175 : dmlr_l=temp_l*damr_l
176 : dmrr_l=temp_r*damr_l
177 : c--------
178 : dmlu_l=n_x/am+temp_l*damu_l
179 : dmru_l=temp_r*damu_l
180 : c--------
181 : dmlv_l=n_y/am+temp_l*damv_l
182 : dmrv_l=temp_r*damv_l
183 : c--------
184 : dmlp_l=temp_l*damp_l
185 : dmrp_l=temp_r*damp_l
186 : c-----------------------------------------------------------
187 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
188 : c-----------------------------------------------------------
189 : mmid=Mplus+Mmin
190 : c-----------------------------------------------------------
191 : c Computing the derivatives of M+ and M-
192 : c-----------------------------------------------------------
193 : if(ml .ge. 1.0d0) then
194 : dMpr_l=dmlr_l
195 : dMpu_l=dmlu_l
196 : dMpv_l=dmlv_l
197 : dMpp_l=dmlp_l
198 : c--------------------
199 : else
200 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
201 : temph=(ml+1.0d0)/2.0d0
202 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
203 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
204 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
205 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
206 : c--------------------
207 : else
208 : dMpr_l=0.0d0
209 : dMpu_l=0.0d0
210 : dMpv_l=0.0d0
211 : dMpp_l=0.0d0
212 : c---------------------
213 : endif
214 : endif
215 : c-----------------------------------------------------------
216 : if(mr .ge. 1.0d0) then
217 : dMmr_l=0.0d0
218 : dMmu_l=0.0d0
219 : dMmv_l=0.0d0
220 : dMmp_l=0.0d0
221 : c---------------------
222 : else
223 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
224 : temph=(-mr+1.0d0)/2.0d0
225 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
226 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
227 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
228 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
229 : c--------------------
230 : else
231 : dMmr_l=dmrr_l
232 : dMmu_l=dmru_l
233 : dMmv_l=dmrv_l
234 : dMmp_l=dmrp_l
235 : c--------------------
236 : endif
237 : endif
238 : c-----------------------------------------------------------------
239 : c computing the derivatives of m_{1/2} (notation as in the paper)
240 : c-----------------------------------------------------------------
241 : dmir_l=dMpr_l+dMmr_l
242 : c-------------
243 : dmiu_l=dMpu_l+dMmu_l
244 : c-------------
245 : dmiv_l=dMpv_l+dMmv_l
246 : c-------------
247 : dmip_l=dMpp_l+dMmp_l
248 : c----------------------------------------------------------------
249 : c computing the main convective variables and their derivatives
250 : c mpl_m is m^{+}_{1/2} (paper's notation) and
251 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
252 : c----------------------------------------------------------------
253 : if(mmid .ge. 0.0d0) then
254 : mpl_m = mmid
255 : d2mr_l=dmir_l
256 : d2mu_l=dmiu_l
257 : d2mv_l=dmiv_l
258 : d2mp_l=dmip_l
259 : c------------
260 : else
261 : mpl_m = 0.0d0
262 : d2mr_l=0.0d0
263 : d2mu_l=0.0d0
264 : d2mv_l=0.0d0
265 : d2mp_l=0.0d0
266 : c------------
267 : endif
268 : c---------------------------------------------
269 : if(mmid .le. 0.0d0) then
270 : mmin_m = mmid
271 : d3mr_l=dmir_l
272 : d3mu_l=dmiu_l
273 : d3mv_l=dmiv_l
274 : d3mp_l=dmip_l
275 : c------------
276 : else
277 : mmin_m = 0.0d0
278 : d3mr_l=0.0d0
279 : d3mu_l=0.0d0
280 : d3mv_l=0.0d0
281 : d3mp_l=0.0d0
282 : c------------
283 : endif
284 : c---------------------------------------------------------------
284 : c Computing the calligraphic P+ and P- with their derivatives,
285 : c see (21a) & (21b) on p.368.
286 : c---------------------------------------------------------------
287 : if(ml .ge. 1.0d0) then
288 : Pplus = 1.0d0
289 : else
290 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
291 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
292 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
293 : else
294 : Pplus = 0.0d0
295 : endif
296 : endif
297 : c---------------------------------------------------------------
298 : if(mr .ge. 1.0d0) then
299 : Pmin = 0.0d0
300 : else
301 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
302 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
303 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
304 : else
305 : Pmin = 1.0d0
306 : endif
307 : endif
308 : c---------------------------------------------------------------
309 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
310 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
311 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
312 : c--------------
313 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
314 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
315 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
316 : c---------------------------------------------------------------
317 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
318 : dPpr_l=brac_l*dmlr_l
319 : dPpu_l=brac_l*dmlu_l
320 : dPpv_l=brac_l*dmlv_l
321 : dPpp_l=brac_l*dmlp_l
322 : else
323 : dPpr_l=0.0d0
324 : dPpu_l=0.0d0
325 : dPpv_l=0.0d0
326 : dPpp_l=0.0d0
327 : endif
328 : c---------------------------------------------------------------
329 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
330 : dPmr_l=brac_r*dmrr_l
331 : dPmu_l=brac_r*dmru_l
332 : dPmv_l=brac_r*dmrv_l
333 : dPmp_l=brac_r*dmrp_l
334 : else
335 : dPmr_l=0.0d0
336 : dPmu_l=0.0d0
337 : dPmv_l=0.0d0
338 : dPmp_l=0.0d0
339 : endif
340 : c-------------------------------------------------------------------
341 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
342 : c-------------------------------------------------------------------
343 : pmid=Pplus*pold_l + Pmin*pold_r
344 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
345 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
346 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
347 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
348 : c---------------------------------------------------------------------
349 : c computing JACOBIAN as a derivative of the numerical flux function
350 : c with respect to the primitive variables.
351 : c Notation: jl(i,j) --- is the derivative of the i-component of the
352 : c flux function with respect to the j-component of the
353 : c vector of primitive variables of the left state.
354 : c jr(i,j) --- is the derivative of the i-component of the
355 : c flux function with respect to the j-component of the
356 : c vector of primitive variables of the right state.
357 : c---------------------------------------------------------------------
358 : f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
359 : c---------------------------------------------------------------------
360 : jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
361 : jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
362 : jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
363 : jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
364 : jl(1,4)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
365 : c------------------------------------
366 : c------------------------------------
367 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
368 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
369 : f(2)=n_x*(am*br1+pmid)+am*t_x*br2
370 : c---------------------------------------------------------
371 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
372 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
373 : jl(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
374 : jl(2,1)=jl(2,1)+t_x*damr_l*br2+am*t_x*br4
375 : c-------------------
376 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*n_x
377 : br3=br3+d3mu_l*rold_r*un_r
378 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*t_x
379 : br4=br4+d3mu_l*rold_r*ut_r
380 : jl(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
381 : jl(2,2)=jl(2,2)+t_x*damu_l*br2+am*t_x*br4
382 : c-------------------
383 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*n_y
384 : br3=br3+d3mv_l*rold_r*un_r
385 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*t_y
386 : br4=br4+d3mv_l*rold_r*ut_r
387 : jl(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
388 : jl(2,3)=jl(2,3)+t_x*damv_l*br2+am*t_x*br4
389 : c-------------------
390 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
391 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
392 : jl(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
393 : jl(2,4)=jl(2,4)+t_x*damp_l*br2+am*t_x*br4
394 : c-------------------------------------------------------------
395 : c-------------------------------------------------------------
396 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
397 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
398 : f(3)=n_y*(am*br1+pmid)+am*t_y*br2
399 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
400 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
401 : jl(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
402 : jl(3,1)=jl(3,1)+t_y*damr_l*br2+am*t_y*br4
403 : c-------------------
404 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*n_x
405 : br3=br3+d3mu_l*rold_r*un_r
406 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*t_x
407 : br4=br4+d3mu_l*rold_r*ut_r
408 : jl(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
409 : jl(3,2)=jl(3,2)+t_y*damu_l*br2+am*t_y*br4
410 : c-------------------
411 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*n_y
412 : br3=br3+d3mv_l*rold_r*un_r
413 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*t_y
414 : br4=br4+d3mv_l*rold_r*ut_r
415 : jl(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
416 : jl(3,3)=jl(3,3)+t_y*damv_l*br2+am*t_y*br4
417 : c-------------------
418 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
419 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
420 : jl(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
421 : jl(3,4)=jl(3,4)+t_y*damp_l*br2+am*t_y*br4
422 : c-------------------------------------------------------------
423 : c ------ f44444444444444444444444444444 ---------
424 : c-------------------------------------------------------------
425 : hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1
426 : hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1
427 : f(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
428 : c---------------------
429 : br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
430 : br1=br1+d3mr_l*hr_r
431 : jl(4,1)=damr_l*f(4)/am+am*br1
432 : c---------------------
433 : br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
434 : br1=br1+d3mu_l*hr_r
435 : jl(4,2)=damu_l*f(4)/am+am*br1
436 : c---------------------
437 : br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
438 : br1=br1+d3mv_l*hr_r
439 : jl(4,3)=damv_l*f(4)/am+am*br1
440 : c---------------------
441 : br1=d2mp_l*hr_l+mpl_m*ga/gm1
442 : br1=br1+d3mp_l*hr_r
443 : jl(4,4)=damp_l*f(4)/am+am*br1
444 : c----------------------------------------------------------
445 : c----------------------------------------------------------
446 : do 1 i=1,4
447 : dfro(i)=jl(1,i)
448 : dfrun(i)=jl(2,i)*n_x+jl(3,i)*n_y
449 : dfrut(i)=jl(2,i)*t_x+jl(3,i)*t_y
450 : dfret(i)=jl(4,i)
451 : 1 continue
452 : c----------------------------------------------------------------------
453 : return
454 : end
455 :
456 :
457 :
458 :
459 :
460 :
461 :
462 :
463 :
464 :
465 :
466 :
467 :
468 :
469 :
470 :
471 :
472 :
473 :
© Cast3M 2003 - All rights reserved.
Disclaimer