Interface sharpening methods for compressible multiphase flow problems: Overview & look ahead

Interface sharpening methods for

compressible multiphase flow problems:

Overview & look ahead

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

Joint work: M. Pelanti (Paris Tech) & F. Xiao (Tokyo Tech) 11:00-12:00, September 17, 2015, IUSTI, AMU

Outline

Interface sharpening methods for model systems:

1. Passive tracer transport (motivation) 2. Compressible single-phase flow (inviscid)

3. 5-equation 1-velocity, 1-pressure model for compressible two-phase flow

4. 6-equation 1-velocity, 2-pressure model for compressible two-phase flow

5. Model for compressible two-phase flow with drift-flux approximation

How interface sharpening ?

Higher order method for sharpening interface ?

How interface sharpening ?

Higher order method for sharpening interface ?

Benchmark test for shock in air & R22 bubble interaction

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

With anti-diffusion WENO 5

WENO gives more chaotic interface, positivity-preserving in WENO5 with MG EOS (working open issue)

With anti-diffusion WENO 5

THINC gives more regularized interface

With anti-diffusion With THINC

Standard high-resolution method gives poor interface resolution

Volume tracking (Shyue 2006) 2nd order

Shock-contact interaction: Wave speeds diagnosis

Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf

Experiment 415 240 540 73 90 78 78

Quirk & Karni 420 254 560 74 90 116 82 Kokh & Lagoutiere 411 243 525 65 86 86 64

Ullah et al. 410 246 535 65 86 76 60

Shyue (tracking) 411 243 538 64 87 82 60 Capturing results 410 244 536 65 86 98 76 THINC results 410 244 538 65 86 87 64 Anti-df results 410 244 532 64 85 100 78

Toy problem: Passive tracer transport

Free-surface (or 2-phase) flow modelled by incompressible Navier-Stokes equations read

∂t(ρ~u) + div (ρ~u ⊗ ~u) + ∇p = ∇ · τ + ρ~g + ~fσ

∂tα + ~u · ∇α= 0 (volume fraction transport) div(~u) = 0

Material quantities in 2-phase region determined by z =αz1 + (1 −α) z2, z = ρ, ǫ, & σ Source terms are dependent on α also

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

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

 ∇α

|∇α|



Sharply resolved positivity-preserving α is fundamental

Standard interface capturing results for toy problem, observing poor interface resolution

Passive advection Rotating disc

Vortex in cell

How interface sharpening ? Toy problem

Eulerian interface sharpening methods (i.e., use uniform underlying grid) for volume fraction transport include

1. Differential-based approach

Artificial compression: Harten CPAM 1977, Olsson &

Kreiss JCP 2005

Anti-diffusion: So, Hu & Adams JCP 2011 2. 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 No Lagrangian moving grid or volume tracking cut cells

Differential-based interface sharpening

Use modified volume-fraction transport model as basis

∂tα + ~u · ∇α = 1 µI

Dα

µ_I ∈ R+: free parameter, D_α: interface-sharpening operator Compression form (Olsson & Kreiss JCP 2005)

Dα:= ∇ · [(εc∇α · ~n − α (1 − α)) ~n]

~n = ∇α/k∇αk, εc ∈ R+ (order of mesh size) Anti-diffusionform (So, Hu & Adams JCP 2011)

Dα :=−∇ · (εd∇α) ε_d∈ R^N+ (order of velocity magnitude)

Solution of interface-sharpening model

∂tα + ~u · ∇α = 1 µI

Dα

can be approximated by employing fractional step method That is, in each time step,

1. Transport step over time step ∆t for

∂tα + ~u · ∇α = 0 by state-of-the-art interface-capturing solver

2. Interface-sharpening stepover pseudo-time step ∆τ for

∂τα = Dα, τ = t/µI

by explicit finite-difference solver, for example

Passive tracer transport: Anti-diffusion results

Passive advection Rotating disc

Vortex in cell

Passive tracer transport: THINC results

Passive advection Rotating disc

Vortex in cell

Vortex-in-cell: Comparisons of results

With THINC

With anti-diffusion

Without interface sharpening

THINC-based interface sharpening

In THINC method, originalvolume-fraction equation

∂tα + ~u · ∇α = 0 is used in each time step as

1. Reconstruct piecewise smooth function ˜αi(x, tn) based on THINC reconstruction procedure from cell average {αⁿ_i} at time t_n

2. Constructspatial discretization for ~u · α using interpolated initial data from { ˜αi(x, tn)} obtained in step 1

3. Employsemi-discrete method to update αⁿ from current time to next αⁿ⁺¹ over time step ∆t

THINC reconstruction in step 1 assumes

˜

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

2



1 +γtanh



βx − xi−1/2

∆x −x¯i



αM = M (α_i−1, α_i+1) , M := min, max, γ = sgn (α_i+1− α_i−1) β measures sharpeness (given constant) & x¯i chosen to fulfill

α_i = 1

∆x

Z x_i+1/2 x_i−1/2

˜

α_i(x)dx

For example, cell edges used in step 2 determined by α_i+1/2,L = αmin+ αmax− αmin

2



1 + γ tanh β +C 1 +Ctanh β



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

2 (1 + γC) (C not shown)

THINC reconstruction: Graphical view

-0.5 0 0.5 1 1.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

β=2 β=3

(x − xi+1/2)/(xi+1/2− xi−1/2)

˜αi(x)

αi−1/2,R

α_i+1/2,L

Passive tracer transport: Convergence study

Convergence study of 1-norm errors E1(α) as mesh is refined;

results for passive advectionare shown only

With THINC No sharpening With anti-diffu N × N E1(α) Order E1(α) Order E1(α) Order 50 × 50 9.8840 NaN 91.7486 NaN 4.0436 NaN 100 × 100 5.1746 0.93 60.6698 0.60 2.0558 0.98 200 × 200 2.6455 0.97 39.3623 0.62 0.9921 1.05 400 × 400 1.3373 0.98 25.2699 0.64 0.4414 1.17

Passive tracer transport: CPU timing study

CPU timing in seconds for Passive tracer transport problems (run on HP xw 9400 with AMD Dural-Core Opteron)

Method/Problem Passive advec. Rotating disk Vortex in cell

With THINC 20.5 21.8 280.7

With anti-df 32.1 29.4 383.5

No sharpening 33.2 28.8 344.9

Grid 100 × 100 100 × 100 200 × 200

This validates THINC & anti-diffusion schemes for passive tracer transport

Toy problem: Summary & remarks

1. Anti-diffusioninterface sharpening Employed aspost-processing step

Diffusion coefficientε_d & parameter µ_I controls interface sharpeness

Toy problem: Summary & remarks

1. Anti-diffusioninterface sharpening Employed aspost-processing step

Diffusion coefficientε_d & parameter µ_I controls interface sharpeness

2. THINCinterface sharpening

Employed aspre-processing step

Parameterβ controls interface sharpeness

Toy problem: Summary & remarks

1. Anti-diffusioninterface sharpening Employed aspost-processing step

Diffusion coefficientε_d & parameter µ_I controls interface sharpeness

2. THINCinterface sharpening

Employed aspre-processing step

Parameterβ controls interface sharpeness 3. Hybrid anti-diffusion-THINC is feasible

Toy problem: Summary & remarks

1. Anti-diffusioninterface sharpening Employed aspost-processing step

Diffusion coefficientε_d & parameter µ_I controls interface sharpeness

2. THINCinterface sharpening

Employed aspre-processing step

Parameterβ controls interface sharpeness 3. Hybrid anti-diffusion-THINC is feasible

4. For problems withmore than 2 fluid components,

interface sharpening based onanti-diffusion appears to be more robust than THINC

Compressible 1-phase flow (inviscid)

Assume inviscid, non-heat conducting, 1-phase, compressible flow in Cartesian coordinates:

∂tq + XN

j=1

∂xjfj(q) = 0

with q & fj, j = 1, 2, . . . , N, defined by q = (ρ, ρu1, . . . , ρuN, E)^T

fj = (ρuj, ρu1uj+ pδj1, . . . , ρuNuj+ pδjN, Euj+ puj)^T

Assume Mie-Gr¨uneisen equation of state

p(ρ, e) = p∞(ρ) + Γ(ρ)ρ [e − e∞(ρ)]

Sod Riemann problem: Exact solution

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Density

Rarefaction→

contact→

←Shock

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Pressure (Velocity)

Rarefaction→ ←Shock

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

t = 0 t=0.15

Density

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

t = 0 t=0.15

Pressure

Sod Riemann problem: Interface capturing solution

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Density

Rarefaction→

contact→

←Shock

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Pressure (Velocity)

Rarefaction→ ←Shock

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

Numerical Exact

Density

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

Numerical Exact

Pressure

Sod Riemann problem: THINC solution

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Density

Rarefaction→

contact→

←Shock

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Pressure (Velocity)

Rarefaction→ ←Shock

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

THINC Exact

Density

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

THINC Exact

Pressure

Sod Riemann problem: Anti-diffusion solution

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Density

Rarefaction→

contact→

←Shock

0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2

Pressure (Velocity)

Rarefaction→ ←Shock

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

adf Exact

Density

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

adf Exact

Pressure

How interface sharpening ? Compressible flow

Novel elements (as compared to toy problem for incompressible flow):

1. Robust local interface indicator Physical principlesbased

i.e., Flag interface cells by checking jumps of physical quantities nearby (good in 1D, less effective in 2D) Tracer transportbased

i.e., Flag interface cells based on classical volume-fraction approach (more effective in 2D) 2. Consistent interfacesolution reconstruction

(algebraic-based) or post-sharpening (differential-based) to ρ,ρ~u,E, · · ·

Interface reconstruction: Compressible 1-phase flow

Assume equilibrium pressure p & velocity ~u, motion of interface (contact discontinuity) is governed by

∂ρ

∂t + ~u · ∇ρ = 0

~u

∂ρ

∂t + ~u · ∇ρ



= 0

~u · ~u 2

∂ρ

∂t + ~u · ∇ρ

 +

∂

∂t(ρe) + ~u · ∇ (ρe)



= 0

1. Density ρ is reconstructed by basic THINC scheme 2. Momentum ρ~uis reconstructed by ~u times density in

part 1 (see second equation)

3. Total energy E is reconstructed by corrections on total kinetic & internal energy (see third equation)

THINC interface sharpening: Compressible flow

In each time step, our THINC-based interface-sharpening algorithm for compressible flow consists:

1. Reconstruct piecewise polynomialq˜i(x, tn) based on MUSCL/WENO reconstruction procedure from cell average {Qⁿ} at time tn

2. Modify ˜qi(x, tn) for interface cells using variant of THINC scheme from Qⁿ

3. Construct spatial discretization using interpolated initial data from {˜qi(x, tn)} obtained in steps 1 & 2

4. Employ semi-discrete method to update Qⁿ from current time to next Qⁿ⁺¹ over time step ∆t

Anti-diffusion: Compressible flow

Anti-diffusion model for compressible 1-phaseflow is

∂_tρ + div (ρ~u) = 1 µI

D_ρ

∂t(ρ~u) + div (ρ~u ⊗ u) + ∇p = 1 µI

Dρ~u

∂tE + div (E~u + p~u) = 1 µI

DE

Dρ:= −HI∇ · (εd∇ρ)

Dρu := u Dρ, DE :=hu · u

2 + ∂ρ(ρe)i Dρ

HI : Interface indicator =

 1 if interface cell, 0 otherwise

Anti-diffusion method

Anti-diffusion model in compact form

∂tq + divf (q) = 1 µI

Dq

with q, f, &D_q defined from Euler equations accordingly In each time step, fractional step is used

1. Solvehomogeneous equationwithout source terms

∂tq + divf (q) = 0 by state-of-the-art solver

2. Interface-sharpening stepover pseudo-time

∂τq = Dα, τ = t/µI

by explicit solver, for example

Compressible 2-phase flow: 7-equation model

7-equation non-equilibrium model of Baer & Nunziato (1986)

∂t(αρ)₁+ div (αρ~u)₁ = 0

∂t(αρ~u)₁+ div (αρ~u ⊗ ~u)₁+ ∇(αp)1 =pI∇α1+λ(~u2− ~u1)

∂_t(αE)₁+ div (αE~u + αp~u)₁ =

p_I~u_I · ∇α1−νp_I(p1− p2) +λ~u_I · (~u2− ~u1)

∂t(αρ)2+ div (αρ~u)2 = 0

∂t(αρ~u)2+ div (αρ~u ⊗ ~u)2+ ∇(αp)2 = −pI∇α1−λ(~u2− ~u1)

∂t(αE)₂+ div (αE~u + αp~u)₂ =

−pI~uI · ∇α1+νpI(p1 − p2) −λ~uI· (~u2− ~u1)

∂tα1+~uI · ∇α1 =ν(p1− p2) Saturation condition α1 + α2 = 1 Equation of state p_k(ρ_k, e_k), k = 1, 2

pI & ~uI: interfacial pressure & velocity

Baer & Nunziato (1986): pI = p2, ~uI = ~u1

Saurel & Abgrall (JCP 1999, JCP 2003) pI = α1p1+ α2p2, ~uI = α1ρ1~u1+ α2ρ2~u2

α1ρ1+ α2ρ2

pI = p1/Z1+ p2/Z2

1/Z1+ 1/Z2

, ~uI = ~u1Z1+ ~u2Z2

Z1+ Z2

, Zk = ρkck

ν = SI

Z1+ Z2

, λ = SIZ1Z2

Z1+ Z2

, SI(Interfacial area) ν & λ: relaxation parameters that express rates at which pressure & velocity toward equilibrium respectively

This non-equilibrium model can be used to simulate 1. Mixtures with different phasic pressures, velocities,

temperatures 2. Material interfaces

3. Permeable interfaces (i.e., interfaces separating a cloud of dispersed phases such as liquid drops or gases)

4. Cavitation if it is modeled as a simplified mechanical relaxation process, occurring at infinite rate µ → ∞ &

not modeled as a mass transfer process

Reduced 5-equation model

Kapila et al. 2001, Murrone et al. 2005, & Saurel et al. 2008 showed in asymptotic limits of λ & ν → ∞ i.e., flow towards mechanical equilibrium: ~u1 = ~u2 = ~u&p1 = p2 = p,

7-equation model reduces to 5-equation model

∂t(α1ρ1) + div (α1ρ1~u) = 0

∂t(α2ρ2) + div (α2ρ2~u) = 0

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

∂tE + div (E~u + p~u) = 0

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

 K2− K1

K1/α1+ K2/α2



div(~u), Ki = ρic²_i

Reduced 5-equation model

Kapila et al. 2001, Murrone et al. 2005, & Saurel et al. 2008 showed in asymptotic limits of λ & ν → ∞ i.e., flow towards mechanical equilibrium: ~u1 = ~u2 = ~u&p1 = p2 = p,

7-equation model reduces to 5-equation model

∂t(α1ρ1) + div (α1ρ1~u) = 0

∂t(α2ρ2) + div (α2ρ2~u) = 0

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

∂tE + div (E~u + p~u) = 0

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

 K2− K1

K1/α1+ K2/α2



div(~u), Ki = ρic²_i Equilibrium (mixture) pressure psatisfies

p = ρe − X2

k=1

αkρke∞,k(ρk) + X2 k=1

αk

p∞,k(ρk) Γk(ρk)

! X²

k=1

αk

Γk(ρk)

Reduced 5-equation model is hyperbolicwith non-monotonic equilibrium sound speed ceq:

1

ρc²_eq = α1

ρ1c²1

+ α2

ρ2c²2

(Wood’s formula) yielding stiffness in equations & numerical solver

0 0.2 0.4 0.6 0.8 1

10² 10³

α_water ceq

2-phase (air-water)

5 -equation transport model

For nearly single-phase flow, where α1 ≈ 0 or 1, Allaire, Clerc,

& Kokh (JCP 2002) proposed using

∂t(α1ρ1) + div (α1ρ1~u) = 0

∂t(α2ρ2) + div (α2ρ2~u) = 0

∂t(ρ~u) + div (ρ~u ⊗ ~u) + ∇p = 0

∂tE + div (E~u + p~u) = 0

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

Mixture pressure p computed in same manner as before Phasic entropy Sk satisfies relation

∂p1

∂S1



ρ1

DS1

Dt −

∂p2

∂S2



ρ2

DS2

Dt = ρ1c²1− ρ2c²2

div(~u) 6= 0

Model is hyperbolic, but with monotone frozen sound speed

ρc²_frozen = X2 k=1

αkρkc²_k

0 0.2 0.4 0.6 0.8 1

10¹ 10² 10³ 10⁴

µ−equil frozen

αwater

ceq&cfrozen

Interface reconstruction: 5-equation model

Assume equilibrium pressure p, velocity ~u, &phasic density ρk

for each interface cell

1. Volume fraction α1 is reconstructed by basic THINC scheme, denoted byα˜1

2. Reconstructphasic & mixture density α_iρ_i & ρ by g

αiρi =α˜iρi, i = 1, 2, ρ˜=α˜1ρ1+ (1 −α˜1)ρ2

3. Reconstructmomentum ρ~u by fρ~u = ˜ρ~u

4. Reconstructtotal energy E by E˜ = 1

2ρ~u · ~u +˜ α˜1ρ1e1 + (1 −α˜1)ρ2e2

5 -equation transport model: Anti-diffusion

5-equation transport model with anti-diffusion

∂t(α1ρ1) + div (α1ρ1~u) = 1 µI

Dα1ρ1

∂t(α2ρ2) + div (α2ρ2~u) = 1 µ_IDα2ρ2

∂_t(ρu) + div (ρ~u ⊗ ~u) + ∇p = 1 µI

D_ρu

∂tE + div (E~u + p~u) = 1 µI

DE

∂_tα1 + ~u · ∇α1 = 1 µI

D_α1

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

1. α-ρ based (Shyue 2011)

Dα1 :=−∇ · (εd∇α1), Dα_kρ_k :=−HI∇ · (εd∇αkρk) D_ρ:=

X2 k=1

D_αkρk, D_ρ~u := ~uD_ρ, K = 1 2~u · ~u DE := KDρ+

X2 k=1

∂α_kρ_k(ρkek)Dα_kρ_k + X2 k=1

ρkekDα_k

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

Dα1 :=−∇ · (εd∇α1), Dαkρk :=ρkDαk, Dα2 := −Dα1

Dρ :=

X2 k=1

Dα_kρ_k, Dρ~u := ~uDρ, DE := KDρ+ X2 k=1

ρkekDα_k

Dα1 := ∇ · [(εc∇α1· ~n − α1(1 − α1)) ~n] applicable

Shock(morb)-contact(moly) interaction

With THINC Without THINC

With THINC Pressure

Without THINC Pressure

With THINC Without THINC

With THINC Pressure

Without THINC Pressure

5 -equation model: Axisymmetric case

Axisymmetric version of reduced 5-equation model reads

∂t(A(x)α1ρ1) + ∂x(A(x)α1ρ1u) = 0

∂t(A(x)α2ρ2) + ∂x(A(x)α2ρ2u) = 0

∂t(A(x)ρu) + ∂x A(x)ρu²

+ A(x)∂xp = 0

∂t(A(x)E) + ∂x(A(x)Eu + A(x)pu) = 0

∂tα1+ ∂x(α1u) = α1K¯ K1

∂xu +

α1K¯ K1

− α1

A^′(x) A(x)∂xu

where 1/ ¯K =P²

i=1αi/Ki, Ki = ρic²_i

Spherical UNDEX (Wardlaw)

−1 −0.5 0 0.5 1 −1

−0.5 0

0.5 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

water explosive

Spherical UNDEX: Initial phase

0 5 10 15 20

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

t=5.859e−6s t=13.12e−6s t=21.90e−6s t= 32.36e−6s t=44.64e−6s t=58.81e−6s

Density

Spherical UNDEX: Shock-contact interaction phase

0 10 20 30 40 50 60

0 0.2 0.4 0.6 0.8 1 1.2 1.4

t=93.05e−6s t=135.9e−6s t=188.4e−6s t=252.0e−6s

Spherical UNDEX: Incompressible phase

0 20 40 60 80 100 120 140 160

0 0.2 0.4 0.6 0.8 1 1.2 1.4

t=420.0e−6s t=528.8e−6s t=658.1e−6s t=811.4e−6s

Spherical UNDEX: bubble collapse & rebound

0 20 40 60 80 100

10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

x

t=0.2970e−1s t=0.3004e−1s t=0.3046e−1s t=0.3122e−1s

Spherical UNDEX test: Diagnosis

0 5 10 15 20 25 30

0 5 10 15 20 25 30 35 40 45 50 55

with THINC with antiD

Bubble radius (cm)

time (ms)

0 5 10 15 20 25 30

10⁵ 10⁶ 10⁷ 10⁸ 10⁹ 10¹⁰ 10¹¹

with THINC with antiD

Interface pressure (dyn)

time (ms)

Spherical UNDEX test: Diagnosis

Table: Quantitative study of maximum bubble radius rmax &

period of bubble oscillation T_b for spherical underwater explosion rmax(cm) error (%) Tb(ms) error (%)

Experiment 48.10 0 29.8 0

Incompressible 66.49 38.2 39.1 31.2

Luo et al. 48.75 1.4 29.7 0.3

Wardlaw 46.40 3.5 29.8 0

THINC 48.17 0.1 29.1 2.3

Anti-diffusion 48.57 0.1 29.5 1.1

Underwater explosions: cylindrical wall

High pressure gaseous explosive in water (cylindrical case)

water

explosive

t = 30µs

shock

reflected rarefaction reflected shock

t = 45µs cavitation

interface

t = 60µs

cavitation

t = 120µs

cavitation

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

6 -equation model: Alternative to 5-equation model

Non-equilibrium 6-equation model of Saurel et al. (JCP 2009):

∂_t(α1ρ1) + div (α1ρ1~u) = 0

∂_t(α2ρ2) + div (α2ρ2~u) = 0

∂t(ρ~u) + div (ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0

∂t(α1ρ1e1) + div (α1ρ1e1~u) + α1p1∇ · ~u =−pIν (p1− p2)

∂t(α2ρ2e2) + div (α2ρ2e2~u) + α2p2∇ · ~u =pIν (p1− p2)

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

Model is hyperbolic with monotone frozen sound speed & is equivalent to reduced 5-equation model asymptotically as ν → ∞ (i.e., p1 → p2)

Pelanti & Shyue (2012) proposed

∂t(α1ρ1) + div(α1ρ1~u) = 0

∂t(α2ρ2) + div(α2ρ2~u) = 0

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0

∂t(α1E1) + div (α1E1~u + α1p1~u) +B (q, ∇q)=−νpI (p1− p2)

∂t(α2E2) + div (α2E2~u + α2p2~u) −B (q, ∇q)=νpI (p1− p2)

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

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

Use phasic total energy instead of phasic internal energy;

numerically easy to retain mixture total energy consistency

Interface reconstruction: 6-equation model

Assume equilibrium pressure p, velocity ~u, &phasic density ρ_k for each interface cell again

1. Reconstructvolume fraction α1, phasic density α_iρ_i, total density ρ, momentum ρ~u, in same manner as for

5-equation model

2. Reconstructphasic total energy αiEi by αgiEi = 1

2α˜iρi~u · ~u +α˜iρiei

6 -equation model: Anti-diffusion

∂t(α1ρ1) + div(α1ρ1~u) = 1 µI

Dα1ρ1

∂_t(α2ρ2) + div(α2ρ2~u) = 1 µI

D_α2ρ2

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 1 µI

Dρu

∂t(α1E1) + div (α1E1~u + α1p1~u) + B (q, ∇q) =

− νpI(p1− p2) + 1 µI

Dα1E1

∂t(α2E2) + div (α2E2~u + α2p2~u) − B (q, ∇q) = νpI (p1− p2) + 1

µI

Dα2E2

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

D_α1

Write 6-equation model in compact form as

∂tq + div(f (q)) + w(q, ∇q) = ψµ_I(q) + ψν(q)

Compute approximate solution based on fractional step:

1. Homogeneous hyperbolic step

∂tq + div(f (q)) + w(q, ∇q) = 0 2. Source-terms relaxationstep

∂tq = ψµ_I(q) + ψν(q)

Pressure relaxation ν → ∞

Continue by considering pressure relaxationwith

∂tq = ψν(q), as ν → ∞ Current ODE system is

∂t(αρ)1 = 0

∂t(αρ~u)1 = 0

∂t(αρE)₁ =−νpI(p1− p2)

∂t(αρ)₂ = 0

∂t(αρ~u)₂ = 0

∂t(αρE)₂ =νpI (p1− p2)

∂tα1 =ν (p1− p2)

Combining energy & volume fraction equations, we have

∂t(αρE)_k= −pI∂tαk

yielding (after using mass & momentum equations)

∂t(αρe)_k = −pI∂tαk or αkρk∂tek = −pI∂tαk

Integration wrt t over ∆t, we have αkρk(ek− ek0) = −

Z αk

αk0

pIdαk = −¯pI(αk− αk0) or

e_k = e_k0− ¯p_I

 1 ρk

− 1 ρk0



, k = 1, 2

Combining EOS ek(¯pI, ρk)with ek = ek0− ¯pI

 1 ρ_k − 1

ρ_k0

 , we find phasic density ρk as a function of ¯pI, i.e.,

ρk(¯pI) , k = 1, 2 Use saturation condition

1 = α1ρ1

ρ1(¯p_I) + α2ρ2

ρ2(¯p_I)

leads to algebraic equation for p¯I (relaxed pressure)

Having relaxed p¯I = p1 = p2 & so ρk, αk, conservative vector q can be updated (EOS should be imposed)

Future perspective I

Consider 6-equation model with heat & mass transfer of form

∂t(α1ρ1) + div(α1ρ1~u) =m˙

∂t(α2ρ2) + div(α2ρ2~u) = −m˙

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0

∂t(α1E1) + div (α1E1~u + α1p1~u) + B (q, ∇q) =

− νpI(p1− p2) +Q + eIm˙

∂t(α2E2) + div (α2E2~u + α2p2~u) − B (q, ∇q) = νpI(p1− p2)−Q − eIm˙

∂tα1+ ~u · ∇α1 = µ (p1− p2) + m˙ ρ_I

Assume µ, θ, ν → ∞: instantaneous relaxation effects 1. Volumetransfer via pressure relaxation: µ (p1− p2)

µexpresses rate toward mechanical equilibrium p1 → p2,

& isnonzero in all flow regimes of interest 2. Heattransfer via temperature relaxation:

Q := θ (T2− T1)

θexpresses rate towards thermal equilibrium T1→ T2, 3. Mass transfer viathermo-chemical relaxation:

˙

m := ν (g2− g1)

ν expresses rate towards diffusive equilibrium g1 → g2, &

isnonzero only at 2-phase mixture& metastable state, i.e.,

ν =

∞ ǫ1≤ α1 ≤ 1 − ǫ1 &T_liquid> T_sat

0 otherwise

Assume µ, θ, ν → ∞: instantaneous relaxation effects 1. Volumetransfer via pressure relaxation: µ (p1− p2)

µexpresses rate toward mechanical equilibrium p1 → p2,

& isnonzero in all flow regimes of interest 2. Heattransfer via temperature relaxation:

Q := θ (T2− T1)

θexpresses rate towards thermal equilibrium T1→ T2, 3. Mass transfer viathermo-chemical relaxation:

˙

m := ν (g2− g1)

ν expresses rate towards diffusive equilibrium g1 → g2, &

isnonzero only at 2-phase mixture& metastable state, i.e.,

ν =

∞ ǫ1≤ α1 ≤ 1 − ǫ1 &T_liquid> T_sat

0 otherwise

Liquid-vapor phase change depends strongly on numerical resolution of αk

Dodecane 2-phase Riemann problem

Saurel et al. (JFM 2008) & Zein et al. (JCP 2010):

Liquid phase: Left-hand side (0 ≤ x ≤ 0.75m)

(ρ_v, ρ_l, u, p, α_v)_L = 2kg/m³, 500kg/m³, 0, 10⁸Pa, 10⁻⁸
Vapor phase: Right-hand side (0.75m < x ≤ 1m)

(ρv, ρl, u, p, αv)_R = 2kg/m³, 500kg/m³, 0, 10⁵Pa, 1 − 10⁻⁸
,

Liquid Vapor

← Membrane

Dodecane 2-phase problem: Phase diagram

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³

10¹ 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷ 10⁸

Liquid Vapor

2-phase mixture

Saturation curve →

v

p

Dodecane 2-phase problem: Phase diagram

Wave path in p-v phase diagram

10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹ 10² 10³

10¹ 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷ 10⁸

Liquid Vapor

2-phase mixture

Saturation curve →

v

p

Isentrope

Hugoniot locus

Dodecane 2-phase problem: Sample solution

0 0.2 0.4 0.6 0.8 1

10⁰ 10¹ 10² 10³

t = 0 t=473µs

Density (kg/m³)

0 0.2 0.4 0.6 0.8 1

−50 0 50 100 150 200 250 300 350

t = 0 t=473µs

Velocity (m/s)

0 0.2 0.4 0.6 0.8 1

10⁴ 10⁵ 10⁶ 10⁷ 10⁸

t = 0 t=473µs

Pressure(bar)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

t = 0 t=473µs

Vapor volume fraction

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

t = 0 t=473µs

Vapor mass fraction

4-wave structure:

Rarefaction, phase, contact, &

shock

Dodecane 2-phase problem: Sample solution

Allphysical quantities arediscon- tinuous across phase boundary

Future perspective II

Consider barotropic1-pressure, 1-velocitycompressible 2-phase flow model with drift flux approximation

∂t(α1ρ1) + div(α1ρ1~u) = div (ρY1Y2~uR) (Continuity α1ρ1)

∂t(α2ρ2) + div(α2ρ2~u) = −div (ρY1Y2~uR) (Continuity α2ρ2)

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇p= 0 (Momentum)

Future perspective II

Consider barotropic1-pressure, 1-velocitycompressible 2-phase flow model with drift flux approximation

∂t(α1ρ1) + div(α1ρ1~u) = div (ρY1Y2~uR) (Continuity α1ρ1)

∂t(α2ρ2) + div(α2ρ2~u) = −div (ρY1Y2~uR) (Continuity α2ρ2)

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇p= 0 (Momentum) Equilibrum pressure p computed by solving

α1+ α2 = α1ρ1

ρ1(p) + α1ρ1

ρ2(p) = 1

Future perspective II

Consider barotropic1-pressure, 1-velocitycompressible 2-phase flow model with drift flux approximation

∂t(α1ρ1) + div(α1ρ1~u) = div (ρY1Y2~uR) (Continuity α1ρ1)

∂t(α2ρ2) + div(α2ρ2~u) = −div (ρY1Y2~uR) (Continuity α2ρ2)

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇p= 0 (Momentum) Equilibrum pressure p computed by solving

α1+ α2 = α1ρ1

ρ1(p) + α1ρ1

ρ2(p) = 1 Darcy law model for relative velocity~u_R assumes

~u_R= 1 λα1α2

ρ2− ρ1

ρ



∇p

Future perspective II

Consider barotropic1-pressure, 1-velocitycompressible 2-phase flow model with drift flux approximation

∂t(α1ρ1) + div(α1ρ1~u) = div (ρY1Y2~uR) (Continuity α1ρ1)

∂t(α2ρ2) + div(α2ρ2~u) = −div (ρY1Y2~uR) (Continuity α2ρ2)

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇p= 0 (Momentum) Equilibrum pressure p computed by solving

α1+ α2 = α1ρ1

ρ1(p) + α1ρ1

ρ2(p) = 1 Darcy law model for relative velocity~u_R assumes

~u_R= 1 λα1α2

ρ2− ρ1

ρ



∇p Accurate resolution of dissipative source terms& so mathematical model requires good approximation of αk

Rather than solving saturation condition for p, we may consider model that includes volume fraction equation explicitly as

∂t(α1ρ1) + div(α1ρ1~u) = div (ρY1Y2~uR)

∂t(α2ρ2) + div(α2ρ2~u) = −div (ρY1Y2~uR)

∂t(ρ~u) + div(ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0

∂_tα1 + ~u · ∇α1− ρc²_w α1α2

ρ1c²1ρ2c²2

ρ2c²2− ρ1c²1

div(~u) =

ρc²_w α1α2

ρ1c²1ρ2c²2

ρ1c²1

α1ρ1

+ ρ2c²2

α2ρ2



div(ρY1Y2~uR)

Here cw is Wood sound speed defined by 1

ρc²_w = α1

ρ1c²1

+ α2

ρ2c²2

(Wood’s formula)

Water wave problem

Existence or non-existence of 2-bump solution in water wave via computer aided proof in 2-phase (air-water) direct numerical simulation

Future perspective III

Hyperelasticity flow · · ·

Future perspective III

Hyperelasticity flow · · ·

Thank you