II: Eulerian interface sharpening approach

Volume of fluid methods for

compressible multiphase flow

II: Eulerian interface sharpening approach

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

Objective

Recall that our aim is to discuss a class of volume of fluid (vs. level set, MAC, particles) methods for interface

problems with application to compressible multiphase flow 1^. Adaptive moving grid approach (last lecture)

Cartesian grid embedded volume tracking Moving mapped grid interface capturing

2^. Eulerian interface sharpening approach (this lecture) Artificial interface compression method

Anti-diffusion method

Outline

Review interface sharpening techniques for viscous incompressibe two-phase flow

Artificial interface compression Anti-diffusion

Extend method to compressible multiphase flow Interface only problem

Problem with shock wave

This is a work in progress since August 2011

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, where

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

, f~^σ = −σκ∇α with κ = ∇ ·

 ∇α 

Interface sharpening techniques

Typical interface sharpening methods for incompressible 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 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 artificial compression proposed by Olsson & Kreiss

∂^tα + ~u · ∇α = 1

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

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

Standard fractional step method may apply as 1^. Advection step over a time step

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

2^. Interface compression step to

Square wave passive advection

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 compress With compress Exact

α

Interface compression: 1D case

To see why this approach works, consider 1D model

∂^tα + u∂^xα = 1

µ∂^x~n · [D (∂^xα · ~n) − α (1 − α)], x ∈ lR, t > 0 with ~n = 1 & initial α(x, 0) = α0(x) = 1/ (1 + exp (−x/D)).

Exact solution for this initial value problem is simply α(x, t) = α0(x − ut)

while solution with perturbed data α˜0(x) = α0(x) + δ(x) is α(ξ, τ ) = ˜α0 (ξ + ξ0) as τ → ∞ (ξ = x − ut, τ = t/µ) If perturbation is zero mass R ^∞

−∞ δ(ξ, 0)dξ = 0 (which is true if model is solved conservatively), we have true solution with ξ0 = 0, see Sattinger (1976) & Goodman (1986)

Interface compression: Multi-D case

Let K = D∇α · ~n − α (1 − α) with ~n = ∇α/|∇α|. In interface-compression step, we solve

∂^τα = ∇ · ~n [D (∇α · ~n) − α (1 − α)] = K∇ · ~n + ~n · ∇K

& reach τ-steady state solution as µ → ∞, yielding K = 0 &

1D profile in coordinate normal to interface n^⊥ as α = 1/ 1 + exp (−n^⊥/D)

= 1 + tanh (n^⊥/2D)
/2

When µ finite, K 6= 0, ^i.e., K∇ · ~n + ~n · ∇K 6= 0, α & so interface would be changed both normally & tangentially depending on both strength & accuracy of curvature ∇ · ~n evaluation numerically

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 runs

Methods used here are very elementary, ^i.e., 1^. Use Clawpack for advection in Step 1

2^. Use simple first order explicit method for interface compression in Step 2

Diffusion coefficient D = ε min^∀i ∆xⁱ 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 (α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(α¹)~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 ∼ ∇α1 · ~n & consider case when

K = D∇α1·~n−α1 (1 − α1) ≈ 0 =⇒ D∇α1 · ~n ≈ α1(1 − α1) yielding density diffusion normal to interface as

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

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

Analogously, we write density equation with correction as

∂^tρ + ∇ · (ρ~u) = 1

µH(α1)~n · (∇ (D∇ρ · ~n) − (1 − 2α1) ∇ρ) H(α1) = tanh (α1(1 − α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 state-of-the-art shock capturing method

2^. Compute promitive variable w = (ρ1, ρ2, ρ, ~u, p, α1) 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 ρ , α1 from step 3

Underwater explosion

Solution adpated from Shukla’s paper (JCP 2010)

Shock in air & water cylinder

Solution adpated from Shukla’s paper (JCP 2010)

Underwater explosion (revisit)

Solution adpated from Shyue’s paper (JCP 2006)

Shukla ^{et al.} algorithm: Remark

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

˜

α1, which in step 4, for total mass & momentum, we set

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

Shukla ^{et al.} algorithm: Remark

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: Remark

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: Remark

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 α^kρ^k, k = 1, 2, via

∂^t (α^kρ^k) + ∇ · (α^kρ^k~u) = 1

µH(α1)~n · (∇ (D∇ (α^kρ^k) · ~n) − (1 − 2α1) ∇ (α^kρ^k)) We set ρⁿ⁺¹ = P²

k=1 (α^kρ^k)ⁿ⁺¹ = P²

k=1 α˜^kρ˜^k

Validation of this approach is required

Anti-diffusion interface sharpening

Alternative interface-sharpening model is anti-diffusion proposed by So, Hu & Adams (JCP 2011)

∂^tα + ~u · ∇α = − 1

µ∇ · (D∇α), D > 0, µ ≫ 1 Standard fractional step method may apply again as 1^. Advection step over a time step

∂^tα + ~u · ∇α = 0 2^. Anti-diffusion step towards “sharp layer”

∂^τα = −∇·(D∇α) or ∂^τα = −∇·~n (D∇α · ~n) , τ = t/µ Numerical regularization is required such as employ ^MINMOD limiter to stabilize ∇α in discretization, Breuß ^{et al.} (’05, ’07)

Square wave passive advection (revisit)

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

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

Passive advection: Conservation

Mass error for sample interface sharpening methods

0 1 2 3 4

−1

−0.5 0 0.5 1

compress anti−diffu

Masserror

time

Anti-diffusion runs

Methods used here are 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

Positivity & accuracy

In compressible multiphase flow, positivity of volume

fraction, ^i.e., α^k ≥ 0, ∀k, is of fundamental importance ; this is because it provides, in particular,

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 & have been mentioned many times in this conference that devise of oscillation-free higher-order

method is still an open problem (if I have stated correctly) In this regards, interface-sharpening of some kind should be a useful tool as opposed to higher-order methods or other volume-of-fluid methods

Anti-diffusion to compressible flow

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

∂^tα1 + ~u · ∇α1 = −1

µ∇ · (D∇α1)

∂^t (α1ρ1) + ∇ · (α1ρ1~u) = − 1

µH(α1)∇ · (D∇α1ρ1)

∂^t (α2ρ2) + ∇ · (α2ρ2~u) = − 1

µH(α1)∇ · (D∇α2ρ2)

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

µH(α1)~u ∇ · (D∇ρ)

∂^t (ρE) + ∇ · (ρE~u + p~u) =

− 1

µH(α1)L

1

2ρ|~u|²



− 1

µH(α1)L (ρe)

Anti-diffusion to compressible flow

To find L ¹₂ρ|~u|²

, assuming |~u|² is constant, we observe

∇

1

2ρ|~u|²



= 1

2|~u|²∇ρ yielding L

1

2ρ|~u|²



= 1

2|~u²|∇ · (D∇ρ) To find L (ρ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 + γ1B1

γ1 − 1 − p + γ2B2

γ2 − 1



∇α1

= β ∇α yielding L (ρe) = β ∇ · (D∇α )

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)

= β0∇α1 +

X2

k=1

β^k∇(α^kρk)

We choose L (ρ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) to τ-steady state until convergence

Circular water column advection

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

Air-Helium Riemann problem

0 0.2 0.4 0.6 0.8 1

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

rho

x

No anti−diffu Anti−diffu Exact

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10 12 14

u

x

0 100 200 300 400 500

u

0 0.2 0.4 0.6 0.8 1

alf

Future perspectives

Effect of local interface identification H(α) Algebraic approach

H(α) =

( tanh (α(1 − α)/D)² (Shukla ^{et al.} ) tanh p

α(1 − α)/D PDE approach

Anti-diffusion on moving mapped grid Extension to other multiphase model

Thank you