Fluid-mixture type algorithm for
compressible multifluid flows in
generalized curvilinear grids
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Taiwan
Overview
Mathematical model for homogeneous multifluid flow Compressible Euler eqs. in generalized coordinates Grid-movement conditions for moving grid system Mixture equations of state
Transport eqs. for multifluid problems of concerns Finite volume numerical method
Godunov-type f -wave formulation of LeVeque et al.
Numerical examples
Underwater explosions, shock-bubble, · · ·
Future direction
Motivations
Some basic facts
Lagrangian method can resolve material or slip lines sharply if there is not too much grid tangling
Generalized curvilinear grid is often superior to
Cartesian grid when they are employed in numerical
methods for complex fixed or moving geometries
Motivations
Some basic facts
Lagrangian method can resolve material or slip lines sharply if there is not too much grid tangling
Generalized curvilinear grid is often superior to
Cartesian grid when they are employed in numerical methods for complex fixed or moving geometries
Some examples done by Cartesian-based method Falling liquid drop problem
Shock-bubble interaction
Flying projectile & ocean surface
Falling rigid object in water tank
Motivations
Some basic facts
Lagrangian method can resolve material or slip lines sharply if there is not too much grid tangling
Generalized curvilinear grid is often superior to
Cartesian grid when they are employed in numerical methods for complex fixed or moving geometries
Some examples done by Cartesian-based method Falling liquid drop problem
Shock-bubble interaction
Flying projectile & ocean surface Falling rigid object in water tank
Search for more robust method (work present here is
preliminary)
Falling Liquid Drop Problem
Interface capturing with gravity
100 200 300 400 500 600 700 800 900 1000
air
t = 0
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900 1000
air
t = 0.2s
Falling Liquid Drop Problem
Interface diffused badily
100 200 300 400 500 600 700 800 900
air
t = 0.4s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.6s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.8s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 1s
Shock-Bubble Interaction
Volume tracking for material interface
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Small moving irregular cells: stability & accuracy
Flying Projectile & Ocean Surface
Moving boundary tracking & interface capturing
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 0
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 12ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 24ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 36ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 48ms
Flying Projectile & Ocean Surface
Small moving irregular cells: stability & accuracy
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water
t = 60ms
Falling Rigid Object in Water Tank
Moving boundary tracking & interface capturing
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
Density Pressure Volume fraction
air
water
t = 0ms
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
Density Pressure Volume fraction
air
water
t = 1ms
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
Density Pressure Volume fraction
air
water
t = 2ms
Falling Rigid Object in Water Tank
Small moving irregular cells: stability & accuracy
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
−2 0 2
−3
−2
−1 0 1 2 3
Density Pressure Volume fraction
air
water
t = 3ms
Euler Eqs. in Generalized Coord.
With gravity effect included, for example, 2 D compressible Euler eqs. in Cartesian coordinates take
∂q
∂t + ∂f (q)
∂x + ∂g(q)
∂y = ψ(q)
where
q =
ρ ρu ρv E
, f (q) =
ρu ρu 2 + p
ρuv Eu + pu
, g(q) =
ρv ρuv ρv 2 + p Ev + pv
, ψ =
0 0 ρg ρgv
ρ : density , (u, v) : vector of particle velocity
p : pressure , E = ρ[e + (u 2 + v 2 )/2] : total energy
e(ρ, p) : internal energy , ψ : gravitational source term
Euler in General. Coord. (Cont.)
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 η −y ξ x τ y ξ − y τ x ξ −x η x ξ
J = x ξ y η − x η y ξ : grid Jacobian
Euler in General. Coord. (Cont.)
With these notations, Euler eqs. in generalized coord. are
∂ ˜ q
∂τ + ∂ ˜ f
∂ξ + ∂˜ g
∂η = ˜ ψ
where
˜
q = J
ρ ρu ρv E
, ˜ f = J
ρU
ρuU + ξ x p ρvU + ξ y p EU + pU − ξ t p
, ˜ g = J
ρV
ρuV + η x p ρvV + η y p EV + pV − η t p
, ˜ ψ = J
0 0 ρg ρgv
with contravariant velocities U & V defined by
U = ξ t + ξ x u + ξ y v & V = η t + η x u + η y v
Grid Movement Conditions
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
Lagrangian case: (x τ , y τ ) = (u, v)
Lagrangian-like case: (x τ , y τ ) = h 0 (u, v) or (h 0 u, k 0 v)
h 0 ∈ [0, 1] & k 0 ∈ [0, 1]
Grid Movement Conditions (Cont.)
General 1 -parameter case: (x τ , y τ ) = h(u, v) Mesh-area preserving condition
∂J
∂τ = ∂
∂τ (x ξ y η − x η y ξ )
= x ξτ y η + x ξ y ητ − x ητ y ξ − x η y ξτ
= · · ·
= Ah ξ + Bh η + Ch = 0 (1st order PDE for h ∈ [0, 1])
with
A = uy η − vx η , B = vx ξ − uy ξ C = u ξ y η + v η x ξ − u η y ξ − v ξ x η
Initial & boundary conditions for h -equation ?
Grid Movement Conditions (Cont.)
General 1 -parameter case: (x τ , y τ ) = h(u, v)
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 ξ )
Initial & boundary conditions for h -equation ?
Grid Movement Conditions (Cont.)
2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions
1. Grid-angle preserving
2. Specialized grid-material line matching (see next)
Grid Movement Conditions (Cont.)
2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions
1. Grid-angle preserving
2. Specialized grid-material line matching (see next)
Good results are shown for steady-state problems
Grid Movement Conditions (Cont.)
2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions
1. Grid-angle preserving
2. Specialized grid-material line matching (see next)
Good results are shown for steady-state problems
Little results for time-dependent problems with rapid
transient solution structures
Grid Movement Conditions (Cont.)
2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions
1. Grid-angle preserving
2. Specialized grid-material line matching (see next) Good results are shown for steady-state problems Little results for time-dependent problems with rapid transient solution structures
Other 2 -parameter case: (x τ , y τ ) = (hu, k v)
Novel imposed conditions for h ∈ [0, 1] & k ∈ [0, 1] ?
Grid Movement Conditions (Cont.)
2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions
1. Grid-angle preserving
2. Specialized grid-material line matching (see next) Good results are shown for steady-state problems Little results for time-dependent problems with rapid transient solution structures
Other 2 -parameter case: (x τ , y τ ) = (hu, k v)
Novel imposed conditions for h ∈ [0, 1] & k ∈ [0, 1] ? Roadmap of current work:
(x τ , y τ ) = h 0 (u, v) → (x τ , y τ ) = h(u, v) → · · ·
Single-Fluid Model
With (x τ , y τ ) = h 0 (u, v) , our model system for single-phase flow reads
∂
∂τ
Jρ Jρu Jρv JE
x ξ y ξ x η y η
+ ∂
∂ξ
JρU
JρuU + y η p JρvU − x η p
JE U + (y η u − x η v)p
−h 0 u
−h 0 v 0 0
+ ∂
∂η
JρV
JρuV − y ξ p Jρv V + x ξ p
JEV + (x ξ v − y ξ u)p 0
0
−h 0 u
−h 0 v
= ˜ ψ
where U = (1 − h 0 )(y η u − x η v) & V = (1 − h 0 )(x ξ v − y ξ u)
Single-Fluid Model: Remarks
Hyperbolicity (under thermodyn. stability cond.)
In Cartesian coordinates, model is hyperbolic
Single-Fluid Model: Remarks
Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic
In generalized-moving coord., model is hyperbolic
when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
Single-Fluid Model: Remarks
Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic
In generalized-moving coord., model is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1 Canonical form
In Cartesian coordinates
∂q
∂t + ∂f (q)
∂x + ∂g(q)
∂y = ψ(q) In generalized coordinates
∂q
∂τ + ∂f (q, Ξ)
∂ξ + ∂g(q, Ξ)
∂η = ψ(q), Ξ : grid metrics
Single-Fluid Model: Remarks
Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic
In generalized-moving coord., model is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1 Canonical form
In Cartesian coordinates
∂q
∂t + ∂f (q)
∂x + ∂g(q)
∂y = ψ(q)
In generalized coordinates : spatially varying fluxes
∂q
∂τ + ∂ f (q, Ξ)
∂ξ + ∂ g(q, Ξ)
∂η = ψ(q), Ξ : grid metrics
Extension to Multifluid
Assume homogeneous ( 1 -pressure & 1 -velocity) flow;
i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι
gas
gas gas
gas
gas gas
gas
gas gas
gas
liquid
Extension to Multifluid
Assume homogeneous ( 1 -pressure & 1 -velocity) flow;
i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι Mathematical model: Fluid-mixture type
Use basic conservation (or balance) laws for single
& multicomponent fluid mixtures
Introduce additional transport equations for
problem-dependent material quantities near
numerically diffused interfaces, yielding direct
computation of pressure from EOS
Extension to Multifluid
Assume homogeneous ( 1 -pressure & 1 -velocity) flow;
i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι Mathematical model: Fluid-mixture type
Use basic conservation (or balance) laws for single
& multicomponent fluid mixtures
Introduce additional transport equations for problem-dependent material quantities near numerically diffused interfaces, yielding direct computation of pressure from EOS
Sample examples
Barotropic 2 -phase flow
Hybrid barotropic & non-barotropic 2 -phase flow
Barotropic 2-Phase Flow
Equations of state
Fluid component 1 & 2 : Tait EOS
p(ρ) = (p 0ι + B ι )
ρ ρ 0ι
γ ι
− B ι , ι = 1, 2
Barotropic 2-Phase Flow
Equations of state
Fluid component 1 & 2 : Tait EOS
p(ρ) = (p 0ι + B ι )
ρ ρ 0ι
γ ι
− B ι , ι = 1, 2
Mixture pressure law (Shyue, JCP 2004)
p(ρ, ρe) =
(p 0ι + B ι )
ρ ρ 0ι
γ ι
− B ι if α = 0 or 1 (γ − 1)
ρe + ρB ρ 0
− γB if α ∈ (0, 1)
Here α denotes volume fraction of one chosen fluid component
Barotropic 2-Phase Flow
Equations of state
Fluid component 1 & 2 : Tait EOS
p(ρ) = (p 0ι + B ι )
ρ ρ 0ι
γ ι
− B ι , ι = 1, 2
Mixture pressure law (Shyue, JCP 2004)
p(ρ, ρe) =
(p 0ι + B ι )
ρ ρ 0ι
γ ι
− B ι if α = 0 or 1 (γ − 1)
ρe + ρB ρ 0
− γB if α ∈ (0, 1)
variant form of p(ρ, S) = A(S) (p 0 + B) ρ ρ 0
γ
− B
A(S) = e [(S−S 0 )/C V ] , S , C V : specific entropy & heat at constant volume
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
Above equations are derived from energy equation
& make use of homogeneous equilibrium flow
assumption together with mass conservation law
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
If we ignore JBρ/ρ 0 term, they are essentially
equations proposed by Saurel & Abgrall (SISC
1999), but are written in generalized coord.
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
α -based equations
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0, with z =
2
X
ι=1
α ι z ι , z = 1
γ − 1 & γB γ − 1
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
α -based equations (Allaire et al. , JCP 2002)
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0 with z =
2
X
ι=1
α ι z ι , z = 1
γ − 1 & γB γ − 1
∂ (Jρ 1 α) + ∂
(Jρ 1 αU ) + ∂
(Jρ 1 αV ) = 0 with z = B
ρ
Barotropic 2-Phase Flow (Cont.)
Transport equations for material quantities γ , B , & ρ 0
γ -based equations
∂
∂τ
1 γ − 1
+ U ∂
∂ξ
1 γ − 1
+ V ∂
∂η
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+ U ∂
∂ξ
γB γ − 1
+ V ∂
∂η
γB γ − 1
= 0
∂
∂τ
J B
ρ 0 ρ
+ ∂
∂ξ
J B
ρ 0 ρU
+ ∂
∂η
J B
ρ 0 ρV
= 0
α -based equations (Kapila et al. , Phys. Fluid 2001)
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = α 1 α 2 ρ 1 c 2 1 − ρ 2 c 2 2 P 2
k=1 α k ρ k c 2 k
!
∇ · ~u
· · · will not be discussed here
Barotropic & Non-Barotropic Flow
Equations of state
Fluid component 1 : Tait EOS
p(ρ) = (p 0 + B) ρ ρ 0
γ
− B
Fluid component 2 : Noble-Abel EOS
p(ρ, ρe) = γ − 1 1 − bρ
ρe (0 ≤ b ≤ 1/ρ)
Barotropic & Non-Barotropic Flow
Equations of state
Fluid component 1 : Tait EOS
p(ρ) = (p 0 + B) ρ ρ 0
γ
− B
Fluid component 2 : Noble-Abel EOS
p(ρ, ρe) = γ − 1 1 − bρ
ρe (0 ≤ b ≤ 1/ρ)
Mixture pressure law (Shyue, Shock Waves 2006)
p(ρ, ρe) =
(p 0 + B) ρ ρ 0
γ
− B if α = 1 (fluid 1)
γ − 1 1 − bρ
(ρe − B) − B if α 6= 1
Baro. & Non-Baro. Flow (Cont.)
Equations of state
Fluid component 1 : Tait EOS
p(V ) = A(S 0 ) (p 0 + B) V 0 V
γ
− B, V = 1/ρ
Fluid component 2 : Noble-Abel EOS
p(V, S) = A(S)p 0 V 0 − b V − b
γ
Mixture pressure law
p(V, S) = A(S) (p 0 + B) V 0 − b V − b
γ
− B
Baro. & Non-Baro. Flow (Cont.)
Equations of state
Fluid component 1 : Tait EOS
p(V ) = A(S 0 ) (p 0 + B) V 0 V
γ
− B, V = 1/ρ
Fluid component 2 : Noble-Abel EOS
p(V, S) = A(S)p 0 V 0 − b V − b
γ
Mixture pressure law
p(V, S) = A(S) (p 0 + B) V 0 − b V − b
γ
− B variant form of p(ρ, ρe) = γ − 1
1 − bρ
(ρe − B) − B
Baro. & Non-Baro. Flow (Cont.)
Transport equations for material quantities γ , b , & B α -based equations
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0
∂
∂τ (Jρ 1 α) + ∂
∂ξ (Jρ 1 αU ) + ∂
∂η (Jρ 1 αV ) = 0 with z = P 2
ι=1 α ι z ι , z = γ−1 1 , γ−1 bρ , & γ−bρ γ−1 B
Baro. & Non-Baro. Flow (Cont.)
Transport equations for material quantities γ , b , & B α -based equations
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0
∂
∂τ (Jρ 1 α) + ∂
∂ξ (Jρ 1 αU ) + ∂
∂η (Jρ 1 αV ) = 0 with z = P 2
ι=1 α ι z ι , z = γ−1 1 , γ−1 bρ , & γ−bρ γ−1 B
Note: 1 − bρ
γ − 1 p + γ − bρ
γ − 1 B = ρe =
2
X
ι=1
α ι ρ ι e ι
=
2
X
ι=1
α ι 1 − b ι ρ ι
γ ι − 1 p ι + γ ι − b ι ρ ι γ ι − 1 B ι
Baro. & Non-Baro. Flow (Cont.)
Transport equations for material quantities γ , b , & B α -based equations
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0
∂
∂τ (Jρ 1 α) + ∂
∂ξ (Jρ 1 αU ) + ∂
∂η (Jρ 1 αV ) = 0 with z = P 2
ι=1 α ι z ι , z = γ−1 1 , γ−1 bρ , & γ−bρ γ−1 B
Note: 1 − bρ
γ − 1 p + γ − bρ
γ − 1 B = ρe =
2
X
ι=1
α ι ρ ι e ι
=
2
X
ι=1
α ι 1 − b ι ρ ι
γ ι − 1 p ι + γ ι − b ι ρ ι γ ι − 1 B ι
Multifluid Model
With (x τ , y τ ) = h 0 (u, v) & sample EOS described above, our α -based model for multifluid flow is
∂
∂τ
Jρ Jρu Jρv JE
x ξ y ξ x η y η Jρ 1 α
+ ∂
∂ξ
JρU
JρuU + y η p JρvU − x η p
JEU + (y η u − x η v)p
−h 0 u
−h 0 v 0 0 Jρ 1 αU
+ ∂
∂η
JρV
JρuV − y ξ p JρvV + x ξ p
JEV + (x ξ v − y ξ u)p 0
0
−h 0 u
−h 0 v Jρ 1 αV
= ˜ ψ
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0
Multifluid Model (Cont.)
For convenience, our multifluid model is written into
∂q
∂τ + f ∂
∂ξ , q, Ξ
+ g ∂
∂η , q, Ξ
= ˜ ψ with
q = [Jρ, Jρu, Jρv, JE, x ξ , y ξ , x η , y η , Jρ 1 α, α] T f = ∂
∂ξ (JρU ), ∂
∂ξ (JρuU + y η p), ∂
∂ξ (JρvU − x η p), ∂
∂ξ (JEU + (y η u − x η v)p),
∂
∂ξ (−h 0 u), ∂
∂ξ (−h 0 v), 0, 0, ∂
∂ξ (Jρ 1 αU ), U ∂α
∂ξ
T
g = ∂
∂η (JρV ), ∂
∂η (JρuV − y ξ p), ∂
∂η (JρvV + x ξ p), ∂
∂η (JEV + (x ξ v − y ξ u)p), 0, 0, ∂
∂η (−h 0 u), ∂
∂η (−h 0 v), ∂
∂η (Jρ 1 αV ), V ∂α
∂η
T
Multifluid model: Remarks
As before, under thermodyn. stability condition, our
multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
Our model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates Our grid system is a time-varying grid
Extension of the model to general non-barotropic
multifluid flow can be made in an analogous manner
Multifluid model: Remarks
As before, under thermodyn. stability condition, our
multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
Our model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates Our grid system is a time-varying grid
Extension of the model to general non-barotropic
multifluid flow can be made in an analogous manner
Numerical approximation ?
Numerical Approximation
Equations to be solved are
∂q
∂τ + f ∂
∂ξ , q, Ξ
+ g ∂
∂η , 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
Introduce limited jumps of fluxes to achieve high
resolution
Numerical Approximation (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
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
q(ξ, 0) =
Q n i−1,j for ξ < ξ i−1/2
Q n ij for ξ > ξ i−1/2
General. Riemann Problem (Cont.)
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
General. Riemann Problem (Cont.)
Exact generalized 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
General. Riemann Problem (Cont.)
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
W 1 = f L (q mL − ) − f L (Q n i−1,j ) W 2 = f R (q mR ) − f R (q mL + )
W 3 = f R (Q n ij ) − f R (q mR )
Numerical Approximation (Cont.)
Basic steps of a dimensional-splitting scheme ξ -sweeps: solve
∂q
∂τ + f ∂
∂ξ , q, Ξ
= 0
updating Q n ij to Q ∗ i,j η -sweeps: solve
∂q
∂τ + g ∂
∂η , q, Ξ
= 0
updating Q ∗ ij to Q n+1 i,j
Numerical Approximation (Cont.)
That is to say,
ξ -sweeps: we use
Q ∗ ij = Q n ij − ∆τ
∆ξ
F i+ − 1
2 ,j − F i− + 1
2 ,j
− ∆τ
∆ξ ˜ F i+ 1
2 ,j − ˜ F i− 1
2 ,j
with F ˜ i− 1
2 ,j = 1 2
m w
X
p=1
sign
λ p i− 1
2 ,j
1 − ∆τ
∆ξ λ
p
i− 1 2 ,j
W ˜ p i− 1
2 ,j
η -sweeps: we use
Q n+1 ij = Q ∗ ij − ∆τ
∆η
G i,j+ − 1
2
− G +
i,j− 1 2
− ∆τ
∆η ˜ G i,j+ 1
2 − ˜ G i,j− 1
2
with G ˜ i,j− 1
2 = 1 2
m w
X
p=1
sign λ p
i,j− 1 2
1 − ∆τ
∆η λ
p i,j− 1 2
W ˜ p
i,j− 1 2
Numerical Approximation (Cont.)
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 First order or high resolution method for geometric
conservation laws ? Their effect to the grid uniformity,
· · ·
Numerical Examples
2 D Riemann problem Underwater explosion Shock-bubble interaction
Helium bubble case
Refrigerant bubble case
2D Riemann Problem
Initial condition for 4 -shock wave pattern
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
