Wave-propagation based methods for
compressible homogeneous two-phase flow
Keh-Ming Shyue
Department of Mathematics
National Taiwan University
Outline
Eulerian formulation
Mathematical models
Wave-propagation based volume tracking method Sample examples
Generalized Lagrangian formulation Mathematical models
Flux-based wave decomposition method Sample examples
Future work
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 2/72
Two-Phase Flow Model (I)
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 )
Two-Phase Flow Model (II)
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 + mE ˙ 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 − mE ˙ 0 − F d ~u 0 − Q 0 (α 2 ) t + ~u 0 · ∇α 2 = µ (p 2 − p 1 )
˙
m : mass transfer, F d : drag force Q 0 : convective heat exchange
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 4/72
Two-Phase Flow 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
Two-Phase Flow Model: Derivation
Standard way to derive these equations is 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ψ∇χ k i & hχ k ψ t i = h(χ k ψ) t i−hψ(χ k ) t i
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 6/72
Two-Phase Flow Model (cont.)
Take product of each conservation law with χ k & 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 (~u 0 : interface velocity ), this leads to mass averaged equation for phase k
hχ ρ i + ∇· < χ ρ ~u >= hρ (~u − ~u ) · ∇χ i
Two-Phase Flow Model (cont.)
In summary, averaged model system, we have, are
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 hχ k i t + h~u k · ∇χ k i = h(~u k − ~u 0 ) · ∇χ k i
Note: existence of various interfacial source terms
Mathematical as well as numerical modelling of these terms are important (but difficult) for general multiphase flow
problems
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 8/72
Reduced Two-Phase Flow Model
Murrone & Guillard (JCP 2005)
Assume λ = λ ′ /ε & µ = µ ′ /ε , λ ′ = O(1) & µ ′ = O(1) Apply formal asymptotic analysis to Baer &
Nunziato’s model, as ε → 0 , gives leading order approximation
(α 1 ρ 1 ) t + ∇ · (α 1 ρ 1 ~u) = 0 (α 2 ρ 2 ) t + ∇ · (α 2 ρ 2 ~u) = 0
(ρ~u) t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 ( mixture momentum )
(ρE) t + ∇ · (ρE~u + p~u) = 0 ( mixture total energy )
Reduced Two-Phase Model (cont.)
Remarks:
1 . In this case, p 1 → p 2 & ~u 1 → ~u 2 , as ε → 0 , which means the flow is homogeneous ( 1 -pressure & 1 -velocity) with p ι = p & ~u ι = ~u , ι = 0, 1, 2 , across interfaces
2 . Mixture equation of state: p = p(α 2 , α 1 ρ 1 , α 2 ρ 2 , ρe) 3 . Isobaric closure: p 1 = p 2 = p
For some EOS, explicit formula for p is available (examples are given next)
For some other EOS, p is found by solving coupled equations
p 1 (ρ 1 , ρ 1 e 1 ) = p 2 (ρ 2 , ρ 2 e 2 ) & α 1 ρ 1 e 1 + α 2 ρ 2 e 2 = ρe
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 10/72
Reduced Two-Phase 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
Reduced Two-Phase 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 = ( 1−b γ k −1
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
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 11/72
Reduced Two-Phase 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 eq.)
Reduced Model: Remarks
4 . 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 5 . Model system is hyperbolic under suitable
thermodynamic stability condition
6 . When α k = 0 , ρ k can not be recovered from α k & α k ρ k , and so take α k ∈ [ε, 1 − ε] , ε ≪ 1
7 . Other model systems exist in the literature that are more robust for homogeneous flow (examples)
8 . When individual pressure law differs in form (see
below), new mixture pressure law should be devised first & construct model equations based on that
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 13/72
Barotropic & Non-Barotropic Flow
Fluid component 1 : Tait EOS
p(ρ) = (p 0 + B) ρ ρ 0
γ
− B
Fluid component 2 : Noble-Abel EOS
p(ρ, e) = γ − 1 1 − bρ
ρe
Mixture pressure law (Shyue, Shock Waves 2006)
(p 0 + B) ρ γ
− B if α = 1
Barotropic Two-Phase Flow
Fluid component ι : Tait EOS
p(ρ) = (p 0 ι + B ι )
ρ ρ 0 ι
γ ι
− B ι , ι = 1, 2
Mixture pressure law (Shyue, JCP 2004)
p =
(p 0 ι + B ι )
ρ ρ 0 ι
γ ι
− B ι if α = α ι (0 or 1) (γ − 1) ρ
e + B ρ 0
− γB if α ∈ (0, 1)
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 15/72
Homogeneous Two-Phase Model
In summary, mathematical model for compressible homogeneous two-phase flow:
Equations of motion
(α 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
(α ) + ~u · ∇α = α α ρ 1 c 2 1 − ρ 2 c 2 2 !
∇ · ~u
Wave Propagation Method
Finite volume formulation of wave propagation method, Q n S gives approximate value of cell average of solution q over cell S at time t n
Q n S ≈ 1 M(S)
Z
S
q(X, t n ) dV
M(S) : measure (area in 2 D or volume in 3 D) of cell S
C E
D F
G H
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 17/72
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
Q n+1 S := Q n+1 S − M (W p ∩ S)
M(S) R p , R p being jump from RP
W ↓ p
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
↓ ↓
W pq W pq
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 19/72
Wave Propagation Method (cont.)
High resolution version: Piecewise linear wave update wave before propagation after propagation
a) b)
c) d)
αpr
p/2
λ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)
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 21/72
Interface Reconstruction Scheme
Given volume fractions on current grid, piecewise linear interface reconstruction (PLIC) method does:
1. Compute interface normal
Gradient method of Parker & Youngs Least squares method of Puckett
2. Determine interface location by iterative bisection
0 0
0 0
0 0
0 0
0
0 .09 0 .51 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)
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 23/72
Surface Moving Procedure
Solve Riemann problem at tracked interfaces & use
resulting wave speed of the tracked wave family over ∆t to find new location of interface at the next time step
o
o
o o
o
↑
↑ old front
old front
new front
Boundary Conditions
For tracked segments representing rigid (solid wall)
boundary (stationary or moving), reflection principle is used to assign states for fictitious subcells in each time step:
z C := z E (z = ρ, p, α)
~u C := ~u E − 2(~u E · ~n)~n + 2(~u 0 · ~n)
~u 0 : moving boundary velocity
C
E F
G H
~n
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 25/72
Interface Conditions
For tracked segments representing material interfaces, pressure equilibrium as well as velocity continuity
conditions across interfaces are fulfilled by 1. Devise of the wave-propagation method
2. Choice of Riemann solver used in the method
Stability Issues
Choose time step ∆t based on uniform grid mesh size
∆x , ∆y as
∆t max p,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 smoothing operator (such as, h -box approach of Berger et al. ) locally for cell averages in irregular cells
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 27/72
Shock-Bubble Interaction
Shock-Bubble Interaction
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72
Shock-Bubble Interaction
Shock-Bubble Interaction
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72
Shock-Bubble Interaction
Shock-Bubble Interaction
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72
Shock-Bubble Interaction
Shock-Bubble Interaction
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72
Shock-Bubble Interaction
Shock-Bubble Interaction
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72
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) V s V R V T V ui V uf V di V df 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
V s ( V R , V T ) Incident (refracted, transmitted) shock speed t ∈ [0, 250]µ s ( t ∈ [0, 202]µ s, t ∈ [202, 250]µ s ) V ui ( V uf ) Initial (final) upstream bubble wall speed t ∈ [0, 400]µ s ( t ∈ [400, 1000]µ s)
V di ( V df ) Initial (final) downstream bubble wall speed t ∈ [200, 400]µ s ( t ∈ [400, 1000]µ s)
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 30/72
Underwater Explosions
Numerical schlieren images for density
Underwater Explosions (cont.)
Approximate locations of interfaces
time=0.2
air
water
time=0.4
air
water
time=0.8
air
water
time=1.2
air
water
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 32/72
Falling Rigid Object in Water Tank
−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
Falling Rigid Object in Water Tank
−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
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 33/72
Falling Rigid Object in Water Tank
−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
Falling Rigid Object in Water Tank
−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
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 33/72
Generalized Lagrangian Model
Introduce transformation (t, x, y) ↔ (τ, ξ, η) via
dt dx dy
=
1 0 0
x τ x ξ x η y τ y ξ y η
dτ dξ dη
or
dτ dξ dη
=
1 0 0 ξ t ξ x ξ y η t η x η y
dt dx dy
Basic grid-metric relations:
1 0 0 ξ t ξ x ξ y η t η x η y
=
1 0 0
x τ x ξ x η y τ y ξ y η
−1
= 1 J
x ξ y η − x η y ξ 0 0
−x τ y η + y τ x η y η −x η x τ y ξ − y τ x ξ −y ξ x ξ
G. Lagrangian Model (cont.)
Homogeneous two-phase model in N d generalized coord.:
∂
∂τ (α 1 ρ 1 J) +
N d
X
j=1
∂
∂ξ j (α 1 ρ 1 JU j ) = 0,
∂
∂τ (α 2 ρ 2 J) +
N d
X
j=1
∂
∂ξ j (α 2 ρ 2 JU j ) = 0,
∂
∂τ (ρJu i ) +
N d
X
j=1
∂
∂ξ j J
ρu i U j + p ∂ξ j
∂x i
= 0 for i = 1, 2, . . . , N d ,
∂
∂τ (JE) +
N d
X
j=1
∂
∂ξ j J
EU j + pU j − p ∂ξ j
∂t
= 0,
∂α 2
∂τ +
N d
X
j=1
U j ∂α 2
∂ξ j = 0, U j = ∂ t ξ j +
N d
X
i=1
u i ∂ x i ξ j
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 35/72
G. Lagrangian Model (cont.)
Continuity on mixed derivatives of grid coordinates gives geometrical conservation laws
∂
∂τ
x ξ
y ξ x η y η
+ ∂
∂ξ
−x τ
−y τ 0 0
+ ∂
∂η
0 0
−x τ
−y τ
= 0
with (x τ , y τ ) to be specified as, for example,
Eulerian case: (x τ , y τ ) = ~0
G. Lagrangian Model (cont.)
General 1 -parameter case: (x τ , y τ ) = h(u, v) , h ∈ [0, 1]
At given time instance, h can be chosen based on
Grid-angle preserving condition (Hui et al. JCP 1999)
∂
∂τ cos −1 ∇ξ
|∇ξ| · ∇η
|∇η|
= ∂
∂τ cos −1
−y η x η − y ξ x ξ q y ξ 2 + y η 2 q
x 2 ξ + x 2 η
= · · ·
= Ah ξ + Bh η + Ch = 0 (1st order PDE )
with
A = q
x 2 η + y η 2 (vx ξ − uy ξ ) , B = q
x 2 ξ + y ξ 2 (uy η − vx η ) C = q
x 2 ξ + y ξ 2 (u η y η − v η x η ) − q
x 2 η + y η 2 (u ξ y ξ − v ξ x ξ )
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 37/72
G. Lagrangian Model (cont.)
General 1 -parameter case: (x τ , y τ ) = h(u, v) , h ∈ [0, 1]
Or alternatively, based on
Mesh-area preserving condition
∂J
∂τ = ∂
∂τ (x ξ y η − x η y ξ )
= x ξτ y η + x ξ y ητ − x ητ y ξ − x η y ξτ
= · · ·
= Ah ξ + Bh η + Ch = 0 (1st order PDE )
with
G. Lagrangian Model (cont.)
Numerics: h - or ˜h -equation constraint geometrical laws
∂
∂τ
x ξ y ξ x η y η
− ∂
∂ξ
hu hv 0 0
− ∂
∂η
0 0 hu hv
= 0
Usability: Mesh-area evolution equation
∂J
∂τ − ∂
∂ξ [h (uy η − vx η )] − ∂
∂η [h (vx ξ − uy ξ )] = 0
Initial & boundary conditions for h - or ˜h -equation ?
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 38/72
G. Lagrangian Model (cont.)
In summary, with (x τ , y τ ) = h 0 (u, v) & EOS, model system for homogeneous two-phase flow reads
∂
∂τ
Jα 1 ρ 1
Jα 2 ρ 2
Jρu Jρv JE
x ξ y ξ x η
+ ∂
∂ξ
Jα 1 ρ 1 U Jα 2 ρ 2 U JρuU + y η p JρvU − x η p
JEU + (y η u − x η v)p
−h 0 u
−h 0 v 0
+ ∂
∂η
Jα 1 ρ 1 V Jα 2 ρ 2 V JρuV − y ξ p JρvV + x ξ p
JEV + (x ξ v − y ξ u)p 0
0
−h u
= 0
G. Lagrangian Model (cont.)
Under thermodyn. stability condition, our multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
Model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates
Grid system is a time-varying grid
Extension of the model to general non-barotropic
multifluid flow can be made in an analogous manner
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 40/72
Flux-based Wave Decomposition
In 2 D, equations to be solved takes the form
∂q
∂τ + f 1
∂
∂ξ , q, ∇~ ξ
+ f 2
∂
∂η , q, ∇~ ξ
= ˜ ψ
A simple dimensional-splitting approach based on f -wave formulation of LeVeque et al. is used
Solve one-dimensional generalized Riemann problem (defined below) at each cell interfaces Use resulting jumps of fluxes (decomposed into
each wave family) of Riemann solution to update cell
averages
F-Waves Decomposition (cont.)
Employ finite volume formulation of numerical solution
Q n ij ≈ 1
∆ξ∆η Z
C ij
q(ξ, η, τ n ) dA
that gives approximate value of cell average of solution q over cell C ij = [ξ i , ξ i+1 ] × [η j , η j+1 ] at time τ n
−1 −0.5 0 0.5 1
−1.5
−1
−0.5 0
−1 −0.5 0 0.5
−2
−1.8
−1.6
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
x
y
ξ
η
computational grid physical grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 42/72
Generalized Riemann Problem
Generalized Riemann problem of our multifluid model at cell interface ξ i−1/2 consists of the equation
∂q
∂τ + F i− 1
2 ,j
∂ ξ , q, ∇~ ξ
= 0 together with flux function
F i− 1
2 ,j =
f i−1,j
∂ ξ , q, ∇~ ξ
for ξ < ξ i−1/2 f ij
∂ ξ , q, ∇~ ξ
for ξ > ξ i−1/2
and piecewise constant initial data
Generalized Riemann Problem
Generalized Riemann problem at time τ = 0
q τ + f i−1,j
∂ ξ , q, ∇~ ξ
= 0 q τ + f i,j
∂ ξ , q, ∇~ ξ
= 0 ξ τ
Q n i−1,j Q n ij
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 44/72
Generalized Riemann Problem
Exact generalized Riemann solution: basic structure
τ
n Q n
Generalized Riemann Problem
Shock-only approximate Riemann solution: basic structure
q τ + f i−1,j
∂ ξ , q, ∇~ ξ
= 0 q τ + f i,j
∂ ξ , q, ∇~ ξ
= 0 ξ
τ
Q n i−1,j Q n ij
λ 1 λ 2
λ 3 q mL − q mL +
q mR
Z 1 = f L (q mL − ) − f L (Q n i−1,j ) Z 2 = f R (q mR ) − f R (q mL + )
Z 3 = f R (Q n ij ) − f R (q mR )
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 44/72
F-Waves Decomposition (cont.)
Basic steps of a dimensional-splitting scheme ξ -sweeps: solve
∂q
∂τ + f 1
∂
∂ξ , q, ∇~ ξ
= 0
updating Q n ij to Q ∗ i,j η -sweeps: solve
∂q
∂τ + f 2
∂
∂η , q, ∇~ ξ
= 0
F-Waves Decomposition (cont.)
That is to say,
ξ -sweeps: we use
Q ∗ ij = Q n ij − ∆τ
∆ξ
F i+ − 1
2 ,j − F i− + 1
2 ,j
− ∆τ
∆ξ ˜ Z i+ 1
2 ,j − ˜ Z i− 1
2 ,j
with Z ˜ i− 1
2 ,j = 1 2
m w
X
p=1
sign
λ p i− 1
2 ,j
1 − ∆τ
∆ξ λ
p
i− 1 2 ,j
Z ˜ p i− 1
2 ,j
η -sweeps: we use
Q n+1 ij = Q ∗ ij − ∆τ
∆η
G i,j+ − 1
2
− G +
i,j− 1 2
− ∆τ
∆η ˜ Z i,j+ 1
2 − ˜ Z i,j− 1
2
with Z ˜ i,j− 1
2 = 1 2
m w
X
p=1
sign λ p
i,j− 1 2
1 − ∆τ
∆η λ
p i,j− 1 2
Z ˜ p
i,j− 1 2
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 46/72
F-Waves Decomposition (cont.)
Flux-based wave decomposition
f i,j − f i−1,j =
m w
X
p=1
Z i−1/2 p =
m w
X
p=1
λ p i−1/2 W i−1/2 p
Some care should be taken on the limited jump of fluxes W ˜ p , for p = 2 (contact wave), in particular to ensure
correct pressure equilibrium across material interfaces
MUSCL-type (slope limited) high resolution extension is
not simple as one might think of for multifluid problems
Splitting of discontinuous fluxes at cell interfaces:
Lax’s Riemann Problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 0.5 1 1.5 2
x
ρ
Exact h 0 =0 h 0 =0.99
0 0.5 1
0 0.5 1 1.5 2
x
u
0 0.5 1
0 1 2 3 4
x
p
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 48/72
Lax’s Riemann Problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
0.05 0.1 0.15
time
Woodward-Colella’s Problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 2 4 6 8
ρ
0 0.5 1
−10 0 10 20
0 0.5 1
0 100 200 300 400
u p
t = 0.016 t = 0.016
t = 0.016
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 50/72
Woodward-Colella’s Problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 5 10 15 20
ρ
Fine grid h 0 =0 h 0 =0.99
−5 0 5 10 15
0 200 400 600
u p
t = 0.032 t = 0.032
t = 0.032
Woodward-Colella’s Problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 2 4 6 8
ρ
0 0.5 1
−5 0 5 10 15
0 0.5 1
0 200 400 600
u p
x x
x
t = 0.038 t = 0.038
t = 0.038
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 52/72
Woodward-Colella’s Problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
0.01 0.02 0.03 0.04
time
2 D Riemann Problem
With initial 4 -shock wave pattern
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ρ u v p
=
1.5
0 0 1.5
0.532 1.206
0 0.3
0.138 1.206 1.206 0.029
0.532 0 1.206
0.3
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 54/72
2 D Riemann Problem
With initial 4 -shock wave pattern Lagrange-like result
Occurrence of simple Mach reflection
0.6 0.8
0.6 0.8
0.6 0.8
Density Pressure Physical grid
2 D Riemann Problem
With initial 4 -shock wave pattern Eulerian result
Poor resolution around simple Mach reflection
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
Density Pressure Physical grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 56/72
More Examples
Two-dimensional case
Radially symmetric problem Underwater explosion
Shock-bubble interaction Helium bubble case
Refrigerant bubble case Three-dimensional case
Underwater explosion Shock-bubble interaction
Helium bubble case
Radially Symmetric Problem
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4
Density Pressure 0.5 Physical grid
gas
liquid a) h 0 = 0.99
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4
Density Pressure 0.5 Physical grid
gas
liquid b) h 0 = 0
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 58/72
Radially Symmetric Prob. (Cont.)
0 0.1 0.2 0.3 0.4 0.5
0.2 0.4 0.6 0.8 1 1.2
one−d h 0 =0 h 0 =0.99
0 0.1 0.2 0.3 0.4 0.5
0 0.1 0.2 0.3 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.5 1 1.5
0 0.1 0.2 0.3 0.4 0.5
0.5 1 1.5 2
ρ (M g /m 3 ) ¯u (k m /s )
p (G P a ) J
gas liquid
Underwater Explosions
Numerical schlieren images h 0 = 0.9 , 800 × 500 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72
Underwater Explosions
Numerical schlieren images h 0 = 0.9 , 800 × 500 grid
Underwater Explosions
Numerical schlieren images h 0 = 0.9 , 800 × 500 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72
Underwater Explosions
Numerical schlieren images h 0 = 0.9 , 800 × 500 grid
Underwater Explosions
Numerical schlieren images h 0 = 0.9 , 800 × 500 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72
Underwater Explosions (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.9
−1
−0.5 0 0.5
time=0ms
Underwater Explosions (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.9
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
−1
−0.5 0 0.5
time=0.2ms
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 61/72
Underwater Explosions (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.9
−1
−0.5 0 0.5
time=0.4ms
Underwater Explosions (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.9
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
−1
−0.5 0 0.5
time=0.8ms
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 61/72
Underwater Explosions (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.9
−1
−0.5 0 0.5
time=1.2ms
Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72
Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72
Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.02
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 63/72
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.08
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.16
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 63/72
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.35
Shock-Bubble (Refrigerant)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72
Shock-Bubble (Refrigerant)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Shock-Bubble (Refrigerant)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72
Shock-Bubble (Refrigerant)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Shock-Bubble (Refrigerant)
Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72
Shock-Bubble (R22) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0
Shock-Bubble (R22) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.02
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 65/72
Shock-Bubble (R22) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.08
Shock-Bubble (R22) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.16
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 65/72
Shock-Bubble (R22) (Cont.)
Grid system (coarsen by factor 5 ) with h 0 = 0.5
−0.1
−0.05 0 0.05 0.1 0.15 0.2
time=0.35
Underwater Explosions
Numerical schlieren images h 0 = 0.6 , 100 3 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72
Underwater Explosions
Numerical schlieren images h 0 = 0.6 , 100 3 grid
Underwater Explosions
Numerical schlieren images h 0 = 0.6 , 100 3 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72
Underwater Explosions
Numerical schlieren images h 0 = 0.6 , 100 3 grid
Underwater Explosions
Numerical schlieren images h 0 = 0.6 , 100 3 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72
3 D Underwater Explosions (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
time = 0
3 D Underwater Explosions (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
−1
−0.5 0
0.5 1
−1
−0.5 0
0.5 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
time = 0.25 ms
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 67/72
3 D Underwater Explosions (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
time = 0.5 ms
3 D Underwater Explosions (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
−1
−0.5 0
0.5 1
−1
−0.5 0
0.5 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
time = 1.0 ms
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 67/72
3 D Underwater Explosions (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
time = 1.5 ms
3 D Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72
3 D Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid
3 D Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72
3 D Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid
3 D Shock-Bubble (Helium)
Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
0.3
0.4
0.5
0.6
0.7
0 0.05 0.1 0.15 0.2 0.25
z
time = 0
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.1 0 0.2
0 0.05 0.1 0.15 0.2 0.25
x y
z
time = 0.02
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 69/72
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
0.3
0.4
0.5
0.6
0.7
0 0.05 0.1 0.15 0.2 0.25
z
time = 0.08
Shock-Bubble (Helium) (Cont.)
Grid system (coarsen by factor 2 ) with h 0 = 0.6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.1 0 0.2
0 0.05 0.1 0.15 0.2 0.25
x y
z
time = 0.16
Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 69/72