CA AA214 Ch5

AA214: NUMERICAL METHODS FOR COMPRESSIBLE FLOWS 1 / 63
AA214: NUMERICAL METHODS FOR

COMPRESSIBLE FLOWS
Representative Model Problems
1 / 63
Outline
1 Scalar Convection-Diffusion Equation
2 Burgers Equation
3 Inviscid Burgers Equation
4 Scalar Conservation Laws
Expansion Waves
Compression and Shock Waves
Contact Discontinuities
5 Riemann Problems
1D Riemann Problems for the Euler Equations
Riemann Problems for the Linearized Euler Equations
6 Roe’s Approximate Riemann Solver for the Euler Equations
Secant Approximations
Roe Averages
Algorithm and Performance
2 / 63
Scalar Convection-Diffusion Equation
Combines the convection (or advection) and diffusion equations to

describe physical phenomena where physical quantities are
transferred inside a physical system due to two processes, namely,
convection and diffusion
Convection is a transport mechanism of a substance or conserved
property by a fluid due to the fluid’s bulk motion
Diffusion is the net movement of a substance from a region of high
concentration to a region of low concentration
Also referred to by different communities as the drift-diffusion,
Smoluchowski, or scalar transport equation
∂c → − →
− →
−
+ ∇ · (~ac) = ∇ · (D ∇c) + S
∂t
where c is the variable of interest (species concentration for mass
transfer, temperature for heat transfer, · · · ), D is the diffusivity (or
diffusion coefficient), ~a is the average velocity of the quantity that is
moving, and S describes sources or sinks of the quantity c
3 / 63
Common simplifications
the diffusion coefficient D is constant, there are no sources or sinks
(S = 0), and the velocity field describes an incompressible flow
−
→ −
→
(∇ · ~
a = ∇ ·~ v = 0)
∂c
a · ∇c = D∇2 c
+~
∂t
in this form, the convection-diffusion equation combines both

parabolic and hyperbolic partial differential equations
stationary convection-diffusion equation
−
→ −
→
∇ · (~
ac) = ∇ · (D∇c) + S
4 / 63
Why is it a good representative model problem (and of what)?

for an incompressible flow (ρ = cst), the Navier-Stokes equations can
be written as

∂(ρ~
v) −
→ µ −
→
+~ v · ∇(ρ~v ) = ∇2 v ) + (f~ − ∇p)
(ρ~ (1)
∂t ρ
where
T
∇2 (ρ~
v) = ∇2 (ρvx ) ∇2 (ρvy ) ∇2 (ρvz )
−
→ − → −
→ − → −
→ − → T
= ∇ · ∇(ρvx ) ∇ · ∇(ρvy ) ∇ · ∇(ρvz )
= (∆(ρvx ) ∆(ρvy ) ∆(ρvz ))T
and f~ is a body force such as gravity
5 / 63
Why is it a good representative model problem? (continue)

for an incompressible flow (ρ = cst), the Navier-Stokes equations can
be written as

∂(ρ~
v) −
→ µ −
→
+~ v · ∇(ρ~v ) = ∇2 v ) + (f~ − ∇p)
(ρ~
∂t ρ
compare with the convection-diffusion equation when D is constant

−
→
and the velocity field describes an incompressible flow ( ∇ · ~
v = 0)
∂c −
→
a · ∇c = ∇2 (Dc) + S
+~
∂t
=⇒ the convection-diffusion equation mimics the incompressible

Navier-Stokes equations
6 / 63
Burgers Equation
Dropping the pressure term from the incompressible Navier-Stokes

equations (1) leads to

∂(ρ~v ) →
− µ
+ ~v · ∇(ρ~v ) = ∇ 2
(ρ~v ) + f~
∂t ρ
In one-dimension and assuming that µ is constant, the above

equation simplifies to Burgers equation (proposed in 1939 by the
dutch scientist Johannes Martinus Burgers)
∂vx ∂vx ∂ 2 vx
+ vx =ν + gx
∂t ∂x ∂x 2
µ fx
where ν = and gx =
ρ ρ
7 / 63
Burgers Equation
∂vx ∂vx ∂ 2 vx
+ vx =ν + gx
∂t ∂x ∂x 2
The above equation can be transformed into a linear parabolic

∂φ
equation using the Hopf-Cole transformation (vx = −2νφ ) then
∂x
solved exactly
This allows one to compare numerically obtained solutions of this
nonlinear equation with the exact one
For all these reasons, the Burgers equation is often used to
investigate the quality of a proposed CFD scheme for viscous (and
inviscid, see next) flows
8 / 63
Inviscid Burgers Equation
For ν = 0 and gx = 0, the Burgers equation simplifies to
∂vx ∂vx
+ vx =0
∂t ∂x
which is known as the inviscid Burgers equation

It is a prototype for equations whose solution can develop
discontinuities (shock waves)
It can be solved by the method of characteristics
It can be written in strong conservation law form as follows
v2
∂vx ∂( 2x )
+ =0
∂t ∂x
9 / 63
Consider the following inviscid Burgers problem
v2
∂vx ∂( 2x )
+ = 0
∂t ∂x
vxL if x < 0
vx (x, 0) = (2)
vxR if x > 0
10 / 63
Consider the following inviscid Burgers problem (continue)

consider now scaling x and t by a constant α > 0
x̄ = αx, t̄ = αt, α>0
since
∂ ∂ ∂ ∂
= α , and =α
∂t ∂ t̄ ∂x ∂ x̄
the inviscid Burgers equation is not affected by this scaling
furthermore, since the initial condition depends only on the sign of x,
it is not affected by the above scaling
=⇒ the inviscid Burgers problem defined above is scale invariant
11 / 63
Scale invariance often implies the risk of multiple solutions

if vx (x, t) is the solution of problem (2), then u(x, t) = vx (αx, αt) is
also a solution of problem (2) for any α > 0
hence, desiring uniqueness of the solution of the above problem is
desiring u ≡ vx — that is
x
vx (x, t) = v̄x
t
this implies that the solution vx (x, t) is constant on the rays
(characteristics) x = ct, and therefore the solution is said to be
self-similar 1
in a homework, it will be shown that more precisely, the solution of
problem (2) is x x
vx (x, t) = v̄x =
t t
this solution is called a rarefaction wave centered at the origin
(x = t = 0)
1 Self-similarity is the property of having a substructure analogous or identical to an
overall structure. For example, a part of a line segment is itself a line segment, and
thus a line segment exhibits self-similarity.
12 / 63
Self-similarity in nature
13 / 63
A rarefaction wave can be attached to a constant solution (for a

proof, look at the form of Eq. (2))
It can also join two constants
14 / 63
In many circumstances, the uniqueness of the solution is enforced by

imposing the condition that characteristics must impinge on a
discontinuity from both sides, which is known as the Lax Entropy
Condition
consider a shock located along the curve x = γ(t) and traveling at
dx dγ
the speed V = =
dt dt
let vx− (t) and vx+ (t) denote the left and right limits of the solution
vx (x, t) of problem (2), respectively
the Lax Entropy Condition states that
vx+ (t) < V < vx− (t)
(recall that the flow before a normal shock wave must be supersonic)
in particular, the Lax Entropy Condition states that the solution
must jump down
for problem (2), it can be shown that for α > 0, the solution jumps
up at the discontinuity: Thus, the only admissible solution — that is,
the solution in which any shock satisfies the Lax Entropy Condition
— is the continuous solution which has no shock
15 / 63
Scalar Conservation Laws
Scalar conservation laws are simple scalar models of the Euler

equations
They can be written in strong conservation form as
∂u ∂f (u)
+ =0 (3)
∂t ∂x
Their integral form in the space-time domain [x1 , x2 ] × [t 1 , t 2 ] is

Z x2 Z t2
[u(x, t 2 ) − u(x, t 1 )] dx + [f (u(x2 , t)) − f (u(x1 , t))] dt = 0
x1 t1
(4)
16 / 63
The solutions of the integral form (4) may contain jump

discontinuities: In this case, the discontinuous solutions are called
weak solutions of the differential form (3)
Jump discontinuities in the differential form (3) must satisfy a jump
condition derived from the integral form: For a jump discontinuity
traveling at a speed V , the jump condition is
f (u+ ) − f (u− ) = V (u+ − u− ) ⇔ Jf (u)K− −

+ = V JuK+ (5)
and therefore is analogous to the Rankine-Hugoniot relations recall

→
− →
− →
− T
J F ? K21 · ∇ ? g = 0, here with F ? (u) = (u f (u)) and
0
g (x, t) = x − xo − V (t − t )
17 / 63
Using chain rule, the non conservation form (or wave speed form) of
a scalar conservation law is
∂u ∂u
+ a(u) = 0
∂t ∂x
where
df
a(u) =
du
a(u) is called the wave speed
18 / 63
Examples
u2
f (u) = ⇒ Burgers equation
2
f (u) = au ⇒ linear advection
u2
f (u) = 2 , where c is a constant ⇒ Bucky-Leverett
u + c(1 − u)2
equation which is a simple model of two-phase flow in a porous
medium
19 / 63
Expansion Waves
Scalar conservation laws support features analogous to simple

expansion waves
For scalar conservation laws, an expansion wave (or a rarefaction
wave) is any region in which the wave speed a(u) increases from left
to right
a (u(x, t)) ≤ a (u(y , t)) , b1 (t) ≤ x ≤ y ≤ b2 (t)
A centered expansion fan is an expansion wave where all

characteristics originate at a single point in the x − t plane
hence,
the solution of problem (2) is a centered expansion fan
Centered expansion fans must originate in the initial
conditions or at intersections between shocks or contacts (see
definitions next)
20 / 63
Expansion Waves
21 / 63
Scalar conservation laws support features analogous to simple

compression and shock waves
For scalar conservation laws, a compression wave is any region in
which the wave speed a(u) decreases from left to right
a (u(x, t)) ≥ a (u(y , t)) , b1 (t) ≤ x ≤ y ≤ b2 (t)
A centered compression fan is a compression wave where all

characteristics converge on a single point in the x − t plane
The converging characteristics in a compression wave must
eventually intersect, creating a shock wave
22 / 63
Mean value theorem: Let f : [a, b] → R be a continuous function on the closed interval
[a, b], and differentiable on the open interval ]a, b[, where a < b. Then, there exists some
0 f (b) − f (a) 0
c ∈ ]a, b[ such that f (c) = ⇐⇒ f (b) − f (a) = f (c)(b − a) 2
b−a
f (b)−f (a) 1
Rb
2 Note that f 0 (c) = b−a
= b−a a f 0 (x)dx = average tangent
23 / 63
Mean value theorem: Let f : [a, b] → R be a continuous function on the closed interval
[a, b], and differentiable on the open interval ]a, b[, where a < b. Then, there exists some
0 f (b) − f (a) 0
c ∈ ]a, b[ such that f (c) = ⇐⇒ f (b) − f (a) = f (c)(b − a) 2
b−a
A shock wave is a jump discontinuity governed by the jump condition

f (u+ ) − f (u− ) = V (u+ − u− ) (see Eq. (5)): From the mean value theorem, it follows that
df
V = (ξ) = a(ξ), u− ≤ ξ ≤ u+
du
f (b)−f (a) 1
Rb
2 Note that f 0 (c) = b−a
= b−a a f 0 (x)dx = average tangent
23 / 63
A shock wave may originate in a jump discontinuity in the initial

conditions or it may form spontaneously from a smooth compression
wave
In addition to the jump condition (5), shock waves must satisfy
(think of the Lax Entropy Condition)
a(u− ) ≥ V ≥ a(u+ )
If wave speeds are interpreted as slopes in the x − t plane, then the

above equation implies that waves (characteristics) terminate on
shocks and never originate in shocks (shocks only absorb waves
— they never emit waves)
24 / 63
Contact Discontinuities
Scalar conservation laws support features analogous to contact

discontinuities – that is, surfaces that separate flow zones of
different density and temperature, but same pressure and velocity
For scalar conservation laws, a contact discontinuity is a jump
discontinuity from u− to u+ such that
a(u− ) = a(u+ )
Like contacts in the Euler equations, contacts in scalar

conservation laws must originate in the initial conditions or at
the intersections of shocks
25 / 63
Riemann Problems
In the theory of hyperbolic equations, a Riemann problem (named

after Bernhard Riemann) consists of a conservation law equipped
with uniform initial conditions on an infinite spatial domain, except
for a single jump discontinuity
In one-dimension (1D), for a hyperbolic problem governing the field
u, the Riemann problem centered on x = x0 and t = t 0 has the
following initial conditions

uL if x < x0
u(x, t 0 ) =
uR if x > x0
For example, problem (2) is a Riemann problem

For convenience, the remainder of this chapter uses x0 = 0 and
t0 = 0
26 / 63
Riemann Problems
In 1D, the Riemann problem has an exact analytical solution for the:
1) Euler equations; 2) scalar conservation laws; and 3) any linear
system of equations
Furthermore, the solution is self-similar (or self-preserving): It
stretches uniformly in space as time increases but otherwise retains
its shape, so that u(x, t 1 ) and u(x, t 2 ) are “similar” to each other
for any two times t 1 and t 2 — in other words, the solution depends
x
on the single variable rather than on x and t separately
t
The Riemann problem is very useful for the understanding of the
Euler equations because shocks and rarefaction waves may
appear as characteristics in the solution
Riemann problems appear in a natural way in finite volume methods
for the solution of equations of conservation laws due to the
discreteness of the grid: They give rise to the Riemann solvers which
are very popular in CFD
27 / 63
Riemann Problems
Shock Tube
Consider a 1D tube containing two regions of stagnant fluid at

different pressures
Suppose that the two regions are initially separated by a rigid
diaphragm
Suppose that this diaphragm is instantly removed (for example, by a
small explosion)
pressure imbalance ⇒ 1D unsteady flow containing a steadily moving
shock, a steadily moving simple centered expansion fan, and a
steadily moving contact discontinuity separating the shock and
expansion
the shock, expansion, and contact separate regions of uniform flow
28 / 63
Riemann Problems
29 / 63
Riemann Problems
Shock Tube
The flow in a shock tube has always zero initial velocity

Removing this restriction, the shock tube problem becomes a
Riemann problem and thus is a special case of the Riemann problem
Major result:
like the shock tube problem, the Riemann problem may give rise to a
steadily moving shock, a steadily moving simple centered expansion
fan, and a steadily moving contact separating the shock and
expansion; and the shock, expansion, and contact separate regions of
uniform flow
unlike the shock tube however, one or two of these waves may be
absent
30 / 63
Riemann Problems
Governing Equations
∂W ∂Fx
+ = 0
∂t ∂x
T
W = (ρ ρvx E )T , Fx = ρvx ρvx2 + p (E + p)vx
(
WL = (ρL ρL vxL EL )T if x < 0
W (x, 0) =
WR = (ρR ρR vxR ER )T if x > 0
(6)
v2

with p = (γ − 1) E − ρ x , and the speed of sound c given by
2
2 γp
c =
ρ
31 / 63
Riemann Problems
Exact Solution
First, consider the shock

in a frame moving with the (steadily moving) shock, the
Rankine-Hugoniot conditions3 can be written as
ρ2 (vx2 − V ) = ρ1 (vx1 − V ) (7)

ρ2 (vx2 − V )2 + p2 = ρ1 (vx1 − V )2 + p1 (8)
(E2 + p2 )(vx2 − V ) = (E1 + p1 )(vx1 − V ) (9)
where V is the speed of the shock

recall the expression of the speed of sound
p
c2 = γ
ρ
3 In −
→
this case, the Rankine-Hugoniot conditions are given by J F x K21 · e~x = JFx K21 = 0
32 / 63
Riemann Problems
Exact Solution
First, consider the shock (continue)

from the Rankine-Hugoniot conditions applied in a frame moving
with the (steadily moving) shock it follows that

2
γ+1 + p2
c2 p2 γ−1 p1
= (10)
c12

p1 1 + γ+1 p2
γ−1 p1

p2
c1 p1
−1
vx2 = vx1 + r (11)
γ γ+1

p2

2γ p1
− 1 + 1
s
γ+1 p2
V = vx1 + c1 −1 +1 (12)
2γ p1
p2
=⇒ functions of the ratio
p1
33 / 63
Riemann Problems
Exact Solution
Next, consider the contact discontinuity

by definition
vx3 = vx2 (13)
p3 = p2 (14)
Finally, consider the simple centered expansion fan
recall that for the 1D Euler equations, an expansion wave is a wave
where the wave speed (vx , or vx ± c) increases monotonically from
left to right
recall that a simple wave is a wave where all states lie on the same
integral curve of one of the characteristic families ⇒ a simple wave
in the entropy characteristic family ξ0 is a wave where
dξ+ = dξ− = 0 and therefore vx = cst and p = cst ⇒ entropy
waves cannot create expansions
it follows that the simple centered expansion fan here is a simple
centered acoustic fan associated with the characteristic curve
dx = (vx − c)dt
34 / 63
Riemann Problems
Finally, consider the simple centered expansion fan (continue)

along the integral curve of a simple centered expansion fan
associated with the characteristic curve dx = (vx − c)dt, the two
Riemann invariants ξ0 and ξ+ are constant
√
Z
dp γ−1 dp 2c
dξ0 = dρ − 2 = 0 ⇒ p = cstργ and c = cstγρ 2 ⇒ =
c ρc γ−1
dp 2c
dξ+ = dvx + = 0 ⇒ ξ+ = vx + = cst for dx = (vx + c)dt
ρc γ−1
2c
(and ξ− = vx − for dx = (vx − c)dt)
γ−1
hence, along the integral curve of a simple centered expansion fan

associated with the characteristic curve dx = (vx − c)dt and on this
characteristic curve
2c 2c
s = cst, vx + = cst, and vx − = cst
γ−1 γ−1
therefore in this flow region, all flow properties are constant and
dx = (vx − c)dt becomes the straight line x = (vx − c)t + cst
35 / 63
Riemann Problems
Finally, consider the simple centered expansion fan (continue)

now, along the integral curve of a simple centered expansion fan
associated with the characteristic curve x = (vx − c)t + cst
2c 2c4
vx + = vx4 +
γ−1 γ−1
hence along this integral curve and on the characteristic curve

x
x = (vx − c)t ⇔ c = vx − the following holds
t

2 x 2c4
vx + vx − = vx4 +
γ−1 t γ−1
γ−1

 2 x

 vx (x, t) = + vx4 + c4


 γ +1 t 2
γ−1


 2 x x
=⇒ c(x, t) = + vx4 + c4 −
 γ+1 t 2 t

 2γ

 c γ−1
 p
 = p4 (from the isentropic relations)
c4
(15)
36 / 63
Riemann Problems
Combine
now the shock, contact, and expansion results to determine
p2 p4 pL
across the shock in terms of the known ratio =
p1 p1 pR
2c
simple wave condition vx + = cst implies
γ−1
2c3 2c4
vx3 + = vx4 + (16)
γ−1 γ−1
from the third of equations (15) and (16) it follows that

" γ−1 #
2c4 p3 2γ
vx3 = vx4 + 1− (17)
γ−1 p4
from (13), (14) and (16) it follows that

   
γ−1 γ−1
2c4 p2 2γ 2c4 p1 p2 2γ
vx2 = vx4 + 1 −  = vx + 1 − 
4
γ−1 p4 γ−1 p4 p1
(18)
37 / 63
Riemann Problems
p4
Solving equation (18) for gives
p1
− 2γ
γ−1

p4 p2 γ−1
= 1+ (vx4 − vx2 ) (19)
p1 p1 2c4

p2
Finally, combining (11) and (19) delivers the nonlinear equation in
p1
  − 2γ
 γ−1

p2
 −1
 
γ−1 

p4 p2 
vx − vx − c1 r p1 
= 1+ 
p1 p1  2c4  4 1
γ γ+1

p2

−1 +1 

 
2γ p 1
(20)

p2
which can be solved by a preferred numerical method to obtain and therefore p2
p1
38 / 63
Riemann Problems
Once p2 is found, equation (11) gives vx2 , equation (10) gives c2 ,

and equation (12) gives the speed of the shock V , which completely
determines the state 2
Then, equations (13) and (14) give vx3 and p3 and equation (16)
gives c3 , which completely determines state 3
Finally, the first, second, and third of equations (15) deliver vx , c,
and p inside the expansion, respectively
In some cases (depending on the values of WL and WR ), the
Riemann problem may yield only one or two waves, instead of three:
To a large extent, the solution procedure described above handles
such cases automatically
39 / 63
Riemann Problems
The exact solution of the Riemann problem (6) is (relatively)

expensive because finding p2 requires solving the nonlinear
equation (20)
To this effect, approximate Riemann problems are often constructed
as surrogate Riemann problems for the Euler equations
Here, the family of approximate Riemann problems based on a
linearization of problem (6) is considered in the general case of m
dimensions
40 / 63
Riemann Problems
Consider the linear Riemann problem
∂W ∂W
+A = 0
∂t ∂x
WL if x < 0
W (x, 0) = (21)
WR if x > 0
where A is an m × m constant matrix whose construction is
discussed in the next section
Assume that A is diagonalizable
A = Q −1 ΛQ, Λ = diag (λ1 , · · · , λm )
where Q and Λ are constant matrices, and that ri and li ,

i = 1, ..., m are its right and left eigenvectors, respectively
Ari = λi ri , AT li = λi li (or liT A = λi liT )
41 / 63
Riemann Problems
In the linear case, the change to characteristic variables dξ = QdW

simplifies to
ξ = QW
and leads to the following characteristic form of problem (21)
∂ξ ∂ξ
+Λ = 0
∂t ∂x
ξL = QWL if x < 0
ξ(x, 0) =
ξR = QWR if x > 0
The individual form of the above problem is

∂ξi ∂ξi
+ λi = 0, i = 1, ..., m
∂t ∂x
ξLi = liT WL

if x < 0
ξi (x, 0) = (22)
ξRi = liT WR if x > 0
42 / 63
Riemann Problems
43 / 63
Riemann Problems
44 / 63
Riemann Problems
Since λi is constant, the solution of problem (22) is trivial: For

m = 3, it can be written as (λ1 > λ2 > λ3 )
 x
 (ξL1 ξL2 ξL3 )T if < λ3
t





 x
 (ξL1 ξL2 ξR3 )T if λ3 < < λ2

x 
t
ξ(x, t) = ξ¯ = x (23)
t 
(ξ ξ ξ )T
if λ < < λ
L R R 2 1
 1 2 3
t




 x
 (ξR1 ξR2 ξR3 )T if

 > λ1
t
45 / 63
Riemann Problems

 x
t





 x
 (ξL1 ξL2 ξR3 )T if λ3 < < λ2

x 
t
ξ(x, t) = ξ¯ = x (23)
t 
(ξ ξ ξ )T
if λ < < λ
L R R 2 1
 1 2 3
t




 x

 > λ1
t
If ∆W = WR − WL , then ∆ξ = Q∆W , and ∆ξi = liT ∆W is often
referred to as the strength or amplitude of the i-th wave
45 / 63
Riemann Problems

 x
t





 x
 (ξL1 ξL2 ξR3 )T if λ3 < < λ2

x 
t
ξ(x, t) = ξ¯ = x (23)
t 
(ξ ξ ξ )T
if λ < < λ
L R R 2 1
 1 2 3
t




 x

 > λ1
t
If ∆W = WR − WL , then ∆ξ = Q∆W , and ∆ξi = liT ∆W is often
referred to as the strength or amplitude of the i-th wave
Let
T T T
∆ξ 1 = (∆ξ1 0 0) , ∆ξ 2 = (0 ∆ξ2 0) , ∆ξ 3 = (0 0 ∆ξ3 )
Note that the superscripts used above are NOT powers: They are
used only to distinguish each of the above vector quantities from the
scalar jump ∆ξi in the i-th characteristic
45 / 63
Riemann Problems
Noting that ∆ξ 1 + ∆ξ 2 + ∆ξ 3 = ∆ξ = ξR − ξL , the solution (23)

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1



 t
x 

ξ(x, t) = ξ̄ =
t 






46 / 63
Riemann Problems

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1


 t
 ξL + ∆ξ 3


x
ξ(x, t) = ξ̄ =
t 






46 / 63
Riemann Problems

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1


 t x
 ξL + ∆ξ 3 = ξR − ∆ξ 2 − ∆ξ 1

x  if λ3 < < λ 2 < λ 1
ξ(x, t) = ξ̄ = t
t 






46 / 63
Riemann Problems

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1


 t x
 ξL + ∆ξ 3 = ξR − ∆ξ 2 − ∆ξ 1

x  if λ3 < < λ 2 < λ 1
ξ(x, t) = ξ̄ = t x (24)
t 
 ξL + ∆ξ 3 + ∆ξ 2 = ξR − ∆ξ 1 if λ3 < λ 2 < < λ 1


 t x
 ξ + ∆ξ 3 + ∆ξ 2 + ∆ξ 1 = ξ

if λ3 < λ 2 < λ 1 <
L R
t
46 / 63
Riemann Problems

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1


 t x
 ξL + ∆ξ 3 = ξR − ∆ξ 2 − ∆ξ 1

x  if λ3 < < λ 2 < λ 1
ξ(x, t) = ξ̄ = t x (24)
t 
 ξL + ∆ξ 3 + ∆ξ 2 = ξR − ∆ξ 1 if λ3 < λ 2 < < λ 1


 t x
 ξ + ∆ξ 3 + ∆ξ 2 + ∆ξ 1 = ξ

if λ3 < λ 2 < λ 1 <
L R
t
And noting that Q −1 ∆ξ i = ∆ξi ri , i = 1, 2, 3, the solution (24) can

be rewritten in terms of the original variables W = Q −1 ξ as follows
46 / 63
Riemann Problems

can be rewritten as
 x
 ξL = ξR − ∆ξ 3 − ∆ξ 2 − ∆ξ 1 if < λ3 < λ2 < λ1


 t x
 ξL + ∆ξ 3 = ξR − ∆ξ 2 − ∆ξ 1

x  if λ3 < < λ 2 < λ 1
ξ(x, t) = ξ̄ = t x (24)
t 
 ξL + ∆ξ 3 + ∆ξ 2 = ξR − ∆ξ 1 if λ3 < λ 2 < < λ 1


 t x
 ξ + ∆ξ 3 + ∆ξ 2 + ∆ξ 1 = ξ

if λ3 < λ 2 < λ 1 <
L R
t
And noting that Q −1 ∆ξ i = ∆ξi ri , i = 1, 2, 3, the solution (24) can

be rewritten in terms of the original variables W = Q −1 ξ as follows
 x
 WL = WR − ∆ξ3 r3 − ∆ξ2 r2 − ∆ξ1 r1 if t < λ3 < λ2 < λ1

 WL + ∆ξ3 r3 = WR − ∆ξ2 r2 − ∆ξ1 r1 if λ3 < x < λ2 < λ1

x 
W = t x (25)
t  WL + ∆ξ3 r3 + ∆ξ2 r2 = WR − ∆ξ 1 r1 if λ3 < λ2 < < λ1


 t x
 W + ∆ξ r + ∆ξ r + ∆ξ r = W if λ < λ < λ <
L 3 3 2 2 1 1 R 3 2 1
t
46 / 63
Riemann Problems
Many CFD methods do not use the solution of a Riemann problem

directly, whether expressed in terms of ξ or W , but use instead only
the flux at x = 0
Here (linear Riemann problem), the flux function at x = 0 is AW (0)
From (25), it follows that

 AWL = AWR − ∆ξ3 λ3 r3 − ∆ξ2 λ2 r2 − ∆ξ1 λ1 r1 if 0 < λ3 < λ2 < λ1
AWL + ∆ξ3 λ3 r3 = AWR − ∆ξ2 λ2 r2 − ∆ξ1 λ1 r1 if λ3 < 0 < λ2 < λ1
AW (0) =
 AWL + ∆ξ3 λ3 r3 + ∆ξ2 λ2 r2 = AWR − ∆ξ1 λ1 r1 if λ3 < λ2 < 0 < λ1
AWL + ∆ξ3 λ3 r3 + ∆ξ2 λ2 r2 + ∆ξ1 λ1 r1 = AWR if λ3 < λ 2 < λ 1 < 0
47 / 63
Riemann Problems

the flux at x = 0

AW (0) =
Let λ− + + −
i = min(0, λi ) and λi = max(0, λi ) ⇒ λi − λi = |λi |
47 / 63
Riemann Problems

the flux at x = 0

AW (0) =
Let λ− + + −
Then, the flux function at x = 0 can be rewritten as
3
X
AW (0) = AWL + λ−
i ∆ξi ri
i=1
47 / 63
Riemann Problems

the flux at x = 0

AW (0) =
Let λ− + + −
3
X 3
X
AW (0) = AWL + λ−
i ∆ξi ri = AWR − λ+
i ∆ξi ri
i=1 i=1
47 / 63
Riemann Problems

the flux at x = 0

AW (0) =
Let λ− + + −
3
X 3
X
AW (0) = AWL + λ−
i ∆ξi ri = AWR − λ+
i ∆ξi ri
i=1 i=1
3
1 1 X
= A(WR + WL ) − |λi |∆ξi ri (26)
2 2
i=1
47 / 63
Riemann Problems
Note that
− −
λ+ +
i − λi = |λi | ⇒ Λ − Λ = |Λ|
− −
λ+ +
i + λi = λi ⇒ Λ + Λ = Λ
(Definitions) A+ = Q −1 Λ+ Q, A− = Q −1 Λ− Q, |A| = Q −1 |Λ|Q
A+ + A− = A, A+ − A− = |A|
48 / 63
Riemann Problems
Note that
− −
λ+ +
i − λi = |λi | ⇒ Λ − Λ = |Λ|
− −
λ+ +
i + λi = λi ⇒ Λ + Λ = Λ
A+ + A− = A, A+ − A− = |A|
It follows that
3
X 3
X
|λi | ∆ξi ri = Q −1 |λi |∆ξ i = Q −1 |Λ| ∆ξ = |A|(WR − WL )
| {z } |{z}
i=1 i=1
Q −1 ∆ξ i Q∆W
48 / 63
Riemann Problems
Note that
− −
λ+ +
i − λi = |λi | ⇒ Λ − Λ = |Λ|
− −
λ+ +
i + λi = λi ⇒ Λ + Λ = Λ
A+ + A− = A, A+ − A− = |A|
It follows that
3
X 3
X
|λi | ∆ξi ri = Q −1 |λi |∆ξ i = Q −1 |Λ| ∆ξ = |A|(WR − WL )
| {z } |{z}
i=1 i=1
Q −1 ∆ξ i Q∆W
Hence, the solution (26) can be written as

AW (0) = AWL + A− (WR − WL ) = AWR − A+ (WR − WL )
1 1
=⇒ AW (0) = A(WR + WL ) − |A|(WR − WL ) (27)
2 2
48 / 63
Roe’s Approximate Riemann Solver for the Euler Equations
Consider first any nonlinear scalar function f (w ), where w is also a

scalar variable, and let
df (w )
a(w ) =
dw
Two linear approximations of this function are
the tangent line approximation(s)
the secant line approximation
49 / 63
Tangent line approximations

about wR : f (w ) ≈ f (wR ) + a(wR ) (w − wR )
about wL : f (w ) ≈ f (wL ) + a(wL ) (w − wL )
These two approximations are more accurate near wR and wL ,
respectively
50 / 63
Secant line approximation
f (w ) ≈ f (wR ) + aRL (w − wR ) ⇔ f (w ) ≈ f (wL ) + aRL (w − wL )
where
f (wR ) − f (wL )
aRL =
(wR − wL )
It is more accurate on average over the entire region between wL
and wR
51 / 63
The mean value theorem connects tangent line and secant line
approximations as follows
aRL = a(η) for η between wL and wR
which essentially states that secant line slopes are average tangent
line slopes
52 / 63
Consider next any nonlinear vector function f (W ), where W is also

a vector
The tangent plane approximation about WL is defined as
f (W ) ≈ f (WL ) + A(WL )(W − WL )
df
where A = is the Jacobian matrix
dW
53 / 63
Consider next any nonlinear vector function f (W ), where W is also

a vector
The tangent plane approximation about WL is defined as
f (W ) ≈ f (WL ) + A(WL )(W − WL )
df
where A = is the Jacobian matrix
dW
A secant plane is any plane containing the line connecting WL and
WR : There are an infinite number of such planes
53 / 63
Secant plane approximations are defined as follows

f (W ) ≈ f (WL ) + ARL (W − WL )
where ARL is any matrix such that
f (WR ) − f (WL ) = ARL (WR − WL ) (28)
54 / 63
Secant plane approximations are defined as follows

f (W ) ≈ f (WL ) + ARL (W − WL )
where ARL is any matrix such that
f (WR ) − f (WL ) = ARL (WR − WL ) (28)
Note that if each of W and f (W ) is a vector with m components,

ARL is a matrix with m2 elements: Hence, equation (28) consists of
m equations with m2 unknowns
54 / 63
Example 1
 f (W ) − f (W ) 
1 R 1 L
0 0
 WR1 − WL1 
f2 (WR ) − f2 (WL )
 
 
 0 

 WR2 − WL2 

 f3 (WR ) − f3 (WL ) 
0 0
WR3 − WL3
Example 2
 f1 (WR ) − f1 (WL ) f1 (WR ) − f1 (WL ) f1 (WR ) − f1 (WL ) 

 WR1 − WL1 WR2 − WL2 WR3 − WL3 

 
f2 (WR ) − f2 (WL ) f2 (WR ) − f2 (WL ) f2 (WR ) − f2 (WL )
 
 
1 
 WR1 − WL1 WR2 − WL2 WR3 − WL3 
3



 

 f3 (WR ) − f3 (WL ) f3 (WR ) − f3 (WL ) f3 (WR ) − f3 (WL ) 

 WR1 − WL1 WR2 − WL2 WR3 − WL3 
55 / 63
f (WR ) − f (WL ) = ARL (WR − WL )
By analogy with the scalar case, suppose that one requires that in
the vector case, secant planes be average tangent planes: In this
case,
ARL = A(WRL )
where WRL is an average between WR and WL , and there are only m
unknowns — the components of WRL — that can be determined by
solving equation (28)
56 / 63
Roe Averages
Consider now the one-dimensional Euler equations: For these

equations, the conservative state vector W , flux vector Fx , and
∂Fx
Jacobian matrix A = can be written as
∂W
T 1 1 2 T
W = (ρ ρvx E ) = (ρ ρvx ρh + (γ − 1)ρvx )
γ 2γ
T
γ−1

2
T γ+1 2
Fx = ρvx ρvx + p (E + p)vx = ρvx ρh + ρvx ρhvx
γ 2γ
0 1 0
 
γ−3 2
A = 2 vx (3 − γ)vx γ−1  (29)
 

−vx h + 12 (γ − 1)vx3 h − (γ − 1)vx2 γvx
H
where h = is the specific enthalpy and H = E + p is the total
ρ
enthalpy per unit volume
57 / 63
Roe Averages
Choose ARL = A(WRL ): In this case, equation (29) leads to the
Roe-average Jacobian matrix
0 1 0
 
 
γ−3 2
(3 − γ)vxRL γ−1 
 
ARL = 2 vxRL

 
 
−vxRL hRL + 12 (γ − 1)vx3 hRL − (γ − 1)vx2 γvxRL
RL RL
(30)
Solving equation (28) using the above Roe-average Jacobian matrix

leads after several algebraic manipulations to
√ √
ρR vxR + ρL vxL
vxRL = √ √
ρR + ρL
H
√R +
H
√L √ √
ρR ρL ρR hR + ρL hL
hRL = √ √ = √ √
ρR + ρL ρR + ρL
Define √
ρRL = ρR ρL
s
1 2
=⇒ cRL = (γ − 1) hRL − v
2 xRL
58 / 63
Roe’s approximate Riemann solver
59 / 63
Roe’s approximate Riemann solver for the Euler equations (f = Fx ) is based on two ideas:
(1) the linear (secant) approximation of the flux vector
Fx (W ) ≈ Fx (WL ) + ARL (W − WL ) = Fx (WR ) + ARL (W − WR ) (31)
where ARL is the Roe-average Jacobian given in (30); and (2) the exact solution of the
linear Riemann problem (21) with A = ARL (see also (25) for A = ARL )
To this end, substituting (31) into the Euler equations (6) gives
∂W ∂ ∂W ∂W
+ {Fx (WL ) + ARL (W − WL )} = + ARL =0
∂t ∂x ∂t ∂x
From (31) and (27), it follows that Roe’s approximate Riemann solver computes the fluxes
at x = 0 as
Fx (W (0)) ≈ Fx (WL ) + ARL (W (0) − WL ) = Fx (WR ) + ARL (W (0) − WR )

1 1
≈ (Fx (WR ) + Fx (WL )) + ARL W (0) − ARL (WR + WL )
2 2
1 1 1
= (Fx (WR ) + Fx (WL )) + ARL (WR + WL ) − |ARL |(WR − WL )
2 2 2
1
− ARL (WR + WL )
2
bx (W (0)) = 1 (Fx (WR ) + Fx (WL )) − 1 |ARL |(WR − WL )

=⇒ F
2 2
60 / 63
Like the true (exact) Riemann solver, Roe’s approximate Riemann

solver yields three equally-spaced waves (see previous Figure)
Unlike in the true Riemann solver however, all three waves in Roe’s
approximate Riemann solver have zero spread (hence, Roe’s
approximate Riemann solver cannot capture the finite spread of the
expansion fan)
Roe’s approximate Riemann solver for the Euler equations is roughly
2.5 times faster than the exact Riemann solver
What about its accuracy?
61 / 63
Suppose that the exact Riemann problem yields a single shock or a

single contact with speed V (recall that unlike in the shock tube
problem, one or two of the shock, contact, and expansion waves may
be absent in the solution of the exact Riemann problem)
The shock or contact must satisfy the Rankie-Hugoniot conditions
Fx (WR ) − Fx (WL ) = V (WR − WL ) (⇒ V = cst)
For Roe’s approximate Riemann solver, ARL must satisfy the secant
plane condition
Fx (WR ) − Fx (WL ) = ARL (WR − WL )
It follows that
ARL (WR − WL ) = V (WR − WL )
which implies that V is a characteristic value (eigenvalue) of ARL

and WR − WL is a right characteristic vector (eigenvector) of ARL
62 / 63
Let V = λj and WR − WL = rj : then, the strength of the i-th wave

is given by

1 if i = j
∆ξi = liT (WR − WL ) = liT rj = δij (QQ −1 = I ) =
0 if i 6= j
In other words, two of the three waves have zero strength and the
single non trivial wave makes the full transition between WL and WR
3
at the speed V (recall ∆W = Q −1 ∆ξ =
P
ri ∆ξi )
i=1
63 / 63

is given by

1 if i = j
0 if i 6= j
3
P
ri ∆ξi )
i=1
It follows that for a single shock or a single contact, Roe’s
approximate Riemann solver — as a matter of fact, any secant plane
approximation — yields the exact solution!
63 / 63

is given by

1 if i = j
0 if i 6= j
3
P
ri ∆ξi )
i=1
It follows that for a single shock or a single contact, Roe’s
approximate Riemann solver — as a matter of fact, any secant plane
approximation — yields the exact solution!
Except in the above case however, Roe’s approximate Riemann
solver deviates substantially from the true Riemann solver: more
specifically, unlike the true nonlinear flux function, Roe’s linear flux
function allows for expansion shocks — that is, jump discontinuities
that satisfy the Rankine-Hugoniot relations but expand rather than
compress the flow and therefore violate the second law of
thermodynamics 63 / 63

CA AA214 Ch5

Uploaded by

Copyright:

Available Formats

CA AA214 Ch5

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

CA AA214 Ch5

Uploaded by

Copyright:

Available Formats

AA214: NUMERICAL METHODS FOR COMPRESSIBLE FLOWS 1 / 63

AA214: NUMERICAL METHODS FOR

Combines the convection (or advection) and diffusion equations to

in this form, the convection-diffusion equation combines both

Why is it a good representative model problem (and of what)?

and f~ is a body force such as gravity

Why is it a good representative model problem? (continue)

compare with the convection-diffusion equation when D is constant

=⇒ the convection-diffusion equation mimics the incompressible

Dropping the pressure term from the incompressible Navier-Stokes

In one-dimension and assuming that µ is constant, the above

The above equation can be transformed into a linear parabolic

For ν = 0 and gx = 0, the Burgers equation simplifies to

which is known as the inviscid Burgers equation

Consider the following inviscid Burgers problem

Consider the following inviscid Burgers problem (continue)

x̄ = αx, t̄ = αt, α>0

Scale invariance often implies the risk of multiple solutions

A rarefaction wave can be attached to a constant solution (for a

In many circumstances, the uniqueness of the solution is enforced by

vx+ (t) < V < vx− (t)

Scalar conservation laws are simple scalar models of the Euler

Their integral form in the space-time domain [x1 , x2 ] × [t 1 , t 2 ] is

The solutions of the integral form (4) may contain jump

f (u+ ) − f (u− ) = V (u+ − u− ) ⇔ Jf (u)K− −

and therefore is analogous to the Rankine-Hugoniot relations recall

Scalar conservation laws support features analogous to simple

a (u(x, t)) ≤ a (u(y , t)) , b1 (t) ≤ x ≤ y ≤ b2 (t)

A centered expansion fan is an expansion wave where all

Scalar conservation laws support features analogous to simple

a (u(x, t)) ≥ a (u(y , t)) , b1 (t) ≤ x ≤ y ≤ b2 (t)

A centered compression fan is a compression wave where all

A shock wave is a jump discontinuity governed by the jump condition

A shock wave may originate in a jump discontinuity in the initial

If wave speeds are interpreted as slopes in the x − t plane, then the

Scalar conservation laws support features analogous to contact

Like contacts in the Euler equations, contacts in scalar

In the theory of hyperbolic equations, a Riemann problem (named

For example, problem (2) is a Riemann problem

Consider a 1D tube containing two regions of stagnant fluid at

The flow in a shock tube has always zero initial velocity

First, consider the shock

ρ2 (vx2 − V ) = ρ1 (vx1 − V ) (7)

where V is the speed of the shock

First, consider the shock (continue)

Next, consider the contact discontinuity

Finally, consider the simple centered expansion fan (continue)

hence, along the integral curve of a simple centered expansion fan

Finally, consider the simple centered expansion fan (continue)

hence along this integral curve and on the characteristic curve

from the third of equations (15) and (16) it follows that

from (13), (14) and (16) it follows that

Once p2 is found, equation (11) gives vx2 , equation (10) gives c2 ,

The exact solution of the Riemann problem (6) is (relatively)

Consider the linear Riemann problem

A = Q −1 ΛQ, Λ = diag (λ1 , · · · , λm )

where Q and Λ are constant matrices, and that ri and li ,

Ari = λi ri , AT li = λi li (or liT A = λi liT )

In the linear case, the change to characteristic variables dξ = QdW