- Email: [email protected]
An efficient flux difference splitting algorithm for unsteady duct flows of a real gas
APPLIED NUMERICAL MATHEMATICS ELSEVIER
Applied
Numerical
Mathematics
1.5 (1994) 27-52
An efficient flux difference splitting algorithm for unsteady duct flows of a real gas P. Glaister Department of Mathematics, P.O. Box 220, University of Reading, Reading, RG6 2Ax, UK
Abstract A numerical scheme is presented for the solution of the Euler equations of compressible flow of a real gas in a single spatial coordinate. This includes flow in a duct of variable cross-section, as well as flow with slab, cylindrical or spherical symmetry, as well as the case of an ideal gas, and can be useful when testing codes for the two-dimensional equations governing compressible flow of a real gas. The resulting scheme requires an average of the flow variables across the interface between cells, and this average is chosen to be the arithmetic mean for computational efficiency, which is in contrast to the usual “square root” averages found in this type of scheme. The scheme is applied with success to five problems with either slab or cylindrical symmetry and for a number of equations of state. The results compare favourably with the results from other schemes.
1. Introduction Recently [3], an approximate Riemann solver was developed for the shallow water which used the arithmetic mean for the average of the flow variables across a cell The purpose of this paper is to present a similar type of Riemann solver for equations for a real gas in a duct of variable cross-section and using the arithmetic computational efficiency. As well as the usual case of slab symmetry, this includes cylindrical or spherical symmetry, as well as the case of an ideal gas. The scheme is five problems with interacting and/or reflecting shocks in slab or cylindrical geometry number of equations of state.
equations interface. the Euler mean for flows with applied to and for a
2. Governing equations The “one-dimensional” Euler written in conservation form as
equations
governing
=0, Pf+ L(w)Pu)r S(r) 0168-9274/94/$07.00
0 1994 Elsevier
SSDI0168-9274(94)00027-E
Science
the flow of a compressible
gas can be
(2.la)
B.V. All rights reserved
28
P. Glaister /Applied
Numerical Mathematics
15 (1994) 27-52
=-p,< (4 + +(r)puZ), s(y)
ef +
L(G-)u(e S(Y)
+p)),
= 0,
(2.lb)
(2.lc)
where the total energy e = pi + $M.L~.
(2.ld)
The quantities p, U, i, and p represent density, velocity, specific internal energy, and pressure, respectively, at a general position Y along the duct (of cross-section S(r)) and at time t. The particular example S(r) = rN represents slab symmetry (N = O), cylindrical symmetry (N = l), or spherical symmetry (N = 2). The coordinate Y is given by Y=x, \lx2=, /Z$?LP+ in the case of central symmetry with N = 0, 1, 2, respectively, where X, y, and z represent Cartesian coordinates. Finally, there is an equation of state relating p, p, and i of the form P =P(P,
(2.le)
i>
together with first partial derivatives. In the case of an ideal gas, (2.lf)
P = (Y - l)Pi,
where y is the ratio of specific heat capacities of the fluid. In the next section we find approximate solutions of Eqs. (2.la)-(2.ld); them in standard “conservation form” as w, + (F(w)),
but first we rewrite (2.2a)
=f(w),
where (2.2b)
w = (P, PU, e)T, P +pu2,
F(w) = (~4
together with (2.ldM2.le) f(w) = -
S’(r)
-( S(r)
u(e +p))T7
(2.2c)
and
PU, PU*,
u(e +I+’
(2.2d)
is a source term arising from the geometry. The special case S(r) = constant is considered first, and the extension to the general case is developed from the special case.
3. Approximate
Riemann
solver
We consider the approximate solution w,! = w(rj, t,) to consist of a set of piecewise constants, and solve approximately the associated Riemann problems at the interface separating adjacent states. An approximate Jacobian needs to be constructed across an interface so
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
that shock capturing is automatic, and this represents evaluated either side of the interface. Consider firstly Eqs. (2.2a) with S(r) = constant,
29
an average of the Jacobian
matrix
3.1. Structure The Jacobian matrix 0 -u2__
aF
A=
Pi
P
_
I
+i-y
P i P
aw
U2
P
I
-+i+T ;;; P
0
1
Plf
‘Pi 2u -P
u2 )
+i-y
-
Pi -
P -+i+ P
U2
--2
P u2Pi P
u+-
‘Pi P
(3.1) of the flux function given by
F(w)
has eigenvalues
Ai with corresponding
1, u +a,
e, =
i
A2=u-a,
e2=
l,u-a,
‘+i+G+uu
e3=
I
I
-+i+f-ua
u2
T ,
P
i
h3=u,
(3.2a)
P
P
u2
l,u,i+----p?
2
ei, i = 1, 2, 3,
7
2 h1=u+a,
eigenvectors
!
pT ,
(3.2~)
Pi I
where the sound speed, a, is given by a2=Pp
+PPi/P2Y
(3.3)
using the derivatives of the equation of state (2.le). 3.2. Shock capturing
Consider two adjacent states wL and wn (left and right) given at either end of the cell (_YL,rR., and consider also the algebraic problem of finding an approximate Jacobian A = A&, wn) in this cell such that AAw = AF,
(3.4)
where A(*) = (.I, - (-IL, w = (p, pu, eIT, and F = (pu, p +pu2, u(e +p)jT. A solution to this problem, for arbitrary jumps Aw, can be used to obtain a conservative scheme with good shock-capturing properties. This technique was originally introduced by Roe [5] for the Euler equations with slab symmetry and for ideal gases.
P. Glaister/Applied NumericalMathematics15 (1994) 27-52
30
3.3. Construction
of A’
The solution
of Eq. (3.4) is straightforward
to determine
and is given by 0
-Pi P
A=
where
the averages u=
P=
(UL+ UR),
fi = /K,
i(PL +PRh
i = +(iL
u2= f(z$ + ug,
P-6)
+ iR),
(3.7a)
p=
(3.7b)
and
(see Appendix A). The quantities equation of state satisfying Ap =fi,Ap (see Appendix
I
q,
A). Th e specific
P,(P> q +
averages
IAPI/F
IApl/p+
we propose
6~
[Ail/,!
where 6p
= Ap
-p,(p,
if I Ap l/p < lo-“,
(3.10a)
otherwise,
(3.10b)
Ap’
E)Ap
(3.1Oc)
6p
I Ap I/jj + I Ai I/i
m is machine-dependent
of the
are
if I Ai I /i G lo-“, /Ail/i
Pi(P> ‘) +
of the derivatives
(3.9)
(Pi(P> ‘)T Pi =
averages
+fiiAi
P,(P,
sp =
cP and ci are appropriate
Ai ’
othemise7
(3.10d)
and -p,(p,
E)Ai.
(3.11)
It is a simple matter to check that the averages given above satisfy (3.9). We note that alternative averages have been suggested for the solution of (3.9); see for example [2,7]. In Section 3.7 we discuss the specific case of an ideal gas giving rise to simplified averages.
P. Glaister /Applied
Numerical Mathematics
15 (1994) 27-52
31
3.4. Approximate eigenualues and eigenvectors
Now, the important quantities that are needed for the scheme are the eigenvalues eigenvectors e’, of A, and it is a simple matter to show that these are given by
h,=u+a’, e’1,2
=
Ai and
h3= ii,
(3.12)
(1, li+6,~/P+I+fE;I+aapAi/p~u~)T,
(3.13a)
/i2=ii4,
t?, = (1, U, i+ ;E;z-pfi,/fii-
fop/&)‘,
(3.13b)
where a’ = (2 + am
+ aApAifii/pZ)1’2.
(We have used the identities these.)
(3.14)
U2 _;2J_-2-
u
-
1 ,(Au12 and pz - ,Gi= $ApAi in determining
3.5. Projection It is necessary to project a general jump Aw onto the eigenvectors Aw = t
gi as
(YiZi
(3.15)
i=l
and by virtue of (3.4) we then have (3.16)
AF = 5 ,&?,. i=l
Solving (3.15) gives (At, f +~AU + +(~~)‘a~) aI,2
&;,=
(3.17a)
>
=
A
p
%f
;&2A~ 4G2
(3.17b)
.
The simplified results in the case of an ideal gas are given in Section 3.7. 3.6. Numerical scheme Thus in (2.2a), with S(r) = constant, we can approximate
F,-
g
_Aw =k j=l Ar Ar
=A_
=
i=l
Ar
’
(3.18)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
32
by virtue of the analysis of this section. Applying upwind differencing to equation (2.2a) using the approximation in (3.18) then gives ,,v+’ - wi”
(3.19)
At
where Ar and At represent the mesh spacings in the Y and t directions, the subscript j - l/2 refers to the cell [rj_i, rj] and (3.20) represent the positive and negative parts of ii. This gives the following first-order algorithm for the solution of (2.1): At_ _ add - zAi~iZi
to wR when hi > 0
(3.21a)
At _ _ add - zAiuiZi
to wL when hi < 0.
(3.21b)
or
Hence the only quantities required for the algorithm are ii, Ci, and Gi, and only one square root is taken in each computational cell, namely that in Eq. (3.15). Thus, we note the direction of flow of information given by the approximate eigenvalues Ai and use this information to update the solution consistent with the theory of characteristics of Eq. (2.2a) with S(r) = constant. In addition, second-order transfers of these first-order increments can be made to achieve higher accuracy, providing they are limited to maintain monotonicity. The use of these “flux-limiters” improves accuracy without introducing nonphysical spurious oscillations, especially at shocks. To allow rarefaction waves to be treated correctly, and hence to avoid entropy-violating solutions, the first-order increment can be considered as two separate increments being sent to either end of the cell. Specifically, the modified version of the scheme can be written cell-wise as wj”-i’ = wi”_
At
1
+
-6iL&idi
Ar At w,!+i = wi”+ -f$Giei, Ar
where
)
i = 1, 2, 3,
(3.22a)
i = 1, 2, 3,
(3.22b)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
33
(3.23~) (3.23d)
The expressions $ and @ are given by j_,,/,hi),
(3.24a)
@ = maX(j_1/*hir j+1/2ii),
(3.24b)
if = min(j_,phi,
where j_l,2ii = hi, and th ese will be different for a rarefaction wave. Finally, in the case of a duct of variable cross-section, or for cylindrically or spherically symmetric flows, there is a source term present:
S’(r) j-= - -( S(r)
PU,
pu2,u(e +P))~.
Upwinding the source term as described in [3] enables Eq. (2.2a) to be solved approximately the case S(r) varying using appropriately modified wavestrengths Gi.
in
3.7. Ideal gas case
In this case, the equation
of state (2.le)
is given by (2.lf), so that p, = (y - 1)i and
pi = (y - 1)~. Thus (3.11) becomes
6p = (y - l)A(pi)
- (y - 1)iAp - (y - 1)pAi = 0
after some simplification. fip =
(y - l)E,
Thus the expressions in (3.10a)-(3.10d)
(3.25) give, in all cases,
fii = (y - l)P.
Therefore, substituting (3.8) into (3.14) and using these averages gives, after simplification use of $= $ + +ApAi,
(3.26) and
(3.27) Similarly, the expression Is/P + i + aApAi/p- in (3.13a) becomes >/(y - 1) where ;;?= yp/p, and t - j?j?,/rYi in (3.13b) is zero, leaving $L2 - $(Au)~/(~ - 1). Finally, in (3.17b), l__=-(Au)’ 4a’*
YPh
G2
=-
aL a’2 a
34
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
4. Test problems and numerical results We consider five test problems for various equations of state to assess the scheme in Section 3. Problem 1. This is the well-known shock tube problem of Sod [6] for the Euler equations with
y = 1.4 and initial data
PY U7 p = i
1, 0, 1,
r< +,
0.125, 0, 0.1,
r> $.
The main features of the exact solution are a shock moving to the right followed by a contact discontinuity, also moving to the right, but more slowly, and an expansion fan moving to the left.
u 0.99..
0.M
..
0.33 . .
0.3 0.9
0.6
0.P
X
.
4.33
3.46
0.w
0.3
0.6
0.9
X
0.3
0.6
0.9
X
.
4.63
..
-0.w..
I 2.61..
1.73
..
0.w
-0.36
an
-I
.oo
-0.119
+
I
-1.73
-2.61..
Fig. 1. Solution
of Problem
1 at t = 0.144 s.
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
I
*: “’
““..__
..‘\
x
s
.,...
..z.
“1
9 0
a
35
36
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
x
e x 9 D
;
a 3
D
”
yl
-
c x
9
.,.,_:... e.
.D
”
----.
L..
x
P
x 9 D
x 9 -
:, j..
.__,
I
P. Glaister/Applied
Numerical Mathematics 15 (1994) 27-52
37
Fig. 1 shows the approximate and exact solutions for p, u, p, and e at t = 0.144 s using 100 mesh points, and a good agreement between the two can be seen. Problem 2. This problem is usually described as “Two Interacting
Blast Waves” and is a shock tube problem for the Euler equations in slab symmetry with an ideal gas with y = 1.4 and initial data P, u, P =
1, 0, 1000, 1, 0, 0.1, 1, 0, 100,
0
The walls of the tube at r = 0 and Y= 1 are assumed to be perfectly reflecting. Two strong blast waves develop and collide producing a complex flow. A detailed description of the time evolution of the flow can be found in [8]. Figs. 2 and 3 show the approximate solution for p and u at t = 0.01, 0.016, 0.026, 0.028, 0.03, 0.032, 0.034 and 0.038 seconds using 400 mesh points, and the results compare well with those given by Woodward and Collella [S] using a more complicated Riemann solver.
JTH SIAB SYnnETRY - Shock Reflection PI
” 0.w
KEY 0.44
p 0.33
”
-
P-
Denslry
Velocity Pl-esEU-e
I - Internal energy
-
Exact solution Approxlmete solutlor
equat,on of state for an Ideal c~as 50 fleshpotnts 89 Time steps Ax = 0.02 At = 0.0065 Pressure rat,0 = 10 'Superbee' llmlter used
Reflected Boundary Condrtlons etr=O
Fig. 4. Solution
of Problem
3: ideal equation
of state.
38
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
Problem 3. This problem is concerned with shock reflection in one dimension of a gas governed by the Euler equations in slab symmetry. We consider a region 0 < Y< 1 divided into fifty equally spaced mesh points and the initial conditions are (p, U, i, p> = (po, uO, i,, p(pO, i,)). This represents a gas of constant density and pressure moving towards Y= 0. The boundary at Y= 0 is a rigid wall and the exact solution describes shock reflection from the wall. Five equations of state are chosen:
(4 ideal gas P = (y -
W,
where we take the ratio of specific heat capacities monatomic gas; (b) stiffened gas (sometimes used for a liquid) P =B(P/P,
- 1) +
(Y -
to be 5/3 corresponding
y
to a
l)pi,
where we take B = 1 and y = 5/3;
I& SYHHFIBT " 0.w
SDLUTION OF THE Fill FR EOUATIONS I
P
- Shock
Reflection
0.66 0.33
p -
Denslry
”
Velocity
-
PreSSWe
P-
I -
q3
0.6
0.P
Internal
energy
x Exact
-
-03..
soLut!on
~OOQIO Approximate
eolutlor
au... av9..
et~ffened
equatlon
of stare 50 flesh po,nts 64 Time
I ’ b
P ?.P..
I
steps
Ax = 0.02 At = 0.0046 Pressure
,.u ..
‘Superbee’
ratlo llmtter
=
IO used
0.1...
0.5
0.b
0.9
x
-0.7. -I.&n Reflected
Boundary atx=O
-2.a..
i_ at
time
t
=
Fig. 5. Solution
0.292
of Problem
3: stiffened
equation
of state.
Condltlons
39
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
(c) copper (due to R.K. Osborne, see [4]) p=
(& IM
a,+%9
+Jqb,+5(b,
+U)
+Jqco+c,6))],
where E =p,i and 5 = p/p,, - 1 and the constants a,, a2, b,, b,, b,, cO, cl, and 4. can be found in [2]; (4 covolume (used in connection with combustion environments) p = (7 - l)Pi 1-pb ’ where we take b = 0.8 and y = 1.4; (e) two molecular vibrating gas p + (Y - GLp = (y - qp’ 1, eP”/P - 1 where we take y = 1.4 and p = 0.5902 which corresponds
ON P
OF
THF
FIIIw
I
I
LAB " 0.w
- Shock
to molecular oxygen at 2000 K.
Reflect,on
KEI 0.66
p 0.33
”
-
P0.3
0.6
0.9
x
(13
0.6
0.9
VetocIty Pr08sUfW
I -
I
oens1ty
Internal
energy
x -
4.33..
Exact
solutron
oooooo Approxlmare
solurlor
-0.M.. 4.w.-
equat,on
of etate
for copper
-
I . 0.72 .-
P 20. 1I.J
O.&Y
6.7
0.24 ..
-
50 Neeh
points
54 Time
steps
Ax = 0.02 At = 0.0053 Pressure 'Superbee'
rat,0
= IO
L,n,ter
used
i 0.1
d.7
4.11
-I,.,
-0.40
-2o.I
4.72
0.6
0.9
x
Reflected
etx=O
-
at t,me t =
Fig. 6. Solution
u.
aI>
of Problem
Boundary
3: equation
of state for copper.
Condlt,ons
40
P. Glaister /Applied
Numerical Mathematics
T
I5 (1994) 27-52
- Shock Reflection
p -
Denslry
Velocity p - Preesure I - Interna 1 energy u -
Exact solut,on ooooo~Approxlmete solutlor
covolume equat Ion of state 50 Mesh pointa 95 Time steps Ax = 0.02 At = 0.0066 Pressure rat10 = IO 'Superbee' llmlter used
Reflected Boundary Condltlons atx=O et time t =
Fig. 7. Solution
0.628
of Problem
3: covolume
equation
of state.
For each equation of state we take u0 = - 1, and take p0 = 8.9 for (c), p0 = 1 for (b), and p0 = 0.1 for (a), (d), and (e). The initial condition for i, (in each case) is chosen to correspond to a shock strength p+/po = 10 where p+ denotes the pressure behind the shock and p. =p(p,,, io) denotes the pressure ahead of the shock. The results for these cases are given in Figs. 4-8, together with the exact solution when the shock has moved a distance of 0.3. Problem 4. This problem is concerned with a converging cylindrical shock and we consider a region 0 < r < 200 for the cylindrically symmetric case given by the Euler equations (2.2a) with S(r) = Y. Initially, a cylindrical diaphragm of radius r = 100 separates two uniform regions of gas at rest. The initial conditions are p = 4, p = 4 in the outer region, and p = 1, p = 1 in the inner region. When the diaphragm is removed at t = 0, a converging shock wave followed by a converging contact discontinuity moves towards the axis, Y = 0, and a diverging rarefaction wave moves outwards. The shock accelerates as it approaches the axis of symmetry and is reflected from the axis, interacts with the contact discontinuity (still converging), which results in a transmitted shock, a converging contact discontinuity and a weak converging reflected shock. Four equations of state are considered, namely (a), (b), Cd), and (e) associated with Problem 3 above.
P. Glaister /Applied
N OF
THF
FillFR
Numerical
1 Ai SY" 0.w
FWAI_IQM
‘ir.
Mathematics
- Shock
15 (1994) 27-52
41
Reflectron
KEY
p-
Density
u -
Veloclry
p - Pressure 1 - Internal -
Exact
energy
solution
ooo~x, Approximate
eq”atl0”
eolutlor
of *fete
for oxygen 50 Mesh
polnrs
142 TIM
steps
Ax = 0.02 At = 0.0068 Pressure 'Superbee'
ret10
= IO
llmlrer
used
0 Reflected
Boundary
Condltlons
atx=O
at
t,me
t =
Fig. 8. Solution
0.769
of Problem
3: equation
of state for oxygen.
Figs. 9-12 refer to the four equations of state using 200 mesh points. Three output times have been featured. For each case, we observe the development of the shock and the contact discontinuity, the reflection of the shock, and the interaction of the shock with the contact discontinuity. The results for equation of state (a) compare favourably with results produced using a more complex algorithm [l]. Problem 5. The final problem is concerned with a converging and a diverging cylindrical shock in 0 < Y G 200. Initially, two cylindrical diaphragms of radii r = 50 and Y = 150 separate three uniform regions of gas at rest. The initial conditions are p = 4, p = 4 in the inner and outer regions, and p = 1, p = 1 in the middle region. When the diaphragms are removed at t = 0, a converging shock wave, followed by a converging contact discontinuity moves towards the axis, r = 0, and a diverging shock wave, followed by a diverging contact discontinuity moves away from the axis. The shocks subsequently interact, resulting in a diverging shock wave weakening in strength, together with a converging shock wave increasing in strength. Each of these shocks then interact with the corresponding contact discontinuity, as in Problem 4, and this results, for each interaction, in a transmitted shock, a weak reflected shock and a contact discontinuity. Four equations of state are again considered, namely (a), (b), (d), and (e).
42
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
~JL. .-..... ..... ...... 54
8.A) E
“‘\
P. Glaister/Applied
Numerical
Mathematics
15 (1994) 27-52
c:, \
8 \
0 \P-’
!
a
I
2
43
.
5
L
I
1
‘.
!
0
“*.‘i
JI
A
N
d
d
d.
.M L
0 8
8
.
E
.
44
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
:
9
i
i
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
45
L
iL..
’ . . .._ -....
s
..
*
8
.
\
”
9
1 R +
Fi +
46
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
J---+.... -a....
.. .
s E
\
s.c .
\ .-........ ..
.‘....._
s
0
.. E .
. \
:
*
.. % / .:.... . .. .._
. . . . . . -.....
3.3
If!
sl d
i
.. s .. XJ
PI
9
+
r
ci
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
J
E d
t d
i% d
L
a-.._.... ...
.-.
s 4
.
\ .
8 *\
L ..
w ..
-a.......* .’
s_
.. 51; . . .
. ‘.
/.
,
f
. . . . .
b
\
8
. . . . . . . . .. 0
2
‘y 9
i
i
48
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
..
-8”*...\.. s
..I %. S’ s i*......_-..., ..... -,,~
L
1*
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
49
L
I s! *: :.
..I .
\ l-L_ 8
/
s
L;
it
n 2
“! 7
.
‘K. *-*..._ -..-.....* i
1
97
d +
i
9
I
1 I
iiN..
-..
s E
.\’ % / . ..::.* -........ -....Q.....
.....
a d
2
fj
I
.
.
.
/ / ,:.... ... ..-..... *.....
s
-.....
-.
\
*
P. Glaister/Applied
50
Numerical Mathematics 15 (1994) 27-52
Figs. 13-16 refer to the four equations of state using 200 mesh points. Three output times have again been featured. For each case, we observe the development of the shocks and contact discontinuities, the interaction of shocks, and the subsequent interaction of each shock with a contact discontinuity, resulting in transmitted and reflected shocks and contact discontinuities.
5. Conclusions
We have presented an efficient shock-capturing algorithm for the Euler equations for real gases in a duct of variable cross-section. We have utilised the arithmetic mean for the averages of the flow variables in computational cells. Numerical experiments typically show a 10% decrease in the computer time required when the square root averages are replaced by the arithmetic mean in the case of an ideal gas. The main features of this scheme are (i) it can be used for slab, cylindrically or spherically symmetric problems with confidence and (ii) it can be used as a comparison with results from two-dimensional schemes by choosing a large number of mesh points for accuracy and not be expensive on computing. The scheme can also be used for general duct flows.
Appendix A
We derive here_the matrix k in Eq. (3.5). To determine A we first write Aw and AF in terms of AU, where w = (P, PU, e)T,
F= (PU, P +Pu2,
u(e +P))~,
24= (P, u, P)‘.
Following the identities AP = PP,
(A4
A
(f-w W.3) W)
= PAu + UAp,
A(pi) = PAi + iAp, A( pu2) = ZAp + pAu2 = ZAp + 2ptiAu, where P = 3PL + PR)? u = gu,
+ UR),
E= +(iL + iR),
the arithmetic mean of left and right states, 2 = ;(ut
+ u;>
(A-8)
the arithmetic mean of the square of the velocity, we can write Aw=tiAu.
(A-9)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
where 1 ii
i=
j+ ;z I
0
P
0.
pu
#z
Similarly, Ap =j?,Ap
1
0
51
(A.lO)
+fiiAi,
(A.ll)
where fiP and ci are appropriate averages of the derivatives satisfy this equation exactly (see Section 3.31,
of the equation
of state chosen
to
A(up) = UAp +pAu,
(A.12)
P=;(PL+PR)
(A.13)
where
and A(pui)
= piAu + UA(pi)
-= piAu + updi
+ EEAp,
(A.14)
where z = +(pLiL + pRiR).
(A.15)
Finally, A( pu3) = u3Ap + ,5Au3 = u3Ap + 3,$Au,
(A.16)
z = +(u;
(A.17)
where + u;),
and ‘?_ 1 u - &L;+ULUR+U;)
(A.18)
is an average of the square of the velocity. (N.B. the choice in (A.14) is made so that the eigenvalues of the approximate Jacobian k have the simplest possible form. In particular, k has fi as one eigenvalue.) Combining (A.21, (A.4), (A.111, (A.121, (A.14) and (A.16) gives AF=eAu,
(A.19)
where
c=
u
P
0
fiP + 2
2pu
fii
ii& + uE + $3
p+ ;pu2+ ! and thus from (A.91 and (A.191 AF = &-‘Arv.
i#+pu
1
(A.20)
(A.21)
P. Glaister/Applied Numerical Mathematics15 (1994) 27-52
52
Therefore a solution of (3.4) is obtained with the approximate the identities
Jacobian A’ = &‘,
and using (A.22)
2U”+z=3>,
(A.23)
2ii3 + z = 3&P,
(A.24)
where ii = ULUR
(A.29
is the geometric mean of left and right states, gives rise to the matrix in (3.5), where (A.26)
References [II
M. Ben-Artzi and J. Falcovitz, A second order Godunov-type scheme for compressible fluid dynamics, .I Comput. Phys. 55 (1984) l-32. Dl P. Glaister, An approximate linearised Riemann solver for real gases, J. Comput. Phys. 74 (1988) 382-408. open channel flows, Appl. Numer. [31 P. Glaister, An efficient numerical method for subcritical and supercritical Math. 11 (1993) 497-508. [41 T.D. Riney, Numerical Evaluation of Hypervelocity Impact Phenomena (Academic Press, New York, 1970). Riemann solvers, parameter vectors and difference schemes, .I. Comput. Phys. 43 (1981) [51 P.L. Roe, Approximate 357-372. laws, .I. 161 G.A. Sod, A. survey of several finite difference methods for systems of nonlinear hyperbolic conservation Comput. Phys. 27 (1978) 1-31. flux-vector splitting and Roe average for an equilibrium real gas, J. Comput. Phys. 89 [71 M. Vinokur, Generalised (1990) 276-300. fluid flow with strong shocks, J. 181 P.R. Woodward and P. Colella, The numerical simulation of two-dimensional Comput. Phys. 54 (1984) 115-173.
Applied
Numerical
Mathematics
1.5 (1994) 27-52
An efficient flux difference splitting algorithm for unsteady duct flows of a real gas P. Glaister Department of Mathematics, P.O. Box 220, University of Reading, Reading, RG6 2Ax, UK
Abstract A numerical scheme is presented for the solution of the Euler equations of compressible flow of a real gas in a single spatial coordinate. This includes flow in a duct of variable cross-section, as well as flow with slab, cylindrical or spherical symmetry, as well as the case of an ideal gas, and can be useful when testing codes for the two-dimensional equations governing compressible flow of a real gas. The resulting scheme requires an average of the flow variables across the interface between cells, and this average is chosen to be the arithmetic mean for computational efficiency, which is in contrast to the usual “square root” averages found in this type of scheme. The scheme is applied with success to five problems with either slab or cylindrical symmetry and for a number of equations of state. The results compare favourably with the results from other schemes.
1. Introduction Recently [3], an approximate Riemann solver was developed for the shallow water which used the arithmetic mean for the average of the flow variables across a cell The purpose of this paper is to present a similar type of Riemann solver for equations for a real gas in a duct of variable cross-section and using the arithmetic computational efficiency. As well as the usual case of slab symmetry, this includes cylindrical or spherical symmetry, as well as the case of an ideal gas. The scheme is five problems with interacting and/or reflecting shocks in slab or cylindrical geometry number of equations of state.
equations interface. the Euler mean for flows with applied to and for a
2. Governing equations The “one-dimensional” Euler written in conservation form as
equations
governing
=0, Pf+ L(w)Pu)r S(r) 0168-9274/94/$07.00
0 1994 Elsevier
SSDI0168-9274(94)00027-E
Science
the flow of a compressible
gas can be
(2.la)
B.V. All rights reserved
28
P. Glaister /Applied
Numerical Mathematics
15 (1994) 27-52
=-p,< (4 + +(r)puZ), s(y)
ef +
L(G-)u(e S(Y)
+p)),
= 0,
(2.lb)
(2.lc)
where the total energy e = pi + $M.L~.
(2.ld)
The quantities p, U, i, and p represent density, velocity, specific internal energy, and pressure, respectively, at a general position Y along the duct (of cross-section S(r)) and at time t. The particular example S(r) = rN represents slab symmetry (N = O), cylindrical symmetry (N = l), or spherical symmetry (N = 2). The coordinate Y is given by Y=x, \lx2=, /Z$?LP+ in the case of central symmetry with N = 0, 1, 2, respectively, where X, y, and z represent Cartesian coordinates. Finally, there is an equation of state relating p, p, and i of the form P =P(P,
(2.le)
i>
together with first partial derivatives. In the case of an ideal gas, (2.lf)
P = (Y - l)Pi,
where y is the ratio of specific heat capacities of the fluid. In the next section we find approximate solutions of Eqs. (2.la)-(2.ld); them in standard “conservation form” as w, + (F(w)),
but first we rewrite (2.2a)
=f(w),
where (2.2b)
w = (P, PU, e)T, P +pu2,
F(w) = (~4
together with (2.ldM2.le) f(w) = -
S’(r)
-( S(r)
u(e +p))T7
(2.2c)
and
PU, PU*,
u(e +I+’
(2.2d)
is a source term arising from the geometry. The special case S(r) = constant is considered first, and the extension to the general case is developed from the special case.
3. Approximate
Riemann
solver
We consider the approximate solution w,! = w(rj, t,) to consist of a set of piecewise constants, and solve approximately the associated Riemann problems at the interface separating adjacent states. An approximate Jacobian needs to be constructed across an interface so
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
that shock capturing is automatic, and this represents evaluated either side of the interface. Consider firstly Eqs. (2.2a) with S(r) = constant,
29
an average of the Jacobian
matrix
3.1. Structure The Jacobian matrix 0 -u2__
aF
A=
Pi
P
_
I
+i-y
P i P
aw
U2
P
I
-+i+T ;;; P
0
1
Plf
‘Pi 2u -P
u2 )
+i-y
-
Pi -
P -+i+ P
U2
--2
P u2Pi P
u+-
‘Pi P
(3.1) of the flux function given by
F(w)
has eigenvalues
Ai with corresponding
1, u +a,
e, =
i
A2=u-a,
e2=
l,u-a,
‘+i+G+uu
e3=
I
I
-+i+f-ua
u2
T ,
P
i
h3=u,
(3.2a)
P
P
u2
l,u,i+----p?
2
ei, i = 1, 2, 3,
7
2 h1=u+a,
eigenvectors
!
pT ,
(3.2~)
Pi I
where the sound speed, a, is given by a2=Pp
+PPi/P2Y
(3.3)
using the derivatives of the equation of state (2.le). 3.2. Shock capturing
Consider two adjacent states wL and wn (left and right) given at either end of the cell (_YL,rR., and consider also the algebraic problem of finding an approximate Jacobian A = A&, wn) in this cell such that AAw = AF,
(3.4)
where A(*) = (.I, - (-IL, w = (p, pu, eIT, and F = (pu, p +pu2, u(e +p)jT. A solution to this problem, for arbitrary jumps Aw, can be used to obtain a conservative scheme with good shock-capturing properties. This technique was originally introduced by Roe [5] for the Euler equations with slab symmetry and for ideal gases.
P. Glaister/Applied NumericalMathematics15 (1994) 27-52
30
3.3. Construction
of A’
The solution
of Eq. (3.4) is straightforward
to determine
and is given by 0
-Pi P
A=
where
the averages u=
P=
(UL+ UR),
fi = /K,
i(PL +PRh
i = +(iL
u2= f(z$ + ug,
P-6)
+ iR),
(3.7a)
p=
(3.7b)
and
(see Appendix A). The quantities equation of state satisfying Ap =fi,Ap (see Appendix
I
q,
A). Th e specific
P,(P> q +
averages
IAPI/F
IApl/p+
we propose
6~
[Ail/,!
where 6p
= Ap
-p,(p,
if I Ap l/p < lo-“,
(3.10a)
otherwise,
(3.10b)
Ap’
E)Ap
(3.1Oc)
6p
I Ap I/jj + I Ai I/i
m is machine-dependent
of the
are
if I Ai I /i G lo-“, /Ail/i
Pi(P> ‘) +
of the derivatives
(3.9)
(Pi(P> ‘)T Pi =
averages
+fiiAi
P,(P,
sp =
cP and ci are appropriate
Ai ’
othemise7
(3.10d)
and -p,(p,
E)Ai.
(3.11)
It is a simple matter to check that the averages given above satisfy (3.9). We note that alternative averages have been suggested for the solution of (3.9); see for example [2,7]. In Section 3.7 we discuss the specific case of an ideal gas giving rise to simplified averages.
P. Glaister /Applied
Numerical Mathematics
15 (1994) 27-52
31
3.4. Approximate eigenualues and eigenvectors
Now, the important quantities that are needed for the scheme are the eigenvalues eigenvectors e’, of A, and it is a simple matter to show that these are given by
h,=u+a’, e’1,2
=
Ai and
h3= ii,
(3.12)
(1, li+6,~/P+I+fE;I+aapAi/p~u~)T,
(3.13a)
/i2=ii4,
t?, = (1, U, i+ ;E;z-pfi,/fii-
fop/&)‘,
(3.13b)
where a’ = (2 + am
+ aApAifii/pZ)1’2.
(We have used the identities these.)
(3.14)
U2 _;2J_-2-
u
-
1 ,(Au12 and pz - ,Gi= $ApAi in determining
3.5. Projection It is necessary to project a general jump Aw onto the eigenvectors Aw = t
gi as
(YiZi
(3.15)
i=l
and by virtue of (3.4) we then have (3.16)
AF = 5 ,&?,. i=l
Solving (3.15) gives (At, f +~AU + +(~~)‘a~) aI,2
&;,=
(3.17a)
>
=
A
p
%f
;&2A~ 4G2
(3.17b)
.
The simplified results in the case of an ideal gas are given in Section 3.7. 3.6. Numerical scheme Thus in (2.2a), with S(r) = constant, we can approximate
F,-
g
_Aw =k j=l Ar Ar
=A_
=
i=l
Ar
’
(3.18)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
32
by virtue of the analysis of this section. Applying upwind differencing to equation (2.2a) using the approximation in (3.18) then gives ,,v+’ - wi”
(3.19)
At
where Ar and At represent the mesh spacings in the Y and t directions, the subscript j - l/2 refers to the cell [rj_i, rj] and (3.20) represent the positive and negative parts of ii. This gives the following first-order algorithm for the solution of (2.1): At_ _ add - zAi~iZi
to wR when hi > 0
(3.21a)
At _ _ add - zAiuiZi
to wL when hi < 0.
(3.21b)
or
Hence the only quantities required for the algorithm are ii, Ci, and Gi, and only one square root is taken in each computational cell, namely that in Eq. (3.15). Thus, we note the direction of flow of information given by the approximate eigenvalues Ai and use this information to update the solution consistent with the theory of characteristics of Eq. (2.2a) with S(r) = constant. In addition, second-order transfers of these first-order increments can be made to achieve higher accuracy, providing they are limited to maintain monotonicity. The use of these “flux-limiters” improves accuracy without introducing nonphysical spurious oscillations, especially at shocks. To allow rarefaction waves to be treated correctly, and hence to avoid entropy-violating solutions, the first-order increment can be considered as two separate increments being sent to either end of the cell. Specifically, the modified version of the scheme can be written cell-wise as wj”-i’ = wi”_
At
1
+
-6iL&idi
Ar At w,!+i = wi”+ -f$Giei, Ar
where
)
i = 1, 2, 3,
(3.22a)
i = 1, 2, 3,
(3.22b)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
33
(3.23~) (3.23d)
The expressions $ and @ are given by j_,,/,hi),
(3.24a)
@ = maX(j_1/*hir j+1/2ii),
(3.24b)
if = min(j_,phi,
where j_l,2ii = hi, and th ese will be different for a rarefaction wave. Finally, in the case of a duct of variable cross-section, or for cylindrically or spherically symmetric flows, there is a source term present:
S’(r) j-= - -( S(r)
PU,
pu2,u(e +P))~.
Upwinding the source term as described in [3] enables Eq. (2.2a) to be solved approximately the case S(r) varying using appropriately modified wavestrengths Gi.
in
3.7. Ideal gas case
In this case, the equation
of state (2.le)
is given by (2.lf), so that p, = (y - 1)i and
pi = (y - 1)~. Thus (3.11) becomes
6p = (y - l)A(pi)
- (y - 1)iAp - (y - 1)pAi = 0
after some simplification. fip =
(y - l)E,
Thus the expressions in (3.10a)-(3.10d)
(3.25) give, in all cases,
fii = (y - l)P.
Therefore, substituting (3.8) into (3.14) and using these averages gives, after simplification use of $= $ + +ApAi,
(3.26) and
(3.27) Similarly, the expression Is/P + i + aApAi/p- in (3.13a) becomes >/(y - 1) where ;;?= yp/p, and t - j?j?,/rYi in (3.13b) is zero, leaving $L2 - $(Au)~/(~ - 1). Finally, in (3.17b), l__=-(Au)’ 4a’*
YPh
G2
=-
aL a’2 a
34
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
4. Test problems and numerical results We consider five test problems for various equations of state to assess the scheme in Section 3. Problem 1. This is the well-known shock tube problem of Sod [6] for the Euler equations with
y = 1.4 and initial data
PY U7 p = i
1, 0, 1,
r< +,
0.125, 0, 0.1,
r> $.
The main features of the exact solution are a shock moving to the right followed by a contact discontinuity, also moving to the right, but more slowly, and an expansion fan moving to the left.
u 0.99..
0.M
..
0.33 . .
0.3 0.9
0.6
0.P
X
.
4.33
3.46
0.w
0.3
0.6
0.9
X
0.3
0.6
0.9
X
.
4.63
..
-0.w..
I 2.61..
1.73
..
0.w
-0.36
an
-I
.oo
-0.119
+
I
-1.73
-2.61..
Fig. 1. Solution
of Problem
1 at t = 0.144 s.
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
I
*: “’
““..__
..‘\
x
s
.,...
..z.
“1
9 0
a
35
36
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
x
e x 9 D
;
a 3
D
”
yl
-
c x
9
.,.,_:... e.
.D
”
----.
L..
x
P
x 9 D
x 9 -
:, j..
.__,
I
P. Glaister/Applied
Numerical Mathematics 15 (1994) 27-52
37
Fig. 1 shows the approximate and exact solutions for p, u, p, and e at t = 0.144 s using 100 mesh points, and a good agreement between the two can be seen. Problem 2. This problem is usually described as “Two Interacting
Blast Waves” and is a shock tube problem for the Euler equations in slab symmetry with an ideal gas with y = 1.4 and initial data P, u, P =
1, 0, 1000, 1, 0, 0.1, 1, 0, 100,
0
The walls of the tube at r = 0 and Y= 1 are assumed to be perfectly reflecting. Two strong blast waves develop and collide producing a complex flow. A detailed description of the time evolution of the flow can be found in [8]. Figs. 2 and 3 show the approximate solution for p and u at t = 0.01, 0.016, 0.026, 0.028, 0.03, 0.032, 0.034 and 0.038 seconds using 400 mesh points, and the results compare well with those given by Woodward and Collella [S] using a more complicated Riemann solver.
JTH SIAB SYnnETRY - Shock Reflection PI
” 0.w
KEY 0.44
p 0.33
”
-
P-
Denslry
Velocity Pl-esEU-e
I - Internal energy
-
Exact solution Approxlmete solutlor
equat,on of state for an Ideal c~as 50 fleshpotnts 89 Time steps Ax = 0.02 At = 0.0065 Pressure rat,0 = 10 'Superbee' llmlter used
Reflected Boundary Condrtlons etr=O
Fig. 4. Solution
of Problem
3: ideal equation
of state.
38
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
Problem 3. This problem is concerned with shock reflection in one dimension of a gas governed by the Euler equations in slab symmetry. We consider a region 0 < Y< 1 divided into fifty equally spaced mesh points and the initial conditions are (p, U, i, p> = (po, uO, i,, p(pO, i,)). This represents a gas of constant density and pressure moving towards Y= 0. The boundary at Y= 0 is a rigid wall and the exact solution describes shock reflection from the wall. Five equations of state are chosen:
(4 ideal gas P = (y -
W,
where we take the ratio of specific heat capacities monatomic gas; (b) stiffened gas (sometimes used for a liquid) P =B(P/P,
- 1) +
(Y -
to be 5/3 corresponding
y
to a
l)pi,
where we take B = 1 and y = 5/3;
I& SYHHFIBT " 0.w
SDLUTION OF THE Fill FR EOUATIONS I
P
- Shock
Reflection
0.66 0.33
p -
Denslry
”
Velocity
-
PreSSWe
P-
I -
q3
0.6
0.P
Internal
energy
x Exact
-
-03..
soLut!on
~OOQIO Approximate
eolutlor
au... av9..
et~ffened
equatlon
of stare 50 flesh po,nts 64 Time
I ’ b
P ?.P..
I
steps
Ax = 0.02 At = 0.0046 Pressure
,.u ..
‘Superbee’
ratlo llmtter
=
IO used
0.1...
0.5
0.b
0.9
x
-0.7. -I.&n Reflected
Boundary atx=O
-2.a..
i_ at
time
t
=
Fig. 5. Solution
0.292
of Problem
3: stiffened
equation
of state.
Condltlons
39
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
(c) copper (due to R.K. Osborne, see [4]) p=
(& IM
a,+%9
+Jqb,+5(b,
+U)
+Jqco+c,6))],
where E =p,i and 5 = p/p,, - 1 and the constants a,, a2, b,, b,, b,, cO, cl, and 4. can be found in [2]; (4 covolume (used in connection with combustion environments) p = (7 - l)Pi 1-pb ’ where we take b = 0.8 and y = 1.4; (e) two molecular vibrating gas p + (Y - GLp = (y - qp’ 1, eP”/P - 1 where we take y = 1.4 and p = 0.5902 which corresponds
ON P
OF
THF
FIIIw
I
I
LAB " 0.w
- Shock
to molecular oxygen at 2000 K.
Reflect,on
KEI 0.66
p 0.33
”
-
P0.3
0.6
0.9
x
(13
0.6
0.9
VetocIty Pr08sUfW
I -
I
oens1ty
Internal
energy
x -
4.33..
Exact
solutron
oooooo Approxlmare
solurlor
-0.M.. 4.w.-
equat,on
of etate
for copper
-
I . 0.72 .-
P 20. 1I.J
O.&Y
6.7
0.24 ..
-
50 Neeh
points
54 Time
steps
Ax = 0.02 At = 0.0053 Pressure 'Superbee'
rat,0
= IO
L,n,ter
used
i 0.1
d.7
4.11
-I,.,
-0.40
-2o.I
4.72
0.6
0.9
x
Reflected
etx=O
-
at t,me t =
Fig. 6. Solution
u.
aI>
of Problem
Boundary
3: equation
of state for copper.
Condlt,ons
40
P. Glaister /Applied
Numerical Mathematics
T
I5 (1994) 27-52
- Shock Reflection
p -
Denslry
Velocity p - Preesure I - Interna 1 energy u -
Exact solut,on ooooo~Approxlmete solutlor
covolume equat Ion of state 50 Mesh pointa 95 Time steps Ax = 0.02 At = 0.0066 Pressure rat10 = IO 'Superbee' llmlter used
Reflected Boundary Condltlons atx=O et time t =
Fig. 7. Solution
0.628
of Problem
3: covolume
equation
of state.
For each equation of state we take u0 = - 1, and take p0 = 8.9 for (c), p0 = 1 for (b), and p0 = 0.1 for (a), (d), and (e). The initial condition for i, (in each case) is chosen to correspond to a shock strength p+/po = 10 where p+ denotes the pressure behind the shock and p. =p(p,,, io) denotes the pressure ahead of the shock. The results for these cases are given in Figs. 4-8, together with the exact solution when the shock has moved a distance of 0.3. Problem 4. This problem is concerned with a converging cylindrical shock and we consider a region 0 < r < 200 for the cylindrically symmetric case given by the Euler equations (2.2a) with S(r) = Y. Initially, a cylindrical diaphragm of radius r = 100 separates two uniform regions of gas at rest. The initial conditions are p = 4, p = 4 in the outer region, and p = 1, p = 1 in the inner region. When the diaphragm is removed at t = 0, a converging shock wave followed by a converging contact discontinuity moves towards the axis, Y = 0, and a diverging rarefaction wave moves outwards. The shock accelerates as it approaches the axis of symmetry and is reflected from the axis, interacts with the contact discontinuity (still converging), which results in a transmitted shock, a converging contact discontinuity and a weak converging reflected shock. Four equations of state are considered, namely (a), (b), Cd), and (e) associated with Problem 3 above.
P. Glaister /Applied
N OF
THF
FillFR
Numerical
1 Ai SY" 0.w
FWAI_IQM
‘ir.
Mathematics
- Shock
15 (1994) 27-52
41
Reflectron
KEY
p-
Density
u -
Veloclry
p - Pressure 1 - Internal -
Exact
energy
solution
ooo~x, Approximate
eq”atl0”
eolutlor
of *fete
for oxygen 50 Mesh
polnrs
142 TIM
steps
Ax = 0.02 At = 0.0068 Pressure 'Superbee'
ret10
= IO
llmlrer
used
0 Reflected
Boundary
Condltlons
atx=O
at
t,me
t =
Fig. 8. Solution
0.769
of Problem
3: equation
of state for oxygen.
Figs. 9-12 refer to the four equations of state using 200 mesh points. Three output times have been featured. For each case, we observe the development of the shock and the contact discontinuity, the reflection of the shock, and the interaction of the shock with the contact discontinuity. The results for equation of state (a) compare favourably with results produced using a more complex algorithm [l]. Problem 5. The final problem is concerned with a converging and a diverging cylindrical shock in 0 < Y G 200. Initially, two cylindrical diaphragms of radii r = 50 and Y = 150 separate three uniform regions of gas at rest. The initial conditions are p = 4, p = 4 in the inner and outer regions, and p = 1, p = 1 in the middle region. When the diaphragms are removed at t = 0, a converging shock wave, followed by a converging contact discontinuity moves towards the axis, r = 0, and a diverging shock wave, followed by a diverging contact discontinuity moves away from the axis. The shocks subsequently interact, resulting in a diverging shock wave weakening in strength, together with a converging shock wave increasing in strength. Each of these shocks then interact with the corresponding contact discontinuity, as in Problem 4, and this results, for each interaction, in a transmitted shock, a weak reflected shock and a contact discontinuity. Four equations of state are again considered, namely (a), (b), (d), and (e).
42
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
~JL. .-..... ..... ...... 54
8.A) E
“‘\
P. Glaister/Applied
Numerical
Mathematics
15 (1994) 27-52
c:, \
8 \
0 \P-’
!
a
I
2
43
.
5
L
I
1
‘.
!
0
“*.‘i
JI
A
N
d
d
d.
.M L
0 8
8
.
E
.
44
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
:
9
i
i
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
45
L
iL..
’ . . .._ -....
s
..
*
8
.
\
”
9
1 R +
Fi +
46
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
J---+.... -a....
.. .
s E
\
s.c .
\ .-........ ..
.‘....._
s
0
.. E .
. \
:
*
.. % / .:.... . .. .._
. . . . . . -.....
3.3
If!
sl d
i
.. s .. XJ
PI
9
+
r
ci
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
J
E d
t d
i% d
L
a-.._.... ...
.-.
s 4
.
\ .
8 *\
L ..
w ..
-a.......* .’
s_
.. 51; . . .
. ‘.
/.
,
f
. . . . .
b
\
8
. . . . . . . . .. 0
2
‘y 9
i
i
48
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
L
..
-8”*...\.. s
..I %. S’ s i*......_-..., ..... -,,~
L
1*
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
49
L
I s! *: :.
..I .
\ l-L_ 8
/
s
L;
it
n 2
“! 7
.
‘K. *-*..._ -..-.....* i
1
97
d +
i
9
I
1 I
iiN..
-..
s E
.\’ % / . ..::.* -........ -....Q.....
.....
a d
2
fj
I
.
.
.
/ / ,:.... ... ..-..... *.....
s
-.....
-.
\
*
P. Glaister/Applied
50
Numerical Mathematics 15 (1994) 27-52
Figs. 13-16 refer to the four equations of state using 200 mesh points. Three output times have again been featured. For each case, we observe the development of the shocks and contact discontinuities, the interaction of shocks, and the subsequent interaction of each shock with a contact discontinuity, resulting in transmitted and reflected shocks and contact discontinuities.
5. Conclusions
We have presented an efficient shock-capturing algorithm for the Euler equations for real gases in a duct of variable cross-section. We have utilised the arithmetic mean for the averages of the flow variables in computational cells. Numerical experiments typically show a 10% decrease in the computer time required when the square root averages are replaced by the arithmetic mean in the case of an ideal gas. The main features of this scheme are (i) it can be used for slab, cylindrically or spherically symmetric problems with confidence and (ii) it can be used as a comparison with results from two-dimensional schemes by choosing a large number of mesh points for accuracy and not be expensive on computing. The scheme can also be used for general duct flows.
Appendix A
We derive here_the matrix k in Eq. (3.5). To determine A we first write Aw and AF in terms of AU, where w = (P, PU, e)T,
F= (PU, P +Pu2,
u(e +P))~,
24= (P, u, P)‘.
Following the identities AP = PP,
(A4
A
(f-w W.3) W)
= PAu + UAp,
A(pi) = PAi + iAp, A( pu2) = ZAp + pAu2 = ZAp + 2ptiAu, where P = 3PL + PR)? u = gu,
+ UR),
E= +(iL + iR),
the arithmetic mean of left and right states, 2 = ;(ut
+ u;>
(A-8)
the arithmetic mean of the square of the velocity, we can write Aw=tiAu.
(A-9)
P. Glaister /Applied Numerical Mathematics 15 (1994) 27-52
where 1 ii
i=
j+ ;z I
0
P
0.
pu
#z
Similarly, Ap =j?,Ap
1
0
51
(A.lO)
+fiiAi,
(A.ll)
where fiP and ci are appropriate averages of the derivatives satisfy this equation exactly (see Section 3.31,
of the equation
of state chosen
to
A(up) = UAp +pAu,
(A.12)
P=;(PL+PR)
(A.13)
where
and A(pui)
= piAu + UA(pi)
-= piAu + updi
+ EEAp,
(A.14)
where z = +(pLiL + pRiR).
(A.15)
Finally, A( pu3) = u3Ap + ,5Au3 = u3Ap + 3,$Au,
(A.16)
z = +(u;
(A.17)
where + u;),
and ‘?_ 1 u - &L;+ULUR+U;)
(A.18)
is an average of the square of the velocity. (N.B. the choice in (A.14) is made so that the eigenvalues of the approximate Jacobian k have the simplest possible form. In particular, k has fi as one eigenvalue.) Combining (A.21, (A.4), (A.111, (A.121, (A.14) and (A.16) gives AF=eAu,
(A.19)
where
c=
u
P
0
fiP + 2
2pu
fii
ii& + uE + $3
p+ ;pu2+ ! and thus from (A.91 and (A.191 AF = &-‘Arv.
i#+pu
1
(A.20)
(A.21)
P. Glaister/Applied Numerical Mathematics15 (1994) 27-52
52
Therefore a solution of (3.4) is obtained with the approximate the identities
Jacobian A’ = &‘,
and using (A.22)
2U”+z=3>,
(A.23)
2ii3 + z = 3&P,
(A.24)
where ii = ULUR
(A.29
is the geometric mean of left and right states, gives rise to the matrix in (3.5), where (A.26)
References [II
M. Ben-Artzi and J. Falcovitz, A second order Godunov-type scheme for compressible fluid dynamics, .I Comput. Phys. 55 (1984) l-32. Dl P. Glaister, An approximate linearised Riemann solver for real gases, J. Comput. Phys. 74 (1988) 382-408. open channel flows, Appl. Numer. [31 P. Glaister, An efficient numerical method for subcritical and supercritical Math. 11 (1993) 497-508. [41 T.D. Riney, Numerical Evaluation of Hypervelocity Impact Phenomena (Academic Press, New York, 1970). Riemann solvers, parameter vectors and difference schemes, .I. Comput. Phys. 43 (1981) [51 P.L. Roe, Approximate 357-372. laws, .I. 161 G.A. Sod, A. survey of several finite difference methods for systems of nonlinear hyperbolic conservation Comput. Phys. 27 (1978) 1-31. flux-vector splitting and Roe average for an equilibrium real gas, J. Comput. Phys. 89 [71 M. Vinokur, Generalised (1990) 276-300. fluid flow with strong shocks, J. 181 P.R. Woodward and P. Colella, The numerical simulation of two-dimensional Comput. Phys. 54 (1984) 115-173.