Towards
a Unified Coordinate Method for
Compressible Multifluid Flows
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Taiwan
Outline
Motivation
Inviscid compressible flow model
Generalized coordinates Euler equations Grid movement conditions
Equations of state
Transport eqs. for multifluid problems Finite volume numerical method
Godunov-type f -wave formulation of LeVeque et al.
Numerical examples
Future direction
Motivation
Some basic facts
Lagrangian method resolve material or slip lines sharply if no grid tangling
Generalized curvilinear grid is often superior to
Cartesian when employed in numerical methods for
complex fixed or moving geometries
Motivation
Some basic facts
Lagrangian method resolve material or slip lines sharply if no grid tangling
Generalized curvilinear grid is often superior to
Cartesian when 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 Aluminum-plate problem
Falling rigid object in water tank
Motivation
Some basic facts
Lagrangian method resolve material or slip lines sharply if no grid tangling
Generalized curvilinear grid is often superior to
Cartesian when 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 Aluminum-plate problem
Falling rigid object in water tank
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 Aluminum-plate problem
Vacuum-Al interface tracking
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
Flying Aluminum-plate problem
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
Flying Aluminum-plate problem
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
Flying Aluminum-plate problem
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
Flying Aluminum-plate problem
Small moving irregular cells: stability & accuracy
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
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
2 2
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 E V + pV − η t p
, ˜ ψ = J
0 0 ρg ρgv
with contravariant velocities U & V defined by
Euler in General. Coord. (Cont.)
Model system in quasi-linear form
∂ ˜ q
∂τ + A ∂ ˜ q
∂ξ + B ∂ ˜ q
∂η = ˜ ψ
A = ∂ ˜ f
∂ ˜ q =
ξ t ξ x ξ y 0
ξ x p ρ − uU ξ x u(1 − p E ) + U ξ y u − ξ x vp E ξ x p E ξ y p ρ − vU ξ x v − ξ y up E ξ y v(1 − p E ) + U ξ y p E (p ρ − H)U ξ x H − uU p E ξ y H − vU p E U + p E U
B = ∂˜ g
∂ ˜ q =
η t η x η y 0
η x p ρ − uV η x u(1 − p E ) + V η y u − η x vp E η x p E η y p ρ − vV η x v − η y up E η y v(1 − p E ) + V η y p E (p ρ − H)V η x H − uV p E η y H − vV p E V + p E V
with H = (E + p)/ρ, U = U − ξ t = ξ x u + ξ y v, V = V − η t = η x u + η y v
Euler in General. Coord. (Cont.)
Eigen-structure of matrix A is
Λ A = diag
U − c q
ξ x 2 + ξ y 2 , U, U, U + c q
ξ x 2 + ξ y 2
R A =
1 1 0 1
u − α 1 c u α 2 u + α 1 c v − α 2 c v −α 1 v + α 2 c H − U 1 c H − c 2 /p E −U 2 H + U 1 c
L A =
(p ρ + cU 1 )/2c 2 −(α 1 c + up E )/2c 2 −(α 2 c + vp E )/2c 2 p E /2c 2 1 − p ρ /c 2 up E /c 2 vp E /c 2 −p E /c 2
U 2 α 2 −α 1 0
(p ρ − cU 1 )/2c 2 (α 1 c − up E )/2c 2 (α 2 c − vp E )/2c 2 p E /2c 2
Euler in General. Coord. (Cont.)
Eigen-structure of matrix B is
Λ B = diag
V − c q
η x 2 + η y 2 , V, V, V + c q
η x 2 + η y 2
R B =
1 1 0 1
u − β 1 c u β 2 u + β 1 c v − β 2 c v −β 1 v + β 2 c H − V 1 c H − c 2 /p E −V 2 H + V 1 c
L B =
(p ρ + cV 1 )/2c 2 −(β 1 c + up E )/2c 2 −(β 2 c + vp E )/2c 2 p E /2c 2 1 − p ρ /c 2 up E /c 2 vp E /c 2 −p E /c 2
V 2 β 2 −β 1 0
(p ρ − cV 1 )/2c 2 (β 1 c − up E )/2c 2 (β 2 c − vp E )/2c 2 p E /2c 2
with (β 1 , β 2 ) = (η x , η y )/ q
η x 2 + η y 2 , V 1 = β 1 u + β 2 v, V 2 = −β 2 u + β 1 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)
Grid Movement (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 ξ )
Grid Movement (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
Grid Movement (Cont.)
To ensure h ∈ [0, 1] , transformed variable ˜h = κ(h) is used,
e.g. , Hui et al. employed κ = ln (εh|~u|) , ε normalized constant, yielding
A˜ ˜ h ξ + ˜ B˜ h η + ˜ C = 0
Grid-angle preserving case
A = ˜ q
x 2 η + y η 2 (x ξ sin θ − y ξ cos θ) , B = ˜ q
x 2 ξ + y ξ 2 (y η cos θ − x η sin θ) C = ˜ q
x 2 ξ + y ξ 2 [y η (cos θ) η − x η (sin θ) η ] − q
x 2 η + y η 2 [y ξ (cos θ) ξ − x ξ (sin θ) ξ ]
Mesh-area preserving case
A = y ˜ η cos θ − x η sin θ, B = x ˜ ξ sin θ − y ξ cos θ
C = y ˜ η (cos θ) ξ − x η (sin θ) ξ + x ξ (sin θ) η − y ξ (cos θ) η
where ~u = (u, v) = |~u|(cos θ, sin θ)
Grid Movement: Remarks
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 ?
Grid Movement: 2 Free Degrees
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: 2 Free Degrees
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: 2 Free Degrees
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: 2 Free Degrees
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: 2 Free Degrees
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) → · · ·
Novel Conditions for h & k
Mesh-area preserving case
∂J
∂τ = ∂
∂τ (x ξ y η − x η y ξ )
= x ξτ y η + x ξ y ητ − x ητ y ξ − x η y ξτ
= · · ·
= (A 1 h ξ + B 1 h η + C 1 h) + (A 2 k ξ + B 2 k η + C 2 k) = 0, yielding, for example,
A 1 h ξ + B 1 h η + C 1 h = 0 A 2 k ξ + B 2 k η + C 2 k = 0 with
A = uy , B = uy , C = u y − u y
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, Ξ)
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
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 α ι z ι , z = 1
γ − 1 & γB
γ − 1
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 − b ι ρ ι γ ι − 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, plus α-averaged material quantities
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),
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
−1.6
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
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 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 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 )
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
− ∆τ
∆ξ ˜ 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 = 1 X m w
sign
λ p 1
1 − ∆τ p
1
Z ˜ p 1
Numerical Approx.: Remarks
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:
significance ?
First order or high resolution method for geometric
conservation laws: significance to grid uniformity ?
Numerical Examples: 2D
2 D Riemann problem Underwater explosion Shock-bubble interaction
Helium bubble case
Refrigerant bubble case
2 D 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
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
