1 : C FAUUPT SOURCE BECC 10/06/11 21:15:06 6690
2 : SUBROUTINE FAUUPT(nesp,nsca,idim,
3 : & gaml,rhol,pl,unl,utl,uvl,etherl,
4 : & gamr,rhor,pr,unr,utr,uvr,etherr,
5 : & yl,yr,sl,sr,v_inf,f,
6 : & cellt)
7 : C************************************************************************
8 : C
9 : C PROJET : CASTEM 2000
10 : C
11 : C NOM : FAUUPT ('FLUX for Bas Mach' and thermally perfect)
12 : C
13 : C DESCRIPTION :
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, nsca, idim,iesp
47 : real*8 gaml, rhol, pl, unl, utl, uvl, etherl
48 : & , gamr, rhor, pr, unr, utr, uvr, etherr
49 : & , yl(*), yr(*), v_inf, f(*)
50 : & , cellt, sl(*), sr(*)
51 : C Some functions and parameters
52 : real*8 fau_m4, fau_p5, kp, ku, sigma
53 : C
54 : real*8 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 : c Computation of the specific total enthalpies.
63 : c------------------------------------------------------------------
64 : htotl = ((unl*unl) + (utl*utl) + (uvl*uvl)) / 2.0d0
65 : htotl = htotl + etherl
66 : htotl = htotl + (pl / rhol)
67 : htotr = ((unr*unr) + (utr*utr) + (uvr*uvr)) / 2.0d0
68 : htotr = htotr + etherr
69 : htotr = htotr + (pr / rhor)
70 : * write(*,*) htotl, htotr
71 : c------------------------------------------------------------------
72 : c Computing reference velocity
73 : c see Eq.(70) of the Ref.3).
74 : c------------------------------------------------------------------
75 : asonl=sqrt(gaml*pl/rhol)
76 : asonr=sqrt(gamr*pr/rhor)
77 : qq=0.5D0*sqrt((unl*unl)+(utl*utl)+(uvl*uvl))
78 : c------------------------------------------------------------------
79 : if(qq .lt. v_inf) then
80 : uco_l = v_inf
81 : else
82 : uco_l = qq
83 : endif
84 : c-----------------------------
85 : if(uco_l .ge. asonl) then
86 : uco_l = asonl
87 : endif
88 : c-------------------------------------------------------------------
89 : qq=0.5D0*sqrt((unr*unr)+(utr*utr)+(uvr*uvr))
90 : if(qq .lt. v_inf) then
91 : uco_r = v_inf
92 : else
93 : uco_r = qq
94 : endif
95 : c-------------------------
96 : if(uco_r .ge. asonr) then
97 : uco_r = asonr
98 : endif
99 : c-------------------------------------------------------------------
100 : c mco0 is the equivalent of M_{\inf} in Ref 3.
101 : c-------------------------------------------------------------------
102 : mco0 = max((uco_r / asonr), (uco_l / asonl))
103 : c write(*,*) 'mco0 = ', mco0
104 : c-------------------------------------------------------------------
105 : c Computation of the speed of sound;
106 : c numerical speed of sound at the interface is taken as an average
107 : c of the speeds of sounds of the neighbouring cells
108 : c-------------------------------------------------------------------
109 : asonm=0.5d0*(asonl + asonr)
110 : c-------------------------------------------------------------------
111 : c Computing numerical Mach number
112 : c ml, mr, mbar2 (used in the term for Mu, which does not involves
113 : c the cut-off speed, see equation (70), ref 3.)
114 : c Then, even if mc0 = 1, the term involving kp acts!
115 : c-------------------------------------------------------------------
116 : ml = unl / asonm
117 : mr = unr / asonm
118 : mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm)
119 : * write(*,*) asonm
120 : * write(*,*) ml, mr, mbar2
121 : c-------------------------------
122 : c-------------------------------------------------------------------
123 : c Scaling function for the speed of sound
124 : c We redefine mzero (see (71), ref 3).
125 : c-------------------------------------------------------------------
126 : C
127 : mzero = max((mco0 * mco0), mbar2)
128 : mzero = min(1.0D0, mzero)
129 : mzero = mzero ** 0.5D0
130 : C
131 : fson = mzero * (2.0D0 - mzero)
132 : alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson))
133 : C write(*,*) 'alpha = ', alpha
134 : C-------------------------------------------------------------------
135 : C mhalf (see (73), ref 3)
136 : C First, contribution of Mu; then contribution of ml and mr.
137 : C-------------------------------------------------------------------
138 : C
139 : rhom = 0.5D0 * (rhol + rhor)
140 : mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0)
141 : mhalf = -1.0D0 * mhalf * kp * (pr - pl) /
142 : & (fson * rhom * asonm * asonm)
143 : C
144 : mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0)
145 : c write(*,*) 'mhalf = ', mhalf
146 : uhalf = mhalf * asonm
147 : c write(*,*) 'unl, unr ' , unl, unr
148 : c write(*,*) 'mhalf, uhalf ', mhalf, uhalf
149 : C-------------------------------------------------------------------
150 : C phalf
151 : C-------------------------------------------------------------------
152 : fp5p = fau_p5(ml,1.0D0,alpha)
153 : fp5m = fau_p5(mr,-1.0D0,alpha)
154 : phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm *
155 : & (unr - unl)
156 : phalf = phalf + (fp5p * pl) + (fp5m * pr)
157 : c write(*,*) 'phalf = ', phalf
158 : c---------------------------------------------------------------------
159 : c Computing numerical fluxes
160 : c---------------------------------------------------------------------
161 : if(uhalf .ge. 0.0d0) then
162 : f(1) = uhalf * rhol
163 : f(2) = (f(1) * unl) + phalf
164 : f(3) = f(1) * utl
165 : if (idim .eq. 3) then
166 : f(4) = f(1) * uvl
167 : endif
168 : f(2 + idim) = f(1) * htotl
169 : do iesp = 1, nesp, 1
170 : f(2 + idim + iesp) = f(1) * yl(iesp)
171 : enddo
172 : do iesp = 1, nsca, 1
173 : f(2 + idim + nesp + iesp) = f(1) * sl(iesp)
174 : enddo
175 : else
176 : f(1) = uhalf * rhor
177 : f(2) = (f(1) * unr) + phalf
178 : f(3) = f(1) * utr
179 : if (idim .eq. 3) then
180 : f(4) = f(1) * uvr
181 : endif
182 : f(idim + 2) = f(1) * htotr
183 : do iesp = 1, nesp, 1
184 : f(2 + idim + iesp) = f(1) * yr(iesp)
185 : enddo
186 : do iesp = 1, nsca, 1
187 : f(2 + idim + nesp + iesp) = f(1) * sr(iesp)
188 : enddo
189 : endif
190 : c----------------------------------------------------------------------
191 : cellt = max(abs(unl),abs(unr)) + asonm
192 : c write(*,*) 'cellt = ', cellt
193 : c----------------------------------------------------------------------
194 : return
195 : end
196 :
© Cast3M 2003 - All rights reserved.
Disclaimer