interface-sharpening methods for

Eulerian

interface-sharpening methods for

hyperbolic problems

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 − α) ǫ₂, α ∈ [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² (πx) sin (2πy), sin (2πx) sin² (π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) = α₀(x) = 1/ (1 + exp (−x/D)) , x ∈ lR, t > 0 Exact solution for this problem is simply α(x, t) = α₀(x − ut) When α(x, 0) is perturbed to α₀(x) + δ(x), δ(x) ≪ 1, we have

∂_τα + ∂˜ _ξ α˜²/2

= ∂_ξ D∂_ξα˜

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

Interface compression: Remarks

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

−∞ δ(ξ, 0)dξ = 0 we have true solution with ξ₀ = 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²_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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

∂_t (α₂ρ₂) + ∇ · (α₂ρ₂~u) = 0

∂_t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_t(ρE) + ∇ · (ρE~u + p~u) = 0

∂_tα₁ + ~u · ∇α₁ = 1

µ~n · ∇ (D∇α₁ · ~n − α₁ (1 − α₁))

∂_tρ + ∇ · (ρ~u) = 1

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

Compressible flow: Density correction

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

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

∇ (D∇ρ · ~n) · ~n ≈ ∇ (α₁(1 − α₁)) · ~n = (1 − 2α₁) ∇α₁ · ~n

∼ (1 − 2α₁) ∇ρ · ~n

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

∂_tρ + ∇ · (ρ~u) = 1

µH(α₁)~nK_ρ

H(α₁) = tanh (α₁(1 − α₁)/D)² is localized-interface function

Shukla ^{et al.} interface compression

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 = (ρ₁, ρ₂, ρ, ~u, p, α₁) from conservative variables q = (α₁ρ₁, α₂ρ₂, ρ~u, ρE, α₁)

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 ρ , α₁ 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

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 (α₁ρ₁, α₂ρ₂, α₁)^∗ & retain (~u, p)^∗ = (~u, p) In step 3, ρ^∗ = (α₁ρ₁)^∗ + (α₂ρ₂)^∗ & α^∗₁ are compressed to ρ˜ &

˜

α₁, which in step 4, total mass & momentum are set

(ρ, ρu)ⁿ⁺¹ = (˜ρ, ˜ρ~u^∗) =⇒ ~uⁿ⁺¹ = ˜ρ~u^∗/˜ρ = ~u^∗ as expected

Shukla ^{et al.} algorithm: Remarks

In addition, for total energy, we set

(ρE)ⁿ⁺¹ =

1

2ρ|~u|² + ρe

ⁿ⁺¹

= 1

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

Consider stiffened gas EOS for phasic pressure p_k = (γ_k − 1) (ρe)_k − γ_kB_k, k = 1, 2. We then have

e ρe =

X2 k=1

α^kρ^ke^k =

X2 k=1

˜

α^k p^∗ + γ^kB^k γ^k − 1

= p^∗

X2 k=1

˜ α^k

γ^k − 1 +

X2 k=1

˜

α^k γ^kB^k γ^k − 1

yielding equilibrium pressure pⁿ⁺¹ = p^∗ if

 1 ⁿ⁺¹

=

X2 α˜^k

&

 γB ⁿ⁺¹

=

X2

˜

α^k γ^kB^k

Shukla ^{et al.} algorithm: Remarks

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ρ_ke_k =

X2 k=1

˜

α_kp^∗

γ_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ⁿ⁺¹ = p^∗ if

 1 γ − 1

ⁿ⁺¹

=

X2 k=1

˜ α^k

γ^k − 1 &

(ρ − ρ0)B γ − 1

ⁿ⁺¹

=

X2 k=1

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

Shukla ^{et al.} algorithm: Remarks

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

k=1

(α_kρ_k)ⁿ⁺¹ =

X2 k=1

˜

α_kρ^∗_k 6= ρ˜ = ρⁿ⁺¹

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 ρⁿ⁺¹ = P2

k=1 (α_kρ_k)ⁿ⁺¹ = 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 ^MINMOD

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² (πx) sin (2πy) sin (2πz), − sin (2πx) sin² (πy) sin (2πz),

− sin (2πx) sin (2πy) sin² (π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²_i

Stopping criterion: some measure of interface sharpness

Anti-diffusion to compressible flow

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

∂_tα₁ + ~u · ∇α₁ = − 1

µ∇ · (D∇α₁) = −1

µK_α1

∂_t (α₁ρ₁) + ∇ · (α₁ρ₁~u) = − 1

µH_I∇ · (D∇α₁ρ₁) = −1

µH_IK_α1ρ1

∂_t (α₂ρ₂) + ∇ · (α₂ρ₂~u) = − 1

µH_I∇ · (D∇α₂ρ₂) = −1

µH_IK_α2ρ2

∂_t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = − 1

µH_I~u ∇ · (D∇ρ) = − 1

µH_IK_ρ~u

∂_t (ρE) + ∇ · (ρE~u + p~u) = − 1

µH_I K_ρ|~u|²_/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 assuming |~u|² is constant, we observe

∇

1

2ρ|~u|²



= 1

2|~u|²∇ρ yielding K_ρ|~u|²_/2 = 1

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

stiffened gas case with p_k = (γ_k − 1) (ρe)_k − γ_kB_k,

∇(ρe) = ∇

X2 k=1

α_kρ_ke_k

!

= ∇

X2 k=1

α_k p + γ_kB_k γ_k − 1

!

=

X2 k=1

p + γ_kB_k γ_k − 1



∇α_k =

p + γ₁B₁

γ₁ − 1 − p + γ₂B₂ γ₂ − 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ρ_ke_k

!

= ∇

X2 k=1

α_k p − (ρ_k − ρ_0k)B_k γ_k − 1

!

=

X2 k=1

p + ρ_0kB_k

γ_k − 1 ∇α_k +

X2 k=1

B_k

γ_k − 1∇ (α_kρ_k)

= β₀∇α₁ +

X2 k=1

β_k∇(α_kρk)

We choose K_ρe = β₀ ∇ · (D∇α₁) + P₂

1 β_k ∇ · (D∇α_kρ_k)

Anti-diffusion to compressible flow

Write anti-diffusion model in compact form

∂_tq + ∇ · ~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

∂_tq + ∇ · ~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)² (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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

∂_t (α₂ρ₂) + ∇ · (α₂ρ₂~u) = 0

∂_t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_t (ρE) + ∇ · (ρE~u + p~u) = 0

∂_tα₁ + ~u · ∇ (α₁~u) = α₁ K¯_S

K_S¹ ∇ · ~u + K˜_Γ

K˜_S Q₁ + K˜_c K˜_S ρ ˙Y K¯_S =

 α₁

K_S¹ + α₂ K_S²

−1

, K˜_S =

K_S¹

α₁ + K_S² α₂



, K˜_Γ =

Γ₁

α₁ + Γ₂ α₂



K˜_c =

 c²₁

+ c²₂ 

, K^ι = ρ_ιc²_ι , Q₁ = H(T₂ − T₁), ˙Y (mass trans.)

Thank you