Mathematical Models
and
Numerical Methods
for
Compressible Multicomponent Flow
Keh-Ming Shyue
Department of Mathematics
National Taiwan University
What is Multiphase Flow ?
General features
Existence of more than one phase or component
(fluid-like or not) in an underlying flow environment
separated by interfaces
Problems of interest
Typical Multiphase Flow Patterns
Examples taken from
Fundamentals of Multiphase
Solid Plate Impact
a) phi = 25 b) phi = 30
c) phi = 45 d) phi = 60
The collision of two aluminum plates inclined at an angle of phi degrees to one another. Only the plate on the right is shown, the left hand boundary is the line of impact. Color contours of density are shown: pale blue is vacuum, green is unshocked aluminum, and dark blue is jetted material. Note that the onset of jetting occurs between phi = 25 and phi = 30 degrees. These computations are of the experiments reported by J.M. Walsh, R.G. Shreffler and F.J. Willig in "Limiting Conditions for Jet Formation in High Velocity Collisions", J. Appl. Phys. 24:349--359 (1953).
Shock Wave in Bubbly Flow
Density
a)
Pressure
b)
c)
d)
e)
Two-Phase Flow Model
Saurel & Gallouet (1998)
(α
1
ρ
1
)
t
+ ∇ · (α
1
ρ
1
~u
1
) = ˙
m
(α
1
ρ
1
~u
1
)
t
+ ∇ · (α
1
ρ
1
~u
1
⊗ ~u
1
) + ∇(α
1
p
1
) =
p
0
∇α
1
+ ˙
m
~u
0
+ F
d
(α
1
ρ
1
E
1
)
t
+ ∇ · (α
1
ρ
1
E
1
~u
1
+ α
1
p
1
~u
1
) =
p
0
(α
2
)
t
+ ˙
m
E
0
+ F
d
~u
0
+
Q
0
(α
2
ρ
2
)
t
+ ∇ · (α
2
ρ
2
~u
2
) = − ˙
m
(α
2
ρ
2
~u
2
)
t
+ ∇ · (α
2
ρ
2
~u
2
⊗ ~u
2
) + ∇(α
2
p
2
) =
p
0
∇α
2
− ˙
m
~u
0
− F
d
(α
2
ρ
2
E
2
)
t
+ ∇ · (α
2
ρ
2
E
2
~u
2
+ α
2
p
2
~u
2
) = −
p
0
(α
2
)
t
− ˙
m
E
0
− F
d
~u
0
−
Q
0
(α
2
)
t
+
~u
0
· ∇α
2
= µ (p
2
− p
1
)
α
k
= V
k
/V
,
p
k
= ρ
k
R
k
T
k
,
E
k
= e
k
+ ~u
2
k
/2
, for
k = 1, 2
ρ =
P
k
α
k
ρ
k
,
p =
P
k
α
k
p
k
Two-Phase Flow Model (Cont.)
Model is derived using
averaging theory
of Drew
(Theory of Multicomponent Fluids, D.A. Drew & S. L.
Passman, Springer, 1999)
Namely, introduce
indicator function
χ
k
as
χ
k
(M, t) =
(
1
if
M
belongs to the phase
k
0
otherwise
Denote
< ψ >
as
volume averaged
for flow variable
ψ
,
hψi =
1
V
Z
V
ψ dV
Gauss
&
Leibnitz
rules
Two-Phase Flow Model (Cont.)
Take product of each balance law with
χ
k
, and perform
averaging process. In case of
mass conservation
equation,
for example, we have
hχ
k
ρ
k
i
t
+ ∇· < χ
k
ρ
k
~u
k
>= hρ
k
(χ
k
)
t
+ ρ
k
~u
k
· ∇χ
k
i
Since
χ
k
is governed by
(χ
k
)
t
+
~u
0
· ∇χ
k
= 0,
yielding the averaged equation for mass
hχ
k
ρ
k
i
t
+ ∇· < χ
k
ρ
k
~u
k
>= hρ
k
(~u
k
−
~u
0
)
· ∇χ
k
i
Analogously, we may derive averaged equation for
momentum, energy, & entropy (not shown here)
Two-Phase Flow Model (Cont.)
In summary,
hχ
k
i
t
+ h~u
k
· ∇χ
k
i = h(~u
k
−
~u
0
)
· ∇χ
k
i
hχ
k
ρ
k
i
t
+ ∇· < χ
k
ρ
k
~u
k
>= hρ
k
(~u
k
−
~u
0
)
· ∇χ
k
i
hχ
k
ρ
k
~u
k
i
t
+ ∇· < χ
k
ρ
k
~u
k
⊗ ~u
k
> +∇ hχ
k
p
k
i = hp
k
∇χ
k
i +
hρ
k
~u
k
(~u
k
−
~u
0
)
· ∇χ
k
i
hχ
k
ρ
k
E
k
i
t
+ ∇· < χ
k
ρ
k
E
k
~u
k
+ χ
k
p
k
~u
k
>= hp
k
~u
k
· ∇χ
k
i +
hρ
k
E (~u
k
−
~u
0
)
· ∇χ
k
i
Note existence of various
interfacial source terms
, &
modelling of these terms is essential and is a difficult issues
of multiphase flows
Baer-Nunziato Two-Phase Flow Model
Baer & Nunziato (J. Multiphase Flow 1986)
(α
1
ρ
1
)
t
+ ∇ · (α
1
ρ
1
~u
1
) = 0
(α
1
ρ
1
~u
1
)
t
+ ∇ · (α
1
ρ
1
~u
1
⊗ ~u
1
) + ∇(α
1
p
1
) =
p
0
∇α
1
+ λ(~u
2
− ~u
1
)
(α
1
ρ
1
E
1
)
t
+ ∇ · (α
1
ρ
1
E
1
~u
1
+ α
1
p
1
~u
1
) =
p
0
(α
2
)
t
+ λ
~u
0
· (~u
2
− ~u
1
)
(α
2
ρ
2
)
t
+ ∇ · (α
2
ρ
2
~u
2
) = 0
(α
2
ρ
2
~u
2
)
t
+ ∇ · (α
2
ρ
2
~u
2
⊗ ~u
2
) + ∇(α
2
p
2
) =
p
0
∇α
2
− λ(~u
2
− ~u
1
)
(α
2
ρ
2
E
2
)
t
+ ∇ · (α
2
ρ
2
E
2
~u
2
+ α
2
p
2
~u
2
) = −
p
0
(α
2
)
t
− λ
~u
0
· (~u
2
− ~u
1
)
(α
2
)
t
+
~u
0
· ∇α
2
= µ (p
2
− p
1
)
α
k
= V
k
/V
: volume fraction (
α
1
+ α
2
= 1
),
ρ
k
: density,
~u
k
: velocity,
p
k
: pressure,
E
k
= e
k
+ ~u
2
k
/2
: specific total
Baer-Nunziato Model (Cont.)
p
0
&
~u
0
: interfacial pressure & velocity
Baer & Nunziato (1986)
p
0
= p
2
,
~u
0
= ~u
1
Saurel & Abgrall (1999)
p
0
=
P
2
k=1
α
k
p
k
,
~u
0
=
P
2
k=1
α
k
ρ
k
~u
k
P
2
k=1
α
k
ρ
k
λ
&
µ
(> 0):
relaxation parameters
that determine rates at
which velocities and pressures of two phases reach
equilibrium
Mathematical Structure
Balance Laws with non-conservative products
Stiff relaxation term
Complicated Riemann problem solutions when the
model system is hyperbolic
Reduced model when
λ → 0
,
µ → 0
or
λ → ∞
,
µ → ∞
is
Reduced Two-Phase Flow Model
Murrone & Guillard (JCP 2005)
Assume
λ = λ
′
/ε
&
µ = µ
′
/ε
,
λ
′
= O(1)
&
µ
′
= O(1)
Use
formal asymptotic analysis
to show that when
ε → 0
leading order approximation of two-phase flow
model takes
(α
1
ρ
1
)
t
+ ∇ · (α
1
ρ
1
~u) = 0
(α
2
ρ
2
)
t
+ ∇ · (α
2
ρ
2
~u) = 0
(ρ~u)
t
+ ∇ · (ρ~u ⊗ ~u) + ∇p = 0
(ρE)
t
+ ∇ · (ρE~u + p~u) = 0
(α
2
)
t
+ ~u · ∇α
2
= α
1
α
2
ρ
1
c
2
1
− ρ
2
c
2
2
P
2
k=1
α
k
ρ
k
c
2
k
!
∇ · ~u
Reduced Two-Phase Flow Model (Cont.)
Some remarks about the reduced model:
1
. It can be shown
entropy
of each phase
S
k
now satisfies
DS
k
Dt
=
∂S
k
∂t
+ ~u · ∇S
k
= 0,
for
k = 1, 2
2
. Since product
α
1
α
2
is expected to be small, Shyue (JCP
1998), Allaire, Clerc, & Kokh (JCP 2002) proposed to
use
(α
2
)
t
+ ~u · ∇α
2
= 0
and now phase entropies satisfy
∂p
1
∂S
1
ρ
1
DS
1
Dt
−
∂p
2
∂S
2
ρ
2
DS
2
Dt
= ρ
1
c
2
1
− ρ
2
c
2
2
∇ · ~u
Reduced Two-Phase Flow Model (Cont.)
3
. Equation of State:
p = p(α
2
, α
1
ρ
1
, α
2
ρ
2
, ρe)
4
. Isobaric Closure:
p
1
= p
2
= p
For a class of EOS, explicit formula for
p
is available
(examples are given next)
For complex EOS, from
(α
2
, ρ
1
, ρ
2
, ρe)
in model
equations we recover
p
by solving
p
1
(ρ
1
, ρ
1
e
1
) = p
2
(ρ
2
, ρ
2
e
2
)
2
X
k=1
Reduced Two-Phase Flow Model (Cont.)
Polytropic Ideal Gas:
p
k
= (γ
k
− 1)ρ
k
e
k
ρe =
2
X
k=1
α
k
ρ
k
e
k
=
2
X
k=1
α
k
p
γ
k
− 1
=⇒
p
= ρe
2
X
k=1
α
k
γ
k
− 1
van der Waals Gas:
p
k
= (
γ
k
−
1
1
−b
k
ρ
k
)(ρ
k
e
k
+ a
k
ρ
2
k
) − a
k
ρ
2
k
ρe =
2
X
k=1
α
k
ρ
k
e
k
=
2
X
k=1
α
k
1 − b
k
ρ
k
γ
k
− 1
(
p
+ a
k
ρ
2
k
) − a
k
ρ
2
k
=⇒
p
=
"
ρe −
2
X
k=1
α
k
1 − b
k
ρ
k
γ
k
− 1
− 1
a
k
ρ
2
k
#
2
X
k=1
α
k
1 − b
k
ρ
k
γ
k
− 1
Reduced Two-Phase Flow Model (Cont.)
Two-Molecular Vibrating Gas:
p
k
= ρ
k
R
k
T (e
k
)
,
T
satisfies
e =
RT
γ − 1
+
RT
vib
exp T
vib
/T
− 1
As before, we now have
ρe =
2
X
k=1
α
k
ρ
k
e
k
=
2
X
k=1
α
k
ρ
k
R
k
T
k
γ
k
− 1
+
ρ
k
R
k
T
vib
,k
exp
T
vib
,k
/T
k
− 1
=
2
X
k=1
α
k
p
γ
k
− 1
+
p
vib
,k
exp
p
vib
,k
/
p
− 1
(
Nonlinear
)
Reduced Two-Phase Flow Model (Cont.)
5
. Model system is
hyperbolic
under suitable
thermodynamic stability condition (see below)
6
. It is easy to write the model system in unified
coordinates
7
. In the model
α
2
6= 0
or
1
, otherwise, either
ρ
2
or
ρ
1
is not
able to recover from
α
2
ρ
2
or
α
1
ρ
1
8
. This formulation of model equation would not work
when one fluid component is adiabatic, but the other
fluid component is not
9.
Surely, there are other set of model systems proposed
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −
V
2
(∂
2
p/∂V
2
)
S
(∂p/∂V )
S
,
S :
specific entropy
Assume fluid state satisfy
G > 0
for thermodynamic
stability,
i.e.
,
(∂
2
p/∂V
2
)
S
> 0
&
(∂p/∂V )
S
< 0
(∂
2
p/∂V
2
)
S
> 0
means
convex EOS
(∂p/∂V )
S
< 0
means
real speed of sound
, for
c
2
=
∂p
∂ρ
S
= −V
2
∂p
∂V
S
> 0
Thermodynamic Stability (Cont.)
When
G ≤ 0
, there exist
nonclassical nonlinear wave
&
phase transition
(active research area)
0.5 1 1.5 2 2.5 3 0.6 0.7 0.8 0.9 1 1.1 1.2 G = 0 G < 0 Saturation Curve Liquid−Vapor Coexistence Region
G > 0 C T = 1 T > 1 Vapor Region Liquid Region
