**Eulerian**

**interface-sharpening methods** **for**

**hyperbolic problems**

**Application to compressible multiphase flow**

**Application to compressible multiphase flow**

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

**Outline**

Computing monotone sharp resolution of interfaces is of fundamental importance in many practical problems of interest

Discuss a simple Eulerian interface sharpening approach (vs. Lagrangian, interface tracking, or adaptive moving mesh) for hyperbolic problems Review two PDE-based interface sharpening techniques for solving volume-fraction linear

transport equation that arises, for example, from viscous incompressibe 2-phase flow

Extend method for computing material lines or free surfaces arising from compressible multiphase flow

**Incompressible 2-phase flow: Review**

Consider unsteady, incompressible, viscous, immiscible 2-phase flow with governing equations

∇ · ~u = 0 (Continuity)

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = ∇ · τ + ρ~g + f~_{σ} (Momentum)

∂_{t}α + ~u · ∇α = 0 (Volume-fraction transport )
Material quantities in 2-phase coexistent region are often

computed by α-based weighted average as

ρ = αρ_{1} + (1 − α) ρ_{2}, ǫ = αǫ_{1} + (1 − α) ǫ_{2}, α ∈ [0, 1],
where source terms are volume-fraction dependent

τ = ǫ ∇~u + ∇~u^{T}

, f~_{σ} = −σκ∇α with κ = ∇ ·

∇α

|∇α|

**Interface sharpening techniques**

Typical interface sharpening methods for computing volume fraction in incompressible 2-phase flow include:

Algebraic based approach

CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes): Ubbink & Issa JCP 1999

THINC (Tangent of Hyperbola for INterface

Capturing): Xiao, Honma & Kono Int. J. Numer.

Meth. Fluids 2005

Improved THINC: Xiao ^{et al.}

PDE based approah

Artificial compression: Harten CPAM 1977, Olsson &

Kreiss JCP 2005

Anti-diffusion: So, Hu & Adams JCP 2011

**Artificial interface compression**

Our first interface-sharpening model concerns nonlinear artificial compression of form proposed by Olsson & Kreiss

∂_{t}α + ~u · ∇α = 1

µ∇ · [(D (∇α · ~n) − α (1 − α))~n]

where ~n = ∇α/|∇α|, D > 0, µ ≫ 1

Numerical method based on fractional step may apply
1** ^{.}** Advection step over a time step ∆t to solve

∂_{t}α + ~u · ∇α = 0 or ∂_{t}α + ∇ · (α~u) = 0
since by assumption ∇ · ~u = 0

2** ^{.}** Interface compression step towards τ-steady state

∂ α = ∇ · [(D (∇α · ~n) − α (1 − α))~n] , τ = t/µ

**Square wave passive advection**

Square-wave pluse moving with u = 1 after 4 periodic cycle

0 0.2 0.4 0.6 0.8 1

No compress With compress Exact

α

**Zalesak’s rotating disc**

Contours α = 0.5 at 4 different times within 1 period in that

~u = (1/2 − y, x − 1/2)

No compression With compression

**Vortex in cell**

Contours α = (0.05, 0.5, 0.95) at 6 different times in 1 period

~u = − sin^{2} (πx) sin (2πy), sin (2πx) sin^{2} (πy)

cos (πt/8) No compression

With compression

**Interface compression: Remarks**

Consider 1D model problem with u > 0 of form

∂_{t}α + u∂_{x}α = 1

µ∂_{x} [D (∂_{x}α · ~n) − α (1 − α)]

α(x, 0) = α_{0}(x) = 1/ (1 + exp (−x/D)) , x ∈ lR, t > 0
Exact solution for this problem is simply α(x, t) = α_{0}(x − ut)
When α(x, 0) is perturbed to α_{0}(x) + δ(x), δ(x) ≪ 1, we have

∂_{τ}α + ∂˜ _{ξ} α˜^{2}/2

= ∂_{ξ} D∂_{ξ}α˜

, α(ξ, 0) = ˜˜ α_{0}(ξ)
with ξ = x − ut, τ = t/µ, & α(ξ, τ ) = 1 − α(ξ, τ )˜ , yielding
steady state solution α(ξ, τ ) = ˜α_{0} (ξ + ξ_{0}) as τ → ∞ for
some suitably chosen shift ξ_{0}, see Sattinger (1976)

**Interface compression: Remarks**

If perturbation is zero mass, * ^{i.e.}*, R ∞

−∞ δ(ξ, 0)dξ = 0 we have
true solution with ξ_{0} = 0, see Goodman (1986)

When model is solved by a conservative method, truncation errors will be of zero mass, yielding convergence of

numerical solution to exact one in time we want

In multi-D case, let K_{α} = D∇α · ~n − α (1 − α). We solve

∂_{τ}α = ∇ · [(D (∇α · ~n) − α (1 − α))~n] = K_{α}∇ · ~n + ~n · ∇K_{α}
yielding τ-steady state solution as µ → ∞, when K_{α} = 0 &

1D profile in coordinate normal to interface

When µ finite, K_{α}∇ · ~n + ~n · ∇K_{α} 6= 0, strength & accuracy of
curvature ∇ · ~n plays important role in interface resolution

**Interface compression runs**

Methods used here are very elementary, * ^{i.e.}*,
1

**Use Clawpack for advection in Step 1**

^{.}2** ^{.}** Use simple forward Euler in time, second order in space
for interface compression in Step 2

Diffusion coefficient D = ε min_{∀i} ∆x_{i}
Time step ∆τ

∆τ ≤ 1 2D

Xd

i=1

∆x^{2}_{i}

Stopping criterion: simple 1-norm error measure

**Extension to compressible flow**

Shukla, Pantano & Freund (JCP 2010) proposed extension of interface-compression method for incompressible flow to compressible flow governed by reduced 2-phase model as

∂_{t} (α_{1}ρ_{1}) + ∇ · (α_{1}ρ_{1}~u) = 0

∂_{t} (α_{2}ρ_{2}) + ∇ · (α_{2}ρ_{2}~u) = 0

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_{t}(ρE) + ∇ · (ρE~u + p~u) = 0

∂_{t}α_{1} + ~u · ∇α_{1} = 1

µ~n · ∇ (D∇α_{1} · ~n − α_{1} (1 − α_{1}))

∂_{t}ρ + ∇ · (ρ~u) = 1

µH(α_{1})~n · (∇ (D∇ρ · ~n) − (1 − 2α_{1}) ∇ρ)
Mixture pressure is computed based on isobaric closure

**Compressible flow: Density correction**

To see how density compression term comes from, we
assume ∇ρ · ~n ∼ ∇α_{1} · ~n & consider case when

K_{α}_{1} = D∇α_{1} · ~n − α_{1} (1 − α_{1}) ≈ 0 =⇒ D∇α_{1} · ~n ≈ α_{1}(1 − α_{1})
yielding density diffusion normal to interface at τ-steady as

∇ (D∇ρ · ~n) · ~n ≈ ∇ (α_{1}(1 − α_{1})) · ~n = (1 − 2α_{1}) ∇α_{1} · ~n

∼ (1 − 2α_{1}) ∇ρ · ~n

Define K_{ρ} = ∇(D∇ρ · ~n) − (1 − 2α_{1}) ∇ρ & form

∂_{t}ρ + ∇ · (ρ~u) = 1

µH(α_{1})~nK_{ρ}

H(α_{1}) = tanh (α_{1}(1 − α_{1})/D)^{2} is localized-interface function

**Shukla** ^{et al.} **interface compression**

^{et al.}

In each time step, Shukla’s interface-compression algorithm for compressible 2-phase flow consists of following steps:

1** ^{.}** Solve model equation without interface-compression
terms by WENO method

2** ^{.}** Compute primitive variable w = (ρ

_{1}, ρ

_{2}, ρ, ~u, p, α

_{1}) from conservative variables q = (α

_{1}ρ

_{1}, α

_{2}ρ

_{2}, ρ~u, ρE, α

_{1})

3** ^{.}** Iterate interface & density compression equations to
τ-steady state until convergence

4** ^{.}** Update conserved variables at end of time step from

primitive variables in step 2 & new values of ρ , α_{1} from
step 3

**Underwater explosion (UNDEX)**

Solution adpated from Shukla’s paper (JCP 2010)

**Underwater explosion**

Solution adpated from Shyue’s paper (JCP 2006)

**Shukla** ^{et al.} **algorithm: Remarks**

^{et al.}

In Shukla’s results there are noises in pressure contours for UNDEX means poor calculation of pressure near interface To understand method better, consider simple interface only problem where p & ~u are constants in domain, while ρ &

material quantities in EOS have jumps across interfaces Assume consistent approximation in step 1 for model

equation without interface-compression, yielding

smeared (α_{1}ρ_{1}, α_{2}ρ_{2}, α_{1})^{∗} & retain (~u, p)^{∗} = (~u, p)
In step 3, ρ^{∗} = (α_{1}ρ_{1})^{∗} + (α_{2}ρ_{2})^{∗} & α^{∗}_{1} are compressed to ρ˜ &

˜

α_{1}, which in step 4, total mass & momentum are set

(ρ, ρu)^{n+1} = (˜ρ, ˜ρ~u^{∗}) =⇒ ~u^{n+1} = ˜ρ~u^{∗}/˜ρ = ~u^{∗} as expected

**Shukla** ^{et al.} **algorithm: Remarks**

^{et al.}

In addition, for total energy, we set

(ρE)^{n+1} =

1

2ρ|~u|^{2} + ρe

^{n+1}

= 1

2ρ|~˜ u^{∗}|^{2} + ρe(?)e

Consider stiffened gas EOS for phasic pressure
p_{k} = (γ_{k} − 1) (ρe)_{k} − γ_{k}B_{k}, k = 1, 2. We then have

e ρe =

X2 k=1

α^{k}ρ^{k}e^{k} =

X2 k=1

˜

α^{k} p^{∗} + γ^{k}B^{k}
γ^{k} − 1

= p^{∗}

X2 k=1

˜
α^{k}

γ^{k} − 1 +

X2 k=1

˜

α^{k} γ^{k}B^{k}
γ^{k} − 1

yielding equilibrium pressure p^{n+1} = p^{∗} if

1 ^{n+1}

=

X2 α˜^{k}

&

γB ^{n+1}

=

X2

˜

α^{k} γ^{k}B^{k}

**Shukla** ^{et al.} **algorithm: Remarks**

^{et al.}

Next example concerns linearized Mie-Grüneisen EOS for
phasic pressure p_{k} = (γ_{k} − 1) (ρe)_{k} + (ρ_{k} − ρ_{0k}) B_{k}

e ρe =

X2 k=1

α_{k}ρ_{k}e_{k} =

X2 k=1

˜

α_{k}p^{∗}

γ_{k} − 1 − (α˜_{k}ρ^{∗}_{k} − α˜_{k}ρ_{0k}) B_{k}
γ_{k} − 1

= p^{∗}

X2 k=1

˜
α_{k}

γ_{k} − 1 −

X2 k=1

(α˜_{k}ρ^{∗}_{k} − α˜_{k}ρ_{0k}) B_{k}
γ_{k} − 1
yielding equilibrium pressure p^{n+1} = p^{∗} if

1 γ − 1

^{n+1}

=

X2 k=1

˜
α^{k}

γ^{k} − 1 &

(ρ − ρ0)B γ − 1

^{n+1}

=

X2 k=1

(α˜^{k}ρ^{∗}k − α˜^{k}ρ0k) B^{k}
γ^{k} − 1

**Shukla** ^{et al.} **algorithm: Remarks**

^{et al.}

In Shukla * ^{et al.}* algorithm, there is a consistent problem as
X2

k=1

(α_{k}ρ_{k})^{n+1} =

X2 k=1

˜

α_{k}ρ^{∗}_{k} 6= ρ˜ = ρ^{n+1}

One way to remove this inconsistency is to include

compression terms in partial density α_{k}ρ_{k} directly, k = 1, 2,

∂_{t} (α_{k}ρ_{k}) + ∇ · (α_{k}ρ_{k}~u) =
1

µH(α_{k})~n · (∇ (D∇ (α_{k}ρ_{k}) · ~n) − (1 − 2α_{k}) ∇ (α_{k}ρ_{k}))
We then set ρ^{n+1} = P2

k=1 (α_{k}ρ_{k})^{n+1} = P2

k=1 α˜_{k}ρ˜_{k}
Validation of this approach is required

**Positivity & accuracy**

In compressible multiphase flow, positivity of volume
fraction, * ^{i.e.}*, α

_{k}≥ 0, ∀k, is important due to provision of 1

**information on interface location**

^{.}2** ^{.}** information on thermodynamic states such as ρe & p in
numerical “mixture” region & so ρ

_{k}from α

_{k}ρ

_{k}

It is known that devise of oscillation-free higher-order

method (WENO, DG, or variant) for multiphase flow is still an open problem

In this regards, interface-sharpening of some kind should be a useful tool as opposed to Eulerian higher-order

methods or other adaptive mesh methods

**Anti-diffusion interface sharpening**

Our second interface-sharpening model concerns

anti-diffusion proposed by So, Hu & Adams (JCP 2011)

∂_{t}α + ~u · ∇α = − 1

µ∇ · (D∇α), D > 0, µ ≫ 1 Standard fractional step method may still apply

1** ^{.}** Advection step over a time step

∂_{t}α + ~u · ∇α = 0
2** ^{.}** Anti-diffusion step towards sharp layer

∂_{τ}α = −∇·(D∇α) or ∂_{τ}α = −∇·(D∇α · ~n) ~n, τ = t/µ
Numerical regularization is required such as employ ^{M}^{INMOD}

**Square wave passive advection (revisit)**

Square-wave pluse moving with u = 1 after 4 periodic cycle

0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8 1

No sharpen Compress Anti−diffus Exact

α

time=4

**Vortex in cell (revisit)**

Contours α = (0.05, 0.5, 0.95) at 6 different times in 1 period No interface sharpening (second order)

With interface compression

With anti-diffusion

**Deformation flow in 3D**

In this test, consider velocity field

~u = 2 sin^{2} (πx) sin (2πy) sin (2πz), − sin (2πx) sin^{2} (πy) sin (2πz),

− sin (2πx) sin (2πy) sin^{2} (πz)

cos (πt/3)

**Deformation flow in 3D**

No anti-diffusion With anti-diffusion

**Deformation flow in 3D**

No anti-diffusion With anti-diffusion

**Anti-diffusion runs**

Methods used here are essentially the same as artificial
interface compression runs, * ^{i.e.}*,

1** ^{.}** Use Clawpack for advection in Step 1

2** ^{.}** Use first order explicit method for anti-diffusion in Step 2
Diffusion coefficient D = max |~u|

Time step ∆τ

∆τ ≤ 1 2D

Xd i=1

∆x^{2}_{i}

Stopping criterion: some measure of interface sharpness

**Anti-diffusion to compressible flow**

Reduced 2-phase model with anti-diffusion (Shyue 2011)

∂_{t}α_{1} + ~u · ∇α_{1} = − 1

µ∇ · (D∇α_{1}) = −1

µK_{α}_{1}

∂_{t} (α_{1}ρ_{1}) + ∇ · (α_{1}ρ_{1}~u) = − 1

µH_{I}∇ · (D∇α_{1}ρ_{1}) = −1

µH_{I}K_{α}_{1}_{ρ}_{1}

∂_{t} (α_{2}ρ_{2}) + ∇ · (α_{2}ρ_{2}~u) = − 1

µH_{I}∇ · (D∇α_{2}ρ_{2}) = −1

µH_{I}K_{α}_{2}_{ρ}_{2}

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = − 1

µH_{I}~u ∇ · (D∇ρ) = − 1

µH_{I}K_{ρ~}_{u}

∂_{t} (ρE) + ∇ · (ρE~u + p~u) = − 1

µH_{I} K_{ρ|~}_{u|}^{2}_{/2} + K_{ρe}

Isobaric closure for mixture pressure is used as usual

H_{I} denotes interface indicator, K denotes “diffusion” term

**Anti-diffusion to compressible flow**

To find K_{ρ|~}_{u|}^{2}_{/2} assuming |~u|^{2} is constant, we observe

∇

1

2ρ|~u|^{2}

= 1

2|~u|^{2}∇ρ yielding K_{ρ|~}_{u|}^{2}_{/2} = 1

2|~u^{2}|∇ · (D∇ρ)
To find K_{ρe}, we need to know equation of state. Now in

stiffened gas case with p_{k} = (γ_{k} − 1) (ρe)_{k} − γ_{k}B_{k},

∇(ρe) = ∇

X2 k=1

α_{k}ρ_{k}e_{k}

!

= ∇

X2 k=1

α_{k} p + γ_{k}B_{k}
γ_{k} − 1

!

=

X2 k=1

p + γ_{k}B_{k}
γ_{k} − 1

∇α_{k} =

p + γ_{1}B_{1}

γ_{1} − 1 − p + γ_{2}B_{2}
γ_{2} − 1

∇α_{1}

**Anti-diffusion to compressible flow**

We next consider case with linearized Mie-Grüneisen EOS
p_{k} = (γ_{k} − 1) (ρe)_{k} + (ρ_{k} − ρ_{0k}) B_{k} k = 1, 2, & proceed same
procedure as before

∇(ρe) = ∇

X2 k=1

α_{k}ρ_{k}e_{k}

!

= ∇

X2 k=1

α_{k} p − (ρ_{k} − ρ_{0k})B_{k}
γ_{k} − 1

!

=

X2 k=1

p + ρ_{0k}B_{k}

γ_{k} − 1 ∇α_{k} +

X2 k=1

B_{k}

γ_{k} − 1∇ (α_{k}ρ_{k})

= β_{0}∇α_{1} +

X2 k=1

β_{k}∇(α_{k}ρk)

We choose K_{ρe} = β_{0} ∇ · (D∇α_{1}) + P_{2}

1 β_{k} ∇ · (D∇α_{k}ρ_{k})

**Anti-diffusion to compressible flow**

Write anti-diffusion model in compact form

∂_{t}q + ∇ · ~f + B∇q = −1

µψ(q)
with q, f^{~}, B, & ψ defined (not shown)

In each time step, proposed anti-diffusion algorithm for compressible 2-phase flow consists of following steps:

1** ^{.}** Solve model equation without anti-diffusion terms

∂_{t}q + ∇ · ~f + B∇q = 0

2** ^{.}** Iterate model equation with anti-diffusion terms

∂_{τ}q = −ψ(q)

**Circular water column in uniform flow**

Density surface plot (moving speed ~u = (1, 1/10) )

**Single-phase Riemann problem**

0 0.2 0.4 0.6 0.8 1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

No sharpen Anti−diffu Exact

ρ

x

Sod problem

**Single-phase Riemann problem**

0 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

No sharpen Anti−diffu Exact

ρ

x

Lax problem

**Oscillating water column**

Initially, in closed shock tube, water column moves at u = 1 from left to right, yielding air compression at right

& air expansion at left

Subsequently, pressure difference built up across water column resulting deceleration of column of water to

right, makes a stop, & then acceleration to left; a reverse pressure difference built up across water column redirecting flow from left to right again

Eventually, water column starts to oscillate

air

air water

u →

**Oscillating water column**

0.6 0.8 1 1.2 1.4 1.6 1.8 2

pL

With anti−diffusion No anti−diffusion

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2

pR

### 2 **D Riemann Problem**

With initial 4-slip lines wave pattern

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

ρ u v p

^{=}

1 0.75

−0.5 1

2 0.75

0.5 1

1

−0.75 0.5

1

3

−0.75

−0.5 1

### 2 **D Riemann Problem**

Density on top & pressure on bottom

**Shock in air & water cylinder**

Solution adpated from Shukla’s paper (JCP 2010)

**Shock in air & water cylinder**

Schlieren images for anti-diffusion results at times t = 0.15, 0.4, 0.65 (volume fraction on right)

**Anti-diffusion sharpening: Remarks**

1. Local interface identification Algebraic

H(α) = tanh (α(1 − α)/D)^{2} (Shukla * ^{et al.}* )
Physical-jump (density, volume fraction, . . . )

2. Diffusion coefficient definition

Global D = max |~u| or D_{j} = max |u_{j}|

Local D = max_{M (C)} |~u| or D_{j}^{C} = max_{M (C)} |u_{j}|
3. Stopping criterion

Run 1 − 2 anti-diffusion iteration currently Interface-shapreness measure (?)

**Effect of diffusion coefficient**

Anti-diffusion results for Zalesak’s rotating disc

Global D_{j} Local D_{j}

**Effect of stopping criterion**

Anti-diffusion results for Zalesak’s rotating disc

Local D_{j}: 1 step Local D_{j}: 2 step

**Future perspectives**

Extend method to mapped grid & model with phase

transition proposed by Saurel, Petipas, Abgrall (JFM 2008)

∂_{t} (α_{1}ρ_{1}) + ∇ · (α_{1}ρ_{1}~u) = 0

∂_{t} (α_{2}ρ_{2}) + ∇ · (α_{2}ρ_{2}~u) = 0

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_{t} (ρE) + ∇ · (ρE~u + p~u) = 0

∂_{t}α_{1} + ~u · ∇ (α_{1}~u) = α_{1} K¯_{S}

K_{S}^{1} ∇ · ~u + K˜_{Γ}

K˜_{S} Q_{1} + K˜_{c}
K˜_{S} ρ ˙Y
K¯_{S} =

α_{1}

K_{S}^{1} + α_{2}
K_{S}^{2}

−1

, K˜_{S} =

K_{S}^{1}

α_{1} + K_{S}^{2}
α_{2}

, K˜_{Γ} =

Γ_{1}

α_{1} + Γ_{2}
α_{2}

K˜_{c} =

c^{2}_{1}

+ c^{2}_{2}

, K^{ι} = ρ_{ι}c^{2}_{ι} , Q_{1} = H(T_{2} − T_{1}), ˙Y (mass trans.)