Shock in air & R22 bubble interaction

Eulerian

interface-sharpening methods for

compressible flow

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Objective

Derive consistent interface-sharpening model & method for compressible (single & multiphase) flow problems

PDE-based approach

∂_tq + ∇ · f (q) + B∇q = 1 µD_q

Algebraic-based approach (with F. Xiao, Tokyo Tech.) 1. Sharp cell edge solution reconstruction

2. Solution update: MUSCL or semi-discretize scheme

Shock in air & R22 bubble interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock in air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

Shock air-R22 bubble (cont.)

With anti-diffusion With THINC

Shock air-R22 bubble (cont.)

Volume tracking (Shyue 2006) 2nd order

PDE-based interface sharpening

Incompressible 2-phase flow: PDE-based interface sharpening for volume-fraction transport

∂_tα + u · ∇α = 1

µD_α, µ ∈ lR ≫ 1

Artificial compression: Harten CPAM 1977, Olsson &

Kreiss JCP 2005

D_α := ∇ · [(D(∆x)∇α · ~n − α (1 − α)) ~n]

Anti-diffusion: So, Hu & Adams JCP 2011 D_α := −∇ · (D(u)∇α)

Model interface-only problem

Interface-only problem for compressible 1-phase Euler equations with constant pressure p, velocity u, & jump in density ρ across interfaces

∂_tρ + ∇ · (ρu) = 0 (Mass)

∂_t (ρu) + ∇ · (ρu ⊗ u + pI_N) = 0 (Momentum)

∂_tE + ∇ · (Eu + pu) = 0 (Energy) yielding basic transport equations for interfaces as

∂_tρ + u · ∇ρ = 0 u (∂_tρ + u · ∇ρ) = 0 u · u

2 (∂_tρ + u · ∇ρ) + ∂_t (ρe) + u · ∇ (ρe) = 0

Interface-only anti-diffusion model

Shyue (2011) proposed anti-diffusion model for density

∂_tρ + u · ∇ρ = 1

µD_ρ, D_ρ := −∇ · (D∇ρ)

To ensure velocity equilibrium for momentum

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

µD_ρu, D_ρu := u D_ρ

& velocity-pressure equilibrium, for total energy

u · u

2 (∂_tρ + u · ∇ρ) + ∂_t (ρe) + u · ∇ (ρe) = 1

µD_E, D_E := hu · u

2 + ∂_ρ(ρe)i D_ρ

Mie-Grüneisen EOS p(ρ, e) = p^∞(ρ) + Γ(ρ)ρ [e − e^∞(ρ)]

1 -phase anti-diffusion model

To deal with shock waves, anti-diffusion model for compressible 1-phase flow

∂_tρ + ∇ · (ρu) = 1 µD_ρ

∂_t (ρu) + ∇ · (ρu ⊗ u + pI_N) = 1

µD_ρu

∂_tE + ∇ · (E~u + pu) = 1

µD_E D_ρ := −H_I∇ · (D∇ρ) , D_ρu := u D_ρ, D_E := h_{u · u}

2 + ∂_ρ(ρe)i D_ρ H_I : Interface indicator =

( 1 if near interface, 0 otherwise

Anti-diffusion method

Anti-diffusion model in compact form

∂_tq + ∇ · f = 1 µD_q with q, f, & ^D_q defined (not shown)

Fractional step method:

1^. Solve homogeneous equation without source terms

∂_tq + ∇ · f = 0

2^. Iterate model equation with source terms

∂_τq = D_q

to sharp layer; τ = µt (pseudo time)

Numerical interface-only problem

Consistency of numerical solution for interface-only problem Step 1: assume consistent approximation model equation without anti-diffusion,

smeared ρ^∗ & equilibrium u^∗ = uⁿ, p^∗ = pⁿ

(may be difficult with highly nonlinear MG EOS)

Step 2: assume anti-diffusion for ρ is done consistently &

stably, yielding update ρ^∗ to ρⁿ⁺¹ & for momentum,

(ρu)ⁿ⁺¹ := (ρu)^∗ + D_ρu = (ρu)^∗ + u^∗D_ρ

= (ρu)^∗ + u^∗ ρⁿ⁺¹ − ρ^∗
= ρⁿ⁺¹u^∗ =⇒ uⁿ⁺¹ = ρⁿ⁺¹u^∗

ρⁿ⁺¹ = u^∗

Numerical interface-only (cont.)

For total energy, Eⁿ⁺¹ := E^∗ + D_E, ^i.e.,

(ρK + ρe)ⁿ⁺¹ := (ρK + ρe)^∗ + (K + ∂_ρ(ρe))^∗ D_ρ, K = u · u 2

Splitting kinetic & internal energy with MG EOS

(ρK)ⁿ⁺¹ := (ρK)^∗ + K^∗ ρⁿ⁺¹ − ρ^∗
= ρⁿ⁺¹K^∗ =⇒ Kⁿ⁺¹ = K^∗

=⇒ uⁿ⁺¹ = u^∗ (not true if u has transverse jump)

 p − p^∞

Γ + ρe^∞

n+1

:=  p − p^∞

Γ + ρe^∞

^∗

+ ∂_ρ  p − p^∞

Γ + ρe^∞

^∗

ρⁿ⁺¹ − ρ^∗

=⇒ pⁿ⁺¹ = p^∗ (for linearized MG EOS only)

Numerical shock wave problems

Weak solution for problems with shock & rarefaction waves Interface indicator H_I takes value zero away from interfacs, yielding standard compressible Euler equations in

conservation form

Step 1, use state-of-the-art shock capturing method for entropy-satisfying weak solutions

Implementation issues

Methods used here are very elementary, ^i.e.,

Step 1: Clawpack (or Wenoclaw/Sharpclaw) for

homogeneous equation over CFL-constrained ∆t Step 2: explicit 1-step method for anti-diffusion

Numerical regularization to ^∇ρ & ^D_ρ

Diffusion coefficient D = D(u) (local in space & time) Time step (forward Euler) ∆τ ≤ min 

∆t, ^min

N

i=1 ∆x²ⁱ 2N D^max



Stopping criterion: Run 1 − 2 iteration currently H_I is chosen based on checking jumps in ρ & p

Sod Riemann problem

High-resolution result with anti-diffusion

Ideal gas: p(ρ, e) = (γ − 1)ρe

0 0.5 1

0 0.2 0.4 0.6 0.8 1

0 0.5 1

0 0.2 0.4 0.6 0.8 1

0 0.5 1

0 0.2 0.4 0.6 0.8

Density Velocity 1 Pressure

x x

x

Aluminum impact problem

rest Al ^← moving Al IC:

(ρ, u, p)_L = 

4 × 10³ kg/m³, 0 m/s, 7.93 × 10⁹ Pa (ρ, u, p)_R = 

2785 kg/m³, −2 × 10³ m/s, 0 Pa

Mie-Grüneisen EOS with

Γ(ρ) = Γ0(1 − η)^α, p^∞(ρ) = ρ0c²₀η

(1 − sη)² , e^∞(ρ) = η 2ρ0

(p0 + p^∞(ρ)) ρ0(kg/m³) c0 (m/s) s Γ0 α p0 e0

2785 5328 1.338 2.0 1 0 0

Aluminum impact (cont.)

High-resolution result with anti-diffusion at time t = 50µs

0 0.5 1

2.5 3 3.5 4 4.5

sharpen no sharpen exact

0 0.5 1

−2

−1.5

−1

−0.5 0

0 0.5 1

−5 0 5 10 15 20 25

Density (Mg/m³) Velocity (km/s) 30Pressure (GPa)

x x

x

1 -fluid multiphase anti-diffusion

Compressible multiphase: One-fluid anti-diffusion model

∂_tρ + ∇ · (ρu) = 1 µD_ρ

∂_t (ρu) + ∇ · (ρu ⊗ u + pI_N) = 1

µD_ρu

∂_tE + ∇ · (E~u + pu) = 1 µD_E

∂_tφ_j + u · ∇φ_j = ρd_ρφ_j∇ · u + 1

µD_φj

Mixture pressure modeled by Mie-Grüneisen EOS:

φ1 = 1

Γ(ρ), φ2 = p^∞(ρ)

Γ(ρ) , φ3 = ρe^∞(ρ)

µ → ∞ reduces to fluid-mixture model (Shyue JCP 2001)

5 -equation 2-phase flow model

Unsteady, inviscid, compressible homogeneous 2-phase flow governed by 5-equation model

∂_t (α1ρ1) + ∇ · (α1ρ1u) = 0 (Phasic 1 continuity)

∂_t (α2ρ2) + ∇ · (α2ρ2u) = 0 (Phasic 2 continuity)

∂_t (ρu) + ∇ · (ρu ⊗ u + pI_N) = 0 (Momentum)

∂_tE + ∇ · (Eu + pu) = 0 (Energy)

∂_tα1 + u · ∇α1 = 0 (Volume fraction) Phasic pressure p_k follows Mie-Grüneisen EOS

p_k(ρ_k, e_k) = p^∞_,k(ρ_k) + Γ_k(ρ_k)ρ_k 

e_k − e^∞_,k(ρ_k)

k = 1, 2

5 -equation model (cont.)

Isobaric closure leads to mixture pressure , ^i.e., Substitute p = p1 = p2 in ρe =

2

X

k=1

ρ_ke_k, yielding

p = ρe −

2

X

k=1

α_kρ_ke^∞_,k(ρ_k) +

2

X

k=1

α_k p^∞_,k(ρ_k) Γ_k(ρ_k)

! ² X

k=1

α_k Γ_k(ρ_k)

Model is hyperbolic with mixture acoustic impedance Allaire ^{et al.} (JCP 2002)

ρc² =

2

X

k=1

α_kρ_kc²_k, c_k : phasic sound speed

5 -equation model (cont.)

In cavitated regions, p < psat, cutoff model Non-conservative energy correction

E := Esat =

2

X

k=1

α_k(ρ_ke_k)sat + ρK

cutoff phasic internal energy is

(ρ_ke_k)sat = psat − p^∞_,k(ρsat)

Γ_k(ρsat) + ρsate^∞_,k(ρsat)

α-based energy-preserving correction (Shyue 2012)

α1 := (α1)sat = ρe − (ρ2e2)sat (ρ1e1)sat − (ρ2e2)sat

5 -equation anti-diffusion model

Proposed 5-equation anti-diffusion model

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

µD_α

1ρ₁

∂_t (α2ρ2) + ∇ · (α2ρ2u) = 1

µD_α2ρ2

∂_t (ρu) + ∇ · (ρu ⊗ u + pI_N) = 1

µD_ρu

∂_tE + ∇ · (Eu + pu) = 1

µD_E

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

µD_α1

Exist two ways to set ^D_z, z = α1ρ1, . . . , α1, in literature

5 -equation anti-diffusion (cont.)

α-ρ based (Shyue 2011)

D_α1 := −∇ · (D∇α1), D_αkρk := −H_I∇ · (D∇α_kρ_k), k = 1, 2 D_ρ :=

2

X

k=1

D_αkρk, D_ρu := uD_ρ,

D_E := KD_ρ +

2

X

k=1

∂_αkρk(ρ_ke_k)D_αkρk +

2

X

k=1

ρ_ke_kD_αk

α-based only (So, Hu, & Adams JCP 2012)

D_α1 := −∇ · (D∇α1), D_αkρk := ρ_kD_αk, k = 1, 2, D_α2 := −D_α1 D_ρ :=

2

X

k=1

D_αkρk, D_ρu := uD_ρ, D_E := KD_ρ +

2

X

k=1

ρ_ke_kD_αk

Use ^D := ∇ · [(D∇α ·~n − α (1 − α )) ~n] (Shyue 2012)

5 -equation interface-compression

Shukla, Pantano & Freund (JCP 2010)

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

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

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

∂_tE + ∇ · (Eu + pu) = 0

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

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

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

µH_I(α1)n · (∇ (D∇ρ · n) − (1 − 2α1) ∇ρ)

Nonliner compression & linear diffusion for interface sharpening & stability

Method proposed there is unstable numerically

PDE-based interface sharpening

PDE-based interface-sharpening model in compact form

∂_tq + ∇ · f + B∇q = 1 µD_q with q, f, B, & ^D_q defined (not shown)

Fractional step method is used

1^. Solve homogenous equation without source terms

∂_tq + ∇ · f + B∇q = 0 2^. Iterate model equation with source terms

∂_τq = D_q

to sharp layer; τ = µt (pseudo time)

Numerical interface-only problem

Consistency of numerical solution for interface-only problem Step 1: assume consistent approximation model equation without anti-diffusion,

smeared (α1ρ1)^∗, (α2ρ2)^∗, α₁^∗ & equilibrium u^∗ = uⁿ, p^∗ = pⁿ (can be done even with highly nonlinear MG EOS)

Step 2: update (α1ρ1)^∗, (α2ρ2)^∗, α^∗₁ to (α1ρ1)ⁿ⁺¹, (α2ρ2)ⁿ⁺¹, α₁^{n+1 via} anti-diffusion, & for momentum,

(ρu)ⁿ⁺¹ := (ρu)^∗ + D_ρu = (ρu)^∗ + u^∗D_ρ =⇒ uⁿ⁺¹ = u^∗

Numerical interface-only (cont.)

For total energy, Eⁿ⁺¹ := E^∗ + D_E, ^i.e., α-based

(ρK + ρe)ⁿ⁺¹ := (ρK + ρe)^∗ + K^∗D_ρ +

2

X

k=1

(ρ_ke_k)^∗ D_αkρk

Splitting kinetic & internal energy with MG EOS

(ρK)ⁿ⁺¹ := (ρK)^∗ + K^∗ ρⁿ⁺¹ − ρ^∗
= ρⁿ⁺¹K^∗ =⇒ Kⁿ⁺¹ = K^∗

2

X

k=1

(α_kρ_ke_k)ⁿ⁺¹ :=

2

X

k=1

(α_kρ_ke_k)^∗ +

2

X

k=1

(ρ_ke_k)^∗ αⁿ⁺¹_k − α^∗_k

2

X

k=1



α_k  p − p^∞_,k

Γ_k + ρ_ke^∞_,k

n+1

:=

2

X

k=1

 p − p^∞_,k

Γ_k + ρ_ke^∞_,k

^∗

αⁿ⁺¹_k

=⇒ pⁿ⁺¹ = p^∗ (for general MG EOS)

Piston-induced liquid drops depression

Liquid & vapor governed by stiffened gas EOS Piston velocity u = −100m/s

vapor

liquid liquid

← piston

Liquid drops depression

Vapor mass fraction Mixture pressure

t=2ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=4ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=6ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=8ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=10ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=12ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=14ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=16ms

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=18ms

Underwater explosions: cylindrical wall

High pressure gaseous explosive in water (cylindrical case)

water

explosive

UNDEX: cylindrical wall (cont.)

t = 30µs

shock

reflected rarefaction reflected shock

t = 45µs

cavitation

interface

UNDEX: cylindrical wall (cont.)

t = 60µs

cavitation

t = 120µs

cavitation

UNDEX: cylindrical wall (cont.)

0 20 40 60 80 100 120

0 2000 4000 6000 8000

With anti−difusion Without−anti−diffusion

0 20 40 60 80 100 120

−10 0 10 20 30 40 50 60

With anti−difusion Without−anti−diffusion

Pressure(atm) Velocity(m/s)

Time (µs)

cavi cutoff cavi collapse & reload

Algebraic-based interface sharpening

Incompressible flow: Algebraic-based interface sharpening 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

Compressible flow: THINC-type interface sharpening 1. Cell edge solution reconstruction: THINC

2. Solution update: MUSCL or semi-discretize scheme

THINC interface sharpening

THINC reconstruction assumes

α_i(x) = αmin + αmax − αmin

2



1 + γ tanh



βx − x_i−1/2

∆x − x¯_i



αmin = min (α_i−1, α_i+1) , αmax = max (α_i−1, α_i+1) , γ = sign(α_i+1 − α_i−1)

β measures sharpeness (given constant) & x¯_i chosen by

α_i = 1

∆x

Z x_i+1/2 x_i−1/2

α_i(x) dx

Cell edges determined by

α_i+1/2,L = αmin + αmax − αmin

2



1 + γ tanh β + C 1 + C tanh β



α_i−1/2,R = αmin + αmax − αmin

(1 + γC) (C not shown)

Homogeneous-equilibrium THINC

With α_i+1/2,L & α_i−1/2,R obtained by THINC reconstruction Homogeneous-equilibrium reconstruction (analgously to α-based 5-equation anti-diffusion model)

(α1ρ1)_i+1/2,L = ρ1,i α_1,i+1/2,L − α1,i
(α2ρ2)_i+1/2,L = ρ2,i α_2,i+1/2,L − α2,i

ρ_i+1/2,L = (α1ρ1)_i+1/2,L + (α2ρ2)_i+1/2,L (ρu)_i+1/2,L = u_iρ_i+1/2,L

E_i+1/2,L = K_iρ_i+1/2,L +

2

X

k=1

(ρ_ke_k) _i α_k,i+1/2,L − α_k,i

z_i−1/2,R can be made in a similar manner

6 -equation anti-diffusion model

Anti-diffusion to 6-equation model (Pelanti & Shyue 2012)

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

µD_α1ρ1

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

µD_α2ρ2

∂_t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1 + α2p2) = 1

µD_ρu

∂_t (α1E1) + ∇ · (α1E1~u + α1p1~u) + w (q, ∇q) = −νp_I (p1 − p2) + 1

µD_α1E1

∂_t (α2E2) + ∇ · (α2E2~u + α2p2~u) − w (q, ∇q) = νp_I (p1 − p2) + 1

µD_α2E2

∂_tα1 + ~u · ∇α1 = ν (p1 − p2) + 1

µD_α1

w = −~u ((Y2p1 + Y1p2) ∇α1 + α1Y2∇p1 − α2Y1∇p2)

Future perspective

Shape-preserving interface shapening on mapped grid Interface shapening to flow with physical sources

(usefulness ?)

∂_tq + ∇ · f (q) + B∇q = ψ(q) + 1 µD_q

Extension to low Mach weakly compressible flow Hybrid interface-sharpening & WENO towards high order for compressible turbulent flow

Thank you