Numerical Methods for
Compressible Multicomponent Flow with
Moving Interfaces and Boundaries
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Illustrative Example 1
Flying Aluminum-plate problem Vacuum-Al boundary
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 0µs
Illustrative Example 1
Flying Aluminum-plate problem Vacuum-Al boundary
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 0.5µs
Illustrative Example 1
Flying Aluminum-plate problem Vacuum-Al boundary
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 1µs
Illustrative Example 1
Flying Aluminum-plate problem Vacuum-Al boundary
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 2µs
Illustrative Example 1
Flying Aluminum-plate problem Vacuum-Al boundary
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 4µs
Illustrative Example 2
Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 0ms
Illustrative Example 2
Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 1ms
Illustrative Example 2
Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 2ms
Illustrative Example 2
Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 3ms
Overview
Mathematical model
Fluid-mixture type equations of motion for homogeneous two-phase flow
Mie-Grüneisen EOS for real materials Numerical techniques
Finite volume method based on wave propagation Surface tracking for moving boundaries
Volume tracking for moving interfaces Numerical results
Future work
Two Phase Flow Problem
Ignore physical effects such as viscosity, surface tension, mass diffusion, and so on
Each fluid component k, k = 1, 2, satisties Eulerian conservation laws
ρt + ∇ · (ρ~u) = 0
(ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 Mie-Grüneisen equation of state
p(ρ, e) = pref(ρ) + ρ Γ(ρ)
he − eref(ρ) i
Mie-Grüneisen Equations of State
Typical examples are:
(pref, eref) lies along an isentrope
Mie-Grüneisen Equations of State
Typical examples are:
(pref, eref) lies along an isentrope
1. Jones-Wilkins-Lee EOS for gaseous explosives Γ(V ) = γ − 1
eref(V ) = e0 +
A V0
R1
exp −R1V V0
+ B V0
R2
exp −R2V V0
pref(V ) = p0 + A exp
−R1V V0
+ B exp −R2V V0
Here V = 1/ρ
Mie-Grüneisen Equations of State
Typical examples are:
(pref, eref) lies along an isentrope
2. Cochran-Chan EOS for solid explosives Γ(V ) = Γ0 = γ − 1
eref(V ) = e0 +
−A V0
1 − E1
"
V V0
1−E1
− 1
#
+ B V0
1 − E2
"
V V0
1−E2
− 1
#
pref(V ) = p0 + A
V V0
−E1
− B V V0
−E2
Mie-Grüneisen Equations of State
(pref, eref) lies along a Hugoniot locus
Assume linear shock speed us & particle velocity up us = c0 + s up
We may derive the relations Γ(V ) = Γ0
V V0
α
pref(V ) = p0 +
c02
(V0 − V ) [V0 − s(V0 − V )]2 eref(V ) = e0 +
1 2
hpref(V ) + p0
i(V0 − V )
Material Quantities for Model EOS
JWL EOS ρ0(kg/m3) A(GPa) B(GPa) R1 R2 Γ
TNT1 1630 371.2 3.23 4.15 0.95 0.30
TNT2 1630 548.4 9.375 4.94 1.21 1.28
Water 1004 1582 −4.67 8.94 1.45 1.17
CC EOS ρ0(kg/m3) A(GPa) B(GPa) E1 E2 Γ
TNT 1840 12.87 13.42 4.1 3.1 0.93
Copper 8900 145.67 147.75 2.99 1.99 2
Shock EOS ρ0(kg/m3) c0(m/s) s Γ0 α
Aluminum 2785 5328 1.338 2.0 1
Copper 8924 3910 1.51 1.96 1
Two-Phase Flow Model
Model derivation based on 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 phase k
0 otherwise
Denote < ψ > as volume averaged for flow variable ψ, hψi = 1
V Z
V
ψ dV Gauss & Leibnitz rules
hχk∇ψi = h∇(χkψ)i−hψ∇χki & hχkψti = h(χkψ)ti−hψ(χk)ti
Two-Phase Flow Model (Cont.)
Take product of each conservation (or balance) law with χk
& perform averaging process. In case of mass conservation equation, for example, we have
hχkρkit + ∇· < χkρk~uk >= hρk(χk)t + ρk~uk · ∇χki Since χk is governed by
(χk)t + ~u0 · ∇χk = 0, ~u0: interface velocity yielding averaged equation for mass
hχkρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki Analogously, we may derive averaged equation for momentum, energy, & entropy (not shown here)
Two-Phase Flow Model (Cont.)
In summary, averaged model system, we have, are
hχkρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki
hχkρk~ukit + ∇· < χkρk~uk ⊗ ~uk > +∇ hχkpki = hpk∇χki + hρk~uk (~uk − ~u0) · ∇χki
hχkρkEkit + ∇· < χkρkEk~uk + χkpk~uk >= hpk~uk · ∇χki + hρkE (~uk − ~u0) · ∇χki
hχkit + h~uk · ∇χki = h(~uk − ~u0) · ∇χki Existence of interfacial source terms
Mathematical as well as numerical modelling these terms are essential & difficult for multiphase flow problems
Homogeneous Two-Phase Flow Model
Assume homogeneous flow (i.e., across interfaces:
pk = p & ~uk = ~u, k = 0, 1, 2)
Introduce volume fraction αk = Vk/V (α1 + α2 = 1)
By dropping all interfacial terms, we may obtain a simplified model as
(αkρk)t + ∇ · (αkρk~u) = 0
(αkρk~u)t + ∇ · (αkρk~u ⊗ ~u) + ∇ (αkp) = p∇αk (αkρkEk)t + ∇ · (αkρkEk~u + αkp~u) = p~u · ∇αk (αk)t + ~u · ∇αk = 0
for k = 1, 2
Homogeneous Two-Phase Flow Model
Assume homogeneous flow (i.e., across interfaces:
pk = p & ~uk = ~u, k = 0, 1, 2)
Introduce volume fraction αk = Vk/V (α1 + α2 = 1) Alternatively, a simplified model as
(αkρk)t + ∇ · (αkρk~u) = 0
(ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 (αk)t + ~u · ∇αk = 0
Here ρ = P2
k=1 αkρk, ρE = P2
k=1 αkρkEk
Homogeneous Flow Model (Cont.)
Mixture equation of state: p = p(α2, α1ρ1, α2ρ2, ρe) Isobaric closure: p1 = p2 = p
For a class of EOS, explicit formula for p is available (examples are given next)
For some complex EOS, from (α2, ρ1, ρ2, ρe) in model equations we recover p by solving
p1(ρ1, ρ1e1) = p2(ρ2, ρ2e2) &
2
X
k=1
αkρkek = ρe This homogeneous two-phase model was called a five-equation model by Allaire, Clerc, & Kokh (JCP
2002) or a volume-fraction model by Shyue (JCP 1998)
Homogeneous Flow Model (Cont.)
Polytropic ideal gas: pk = (γk − 1)ρkek ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk p
γk − 1 =⇒ p = ρe
2 X
k=1
αk γk − 1 Van der Waals gas: pk = (1−γkb−1kρk )(ρkek + akρ2k) − akρ2k
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk 1 − bkρk γk − 1
(p + akρ2k) − akρ2k
=⇒
p =
"
ρe −
2
X
k=1
αk 1 − bkρk
γk − 1 − 1
akρ2k
# 2 X
k=1
αk 1 − bkρk γk − 1
Homogeneous Flow Model (Cont.)
Two-molecular vibrating gas: pk = ρkRkT (ek), T satisfies e = RT
γ − 1 + RTvib
exp Tvib/T − 1 As before, we now have
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk
ρkRkTk γk − 1
+ ρkRkTvib,k exp
Tvib,k/Tk
− 1
=
2
X
k=1
αk
p
γk − 1
+ pvib,k exp
pvib,k/p
− 1
(Nonlin. eq.)
Homogeneous Flow Model (Cont.)
It can be shown entropies, Sk, k = 1, 2, satisfy
∂p1
∂S1
ρ1
DS1
Dt − ∂p2
∂S2
ρ2
DS2
Dt = ρ1c21 − ρ2c22 ∇ · ~u Murrone & Guillard (JCP 2005) propsed a reduced two-phase flow model in which
(α2)t + ~u · ∇α2 = α1α2
ρ1c21 − ρ2c22 P2
k=1 αkρkc2k
!
and now phase entropies satisfy DSk
Dt = ∂Sk
∂t + ~u · ∇Sk = 0, for k = 1, 2
Homogeneous Flow Model (Cont.)
Model system is hyperbolic under suitable
thermodynamic stability condition (see below)
In the model, when α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1).
It is not absolutely clear in the model how to compute nonlinear term ρι, ι > 1 from αk & αkρk
This formulation of model equation would not work when one fluid component is adiabatic, but the other fluid component is not
Surely, there are other set of model systems proposed in the literature that are robust for homogeneous flow
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −V 2
(∂2p/∂V 2)S
(∂p/∂V )S , S : specific entropy
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −V 2
(∂2p/∂V 2)S
(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,
(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −V 2
(∂2p/∂V 2)S
(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,
(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0 (∂2p/∂V 2)S > 0 means convex EOS
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −V 2
(∂2p/∂V 2)S
(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,
(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0 (∂2p/∂V 2)S > 0 means convex EOS
(∂p/∂V )S < 0 means real speed of sound, for c2 = ∂p
∂ρ
S
= −V 2 ∂p
∂V
S
> 0
Finite Volume Wave Propagation Method
Finite volume method, QnS gives approximate value of cell average of solution q over cell S at time tn
QnS ≈ 1 M(S)
Z
S
q(X, tn) dV
M(S): measure (area in 2D or volume in 3D) of cell S
C E
D F
G H
Wave Propagation Method (cont.)
First order version: Piecewise constant wave update Godunov-type method: Solve Riemann problem at each cell interface in normal direction & use resulting waves to update cell averages
Qn+1S := Qn+1S −M (Wp ∩ S)
M(S) Rp, Rp being jump from RP
↓
↓
Wp Wp
Wave Propagation Method (cont.)
First order version: Transverse-wave included
Use transverse portion of equation, solve Riemann problem in transverse direction, & use resulting
waves to update cell averages as usual
Stability of method is typically improved, while conservation of method is maintained
↓ ↓
Wpq Wpq
Wave Propagation Method (cont.)
High resolution version: Piecewise linear wave update wave before propagation after propagation
a) b)
c) d)
αpr
p/2
αpr
p/2
λp∆ t
λp∆ t
Volume Tracking Algorithm
1. Volume moving procedure (a) Volume fraction update
Take a time step on current grid to update cell averages of volume fractions at next time step (b) Interface reconstruction
Find new interface location based on volume fractions obtained in (a) using an interface reconstruction scheme. Some cells will be subdivided & values in each subcell must be initialized.
2. Physical solution update
Take same time interval as in (a), but use a method to update cell averages of multicomponent model on new grid created in (b)
Interface Reconstruction Scheme
Given volume fractions on current grid, piecewise linear interface reconstruction (PLIC) method does:
1. Compute interface normal
Gradient method of Parker and Youngs Least squares method of Puckett
2. Determine interface location by iterative bisection
0 0
0 0
0 0
0 0
0 0
0
0 0
0 0
0 0
0 0
0
1 0.29 0.68 0.09
0.51 0.51
0.09 0.68 0.29
↓
↓ interface
interface
Data set Parker & Youngs Puckett
Volume Moving Procedure
(a) Volume fractions given in previous slide are updated with uniform (u, v) = (1, 1) over ∆t = 0.06
(b) New interface location is reconstructed
0
0 0 0 0
0 0 0 0
0 0
0 0
0 0
0 0
0
1 0.01
0.38 0.11
0.25 0.72
0.06 0.85
0.74 5(−3) 1(−3)
↑
↓
old interface
new interface
(a) (b)
Surface Moving Procedure
Solve Riemann problem at tracked interfaces & use resulting waves to find new location of interface at the next time step
o
o
o o
o
o
↑
↑ old front
old front
new front
Interface Conditions
For tracked segments representing rigid (solid wall)
boundary (stationary or moving), appropriate boundary states are assigned for fictitious subcells in each time step
For tracked segments representing material interfaces, jump conditions across interfaces are satisfied only in an approximate manner, & would not be imposed
explicitly in each time step
Stability Issues
Choose time step ∆t based on uniform grid mesh size
∆x, ∆y as
∆t maxp,q (λp, µq)
min(∆x, ∆y) ≤ 1,
λp, µq: speed of p-wave, q-wave from Riemann
problem solution in normal-, transverse-directions Use large time step method of LeVeque (i.e., wave
interactions are assumed to behave in linear manner) to maintain stability of method even in the presence of
small Cartesian cut cells
Apply interpolation operator (such as, h-box approach of Berger et al. ) locally for cell averages in irregular cells
References
(JCP 1998) An efficient shock-capturing algorithm for compressible multicomponent problems
(JCP 1999, 2001) A fluid-mixture type algorithm for
compressible multicomponent flow with van der Waals (Mie-Grüneisen) equation of state
(JCP 2004) A fluid-mixture type algorithm for barotropic two-fluid flow Problems
(JCP 2006) A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions
(Shock Waves 2006) A volume-fraction based algorithm for hybrid barotropic & non-barotropic two-fluid flow
problems
Flying Projectile Problem
0 1 2 3
−1
−0.5 0 0.5 1
Stationary Density
0 1 2 3
−1
−0.5 0 0.5 1
Pressure
0 1 2 3
−1
−0.5 0 0.5 1
Entropy
0 1 2 3
−1
−0.5 0 0.5 1
Moving Density
0 1 2 3
−1
−0.5 0 0.5 1
Pressure
0 1 2 3
−1
−0.5 0 0.5 1
Entropy
Cylinder lift-off Problem
Moving speed of cylinder is governed by Newton’s law Pressure contours are shown with a 1000 × 200 grid
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15 0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15 0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15 0.2
t = 0
t = 0.1641s
t = 0.30085s
Cylinder lift-off Problem
A convergence study of center of cylinder & relative mass loss for at final stopping time t = 0.30085s
Mesh size Center of cylinder Relative mass loss 250 × 50 (0.618181, 0.134456) −0.257528
500 × 100 (0.620266, 0.136807) −0.131474 1000 × 200 (0.623075, 0.138929) −0.066984
Results are comparable with numerical appeared in literature
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock in Air & R22 Bubble Interaction
Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble
Numerical schlieren images for density
Shock-Bubble Interaction (cont.)
Approximate locations of interfaces
time=55µs
air R22
time=115µs time=135µs
time=187µs time=247µs time=200µs
time=342µs time=417µs time=1020µs
Shock-Bubble Interaction (cont.)
Quantitative assessment of prominent flow velocities:
Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf Haas & Sturtevant 415 240 540 73 90 78 78 Quirk & Karni 420 254 560 74 90 116 82 Our result (tracking) 411 243 538 64 87 82 60 Our result (capturing) 411 244 534 65 86 98 76
Vs (VR, VT) Incident (refracted, transmitted) shock speed t ∈ [0, 250]µs (t ∈ [0, 202]µs, t ∈ [202, 250]µs ) Vui (Vuf) Initial (final) upstream bubble wall speed t ∈ [0, 400]µs (t ∈ [400, 1000]µs)
Vdi (Vdf) Initial (final) downstream bubble wall speed t ∈ [200, 400]µs (t ∈ [400, 1000]µs)
Shock wave in molybdenum over MORB
Numerical schlieren images for density
Shock-MORB Interaction (cont.)
Numerical schlieren images for pressure
Shock-MORB Interaction (cont.)
Approximate locations of interfaces
molybdenum
MORB
molybdenum
MORB
time = 50µs time = 100µs
Future Work
Extension to low Mach number flow
Remove sound-speed stiffness by preconditioning techniques or pressure-based method
Extension to include more physics towards real applications
Such as capillary, diffusion, or elastic-plastic effect Extension to 3D volume tracking method
Surface reconstruction
Finite volume method with moving interfaces Stability in presence of small cut cells
Extension to unified coordinates of W.-H. Hui