Diffuse interface models & methods for

Diffuse interface models & methods for

compressible multiphase flows

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

May 8-12, 2017, CSRC, Beijing, China

What is multiphase flow ?

General features: Existence dispersed phases (fluid-like or not) in underlying flow separated by interfaces, e.g.,

Examples from C.E. Brennen: Fundamentals of

Multiphase Flow, Cambridge, 2005

Examples from C.E. Brennen: Fundamentals of

Multiphase Flow, Cambridge, 2005

Examples from C.E. Brennen: Fundamentals of

Multiphase Flow, Cambridge, 2005

Pressure wave in bubbly liquid

Boiling flow (Saurel etal CAF 2016)

Bubble collapse over flat surface

Condensation fronts

Fracturing of solid plate (Ndanou etal JCP 2015)

Diffuse interface models & methods ?

Mathematical models: Average (fluid-mixture) type Numerical methods: Discontinuity-capturing type

Diffuse interface models & methods ?

Mathematical models: Average (fluid-mixture) type Numerical methods: Discontinuity-capturing type Sharp interface models & methods (Prof. Xiaolin Li’s talk)

Mathematical models: Single-phase model with jump conditions across interfaces

Numerical methods: Front tracking

Diffuse interface models & methods ?

Mathematical models: Average (fluid-mixture) type Numerical methods: Discontinuity-capturing type Sharp interface models & methods (Prof. Xiaolin Li’s talk)

Mathematical models: Single-phase model with jump conditions across interfaces

Numerical methods: Front tracking References

D. Jamet: Diffuse interface models in fluid mechanics (technical report)

D. M. Anderson, G. B. McFadden, & A. A. Wheeler:

Diffuse-interface methods in fluid mechanics, Annu. Rev.

Fluid Mech., 1998, 32:139-65

Lecture outline

Single phase → Two phase → Three phase & more

Lecture outline

Single phase → Two phase → Three phase & more

That is,

1 & 2. Fluid-mixture type models & methods for compressible multicomponent flow

3 & 4. Homogeneous relaxation models & methods methods for compressible two-phase flow

5 & 6. Homogeneous relaxation models & methods for compressible three-phase flow & more

Joint work with

F. Xiao (Tokyo Tech., Japan): Eulerian interface-sharpening

M. Pelanti (ENSTA, Paris Tech., France): Two- &

three-phase models & algorithms

Fluid-mixture type models & methods for

compressible multicomponent flow

I: Linear equation of state

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

May 8-12, 2017, CSRC, Beijing, China

Compressible flow model: Single-phase case

Assume continuum mechanics modelling of fluid motion, basic physical principles for pure phase are

1. Mass balance:

∂ρ

∂t + ∇ · (ρ⃗u) = 0 2. Momentum balance:

∂(ρ⃗u)

∂t + ∇ · (ρ⃗u ⊗ ⃗u) = ∇ · T 3. Energy balance:

∂(ρE)

∂t + ∇ · (ρE⃗u) = ∇ ·! T ⃗u"

− ∇ · ⃗q 4. Entropy balance:

∂(ρs)

∂t + ∇ · (ρs⃗u) = −∇ · ⃗qs+ ∆s

Assume Newtonian fluid &Fourier’s, we have ρ: density, ⃗u : flow velocity

E : specific total energy E = e +1

2⃗u · ⃗u, e : specific internal energy T : stress tensor

T = −pI + τ, p: pressure, τ : dissipative stress tensor τ = µ

#

∇⃗u + ∇^T⃗u − 2

3(∇ · ⃗u) I

$

, µ: viscosity

⃗q: heat flux

⃗q = −κ∇T, κ: thermal conductivity, T : temperature s: specific entropy

⃗qs : entropy flux, ⃗qs = ⃗q T

∆s : entropy source, ∆s= κ∇T · ∇T

T² + µ∇⃗u : ∇⃗u T

Thermodynamic laws

Second law of thermodynamics:

ds = dQ

T , Q: heat transfer, T : temperature First law of thermodynamics:

dQ = de + pdV, V = 1/ρ: specific volume This leads to

T ds = de + pdV =⇒ de = T ds − pdV yielding definiton of p & T , i.e.,

p = −# ∂e

∂V

$

s

, T = # ∂e

∂s

$

V

when expression for e(V, s) is known a priori

Equation of state

In compressible fluid flow theory, equation of state such as p = p(ρ, e), p = p(ρ, T ), or p = p(ρ, s)

is commonly used to characterize thermodynamic behavior of underlying medium

If equation of state depends on density only, p = p(ρ), it is callled isentropic(or barotropic) equation

Here with fixed entropy s, equation of state is assumed to be monotonically nondecreasing & convex, i.e.,

# ∂p

∂ρ

$

s

≥ 0, # ∂²p

∂ρ²

$

s

≥ 0 Define speed of sound, denoted by c ∈ R, as

c² =# ∂p

∂ρ

$

s

≥ 0; c² = 0 only if ρ = 0 (vacuum)

Equation of state: Ideal gas

Gas medium is assumed to be ideal gas, if it satisfies laws of Boyle & Gay-Lussac as expressed by equation of state

pV = RT

R = nR0 (n: number of gas moles & R0: univ. gas constant) In ideal gas, internal energy e is function of temperature T only, see Courant-Friedrichs p.8-9

If gas is polytropic, we have

e(T ) = CVT, CV: specific heat at constant volume p(ρ, e) = (γ − 1)ρe, γ: ratio of specific heats

p(ρ, s) = (γ − 1) exp# s − s0

CV

$ ρ^γ e(ρ, s) = exp# s − s0

CV

$ ρ^(γ−1)

Balance equations: Remark

When balance equations are closed, entropy balanceequation

& energy balanceequation are equivalent

In this case, entropy balance equation can be written as ρCp

DT

Dt = ∇ · (κ∇T ) + κTTDp

Dt + τ : ∇⃗u

where Cp & κT denotes specific heat at constant pressure &

coefficient of thermal expansion defined by Cp = T

(∂T /∂s)p

& κT = −ρ# ∂s

∂p

$

T

= −1 ρ

# ∂ρ

∂T

$

p

Recall double dot product of two tensors U & V is defined by U : V = %

i

%

j

UijVji

D/Dt denotes material derivative

Compressible flow: Model summary

System of balance equations for compressible flow is 1. Mass balance:

∂ρ

∂t + ∇ · (ρ⃗u) = 0 2. Momentum balance:

∂(ρ⃗u)

∂t + ∇ · (ρ⃗u ⊗ ⃗u) + ∇! pI"

= ∇ · τ 3. Energy balance^∗:

∂(ρE)

∂t + ∇ ·!

ρE⃗u + pI⃗u"

= ∇ · (τ⃗u) + ∇ · (κ∇T ) 4. Entropy balance^∗:

ρCp

dT

dt = ∇ · (κ∇T ) + κTTdp

dt + τ : ∇⃗u Take only one equation from ∗ & model closes with EOS

Inviscid compressible flow: Riemann problem

For compressible Euler equations in 1D, Riemann problem is Cauchy problem that consists of

∂w

∂t +∂f (w)

∂x = 0, x ∈ R, t > 0 with

w =

⎡

⎣ ρ ρu ρE

⎤

⎦, f (w) =

⎡

⎣ ρu ρu²+ p ρEu + pu

⎤

⎦ as for model equations, & piece-wise constant data

w(x, 0) =

*wL if x < 0 wR if x > 0 as for initial condition

Riemann problem: Hyperbolicity

To close model & Riemann problem, assume ideal gas law p = (γ − 1)ρe

Jacobian matrix of f (w), denoted by A, is

A = ∂f (w)

∂w =

⎡

⎣

0 1 0

(γ − 3)u²/2 −(γ − 1)u γ − 1 (γ − 1)u³/2 − Hu H − (γ − 1)u² γu

⎤

⎦ Its eigen-decomposition AR = RΛ, is with

Λ = diag(u − c, u, u + c)

R =

⎡

⎣

1 1 1

u − c u u + c H − uc ¹₂u² H + uc

⎤

⎦

c = +γp/ρ > 0is speed of sound & H = (e + p)/ρ is specific enthalpy

Riemann problem: Basic solution structure

Elementary waves of Riemann problem in x-t plane

x t

wL wR

rarefaction contact

shock

Riemann problem: Remarks

Basic properties of elementary wavesare 1. Shock wave

Genuinely nonlinear wave

Solution is discontinuous across wave that follows Rankine-Hugoniotjump condition

σ[w] = [f (w)], σ: shock speed & [w] = w_R− w_L 2. Contact discontinuity

Linearly degenerate wave

Solution is inmechanical equilibriumacross wave, i.e., [u] = 0,[p] = 0, while typically [ρ] ̸= 0, [T ] ̸= 0, [s] ̸= 0 3. Rarefaction wave

Genuinely nonlinear wave

Riemann invariantis constant across wave

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

Fluid-mixture type methods

Basic elements:

1. Use diffuse interface model of fluid-mixture type

2. Employ state-of-the-art finite-volume method tocapture discontinuities (shocks & contacts) implicitly

3. Underlying grid may be either static & uniform or time-dependent & nonuniform

Focus on fluid-mixture model first, numerical discretization of model will be discussed later (when time permits)

Sod Riemann problem: High-resolution result

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 result

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 result

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

Extension to two-fluid flow: Model problem

Consider shock-tube problem with immiscible membrane separating two different ideal gases, say, air& He

Assume mechanical equilibrium condition for states across interface, i.e.,

kinematic: [u] = 0 & dynamic: [p] = 0

Breaking of membrane would result in shock/rarafaction waves separated by moving interface

ideal gas 1 (γ₁) ideal gas 2 (γ₂)

Fluid-mixture methods: Numerical issues

1. Oncontinuous level:

Equations of motions & equation of state for fluid mixtures ?

2. Ondiscrete level:

Consistent approximation of model system ?

ideal gas 1 (γ1) ideal gas 2 (γ2) Mixing region (γ =?)

Benchmark test: 2-fluid interface only

Initial Condition

⎛

⎜

⎜

⎝ ρ u p γ

⎞

⎟

⎟

⎠

L

=

⎛

⎜

⎜

⎝ 1 1 1 1.4

⎞

⎟

⎟

⎠

&

⎛

⎜

⎜

⎝ ρ u p γ

⎞

⎟

⎟

⎠

R

=

⎛

⎜

⎜

⎝ 0.125

1 1 1.2

⎞

⎟

⎟

⎠ Exact solution

⎛

⎜

⎜

⎝ ρ u p γ

⎞

⎟

⎟

⎠

(x, t) =

⎛

⎜

⎜

⎝

ρ0(x − t) 1 1 γ0(x − t)

⎞

⎟

⎟

⎠

Model equations ?

1. Basic conservation laws for ρ, ρu, & ρE (fluid mixtures) 2. Pressure computed from

(a) Equation of state p = (γ − 1)[ρE − (ρu)²/(2ρ)] with (i) ∂γ

∂t + u ∂γ

∂x = 0 (ii) ∂

∂t

# 1

γ − 1

$ + u ∂

∂x

# 1

γ − 1

$

= 0 (or others)

(b) Pressure evolution equation directly (Karni 1996)

∂p

∂t + u ∂p

∂x + γp∂u

∂x = 0

Interface only problem: fluid-mixture results

2.54.05.5 2.54.05.5 2.54.05.5 2.54.05.5

1.01.2 1.01.2 1.01.2 1.01.2

1.001.06 1.001.06 1.001.06 1.001.06

0.0 0.4 0.8

1.201.301.40

0.0 0.4 0.8

1.201.301.40

0.0 0.4 0.8

1.201.301.40

0.0 0.4 0.8

1.201.301.40

with ργ with ρ/(γ − 1) with γ with 1/(γ − 1)

x x

x x

ρe

ρe

ρe

ρe uuuu pppp γγγγ

Fluid-mixture model: Ideal gas

Abgrall (1996): Single-density version in multi-D 1. Model (quasi-conservative) system

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·!

ρ⃗u ⊗ ⃗u + pI"

= 0

∂

∂t(ρE) + ∇ ·!

ρE⃗u + pI⃗u"

= 0

∂

∂t

# 1 γ − 1

$

+ ⃗u · ∇

# 1

γ − 1

$

= 0 2. Equation of state

p(ρ, e,γ) = (γ− 1)ρe

Transport equation for 1/(γ − 1) plays role near interfaces only

Fluid-mixture model: Mapped grid version

1. Model system

∂ρ

∂t + 1 J∇ξ·!

ρ⃗U"

= 0

∂

∂t(ρ⃗u) + 1 J∇ξ·!

ρ⃗u ⊗ ⃗U + pIJ"

= 0

∂

∂t(ρE) + 1 J∇ξ·!

ρE ⃗U + pI ⃗U"

= 0

∂

∂t

# 1

γ − 1

$ + 1

JU · ∇⃗ ξ

# 1

γ − 1

$

= 0 2. Equation of state

p(ρ, e, γ) = (γ − 1) ρe Uj =2N

i=1ui∂xiξj: contravariant velocity in ξj-direction

Mapped grid model: Grid metrics

Basic coordinate mapping relations in N = 3 are

⎛

⎜

⎜

⎝

1 0 0 0

0 ∂x1ξ1 ∂x2ξ1 ∂x3ξ1

0 ∂x1ξ2 ∂x2ξ2 ∂x3ξ2

0 ∂x1ξ3 ∂x2ξ3 ∂x3ξ3

⎞

⎟

⎟

⎠

= 1 J

⎛

⎜

⎜

⎝

J 0 0 0

0 J11 J21 J31

0 J12 J22 J32

0 J13 J23 J33

⎞

⎟

⎟

⎠ where J = |∂(x1, x2, x3)/∂(ξ1, ξ2, ξ3)| =

det (∂(x1, x2, x3)/∂(ξ1, ξ2, ξ3)), J11 =

3 3 3 3

∂(x2, x3)

∂(ξ2, ξ3) 3 3 3 3

, J21= 3 3 3 3

∂(x1, x3)

∂(ξ3, ξ2) 3 3 3 3

, J31 = 3 3 3 3

∂(x1, x2)

∂(ξ2, ξ3) 3 3 3 3 , J12 =

3 3 3 3

∂(x2, x3)

∂(ξ₃, ξ₁) 3 3 3 3

, J22= 3 3 3 3

∂(x1, x3)

∂(ξ₁, ξ₃) 3 3 3 3

, J32 = 3 3 3 3

∂(x1, x2)

∂(ξ₃, ξ₁) 3 3 3 3 , J13 =

3 3 3 3

∂(x2, x3)

∂(ξ₁, ξ₂) 3 3 3 3

, J23= 3 3 3 3

∂(x1, x3)

∂(ξ₂, ξ₁) 3 3 3 3

, J33 = 3 3 3 3

∂(x1, x2)

∂(ξ₁, ξ₂) 3 3 3 3 .

Moving cylindrical vessel

Impose uniform flow velocity (u₁, u₂) = (−1, 0) (i.e., in the frame of vessel moving with speed one in x1-direction)

air helium

interface

Moving cylindrical vessel

Solution at time t = 0.25

Moving cylindrical vessel

Solution at time t = 0.5

Moving cylindrical vessel

Solution at time t = 0.75

Moving cylindrical vessel

Solution at time t = 1

Extension to stiffened gas

One simple equation of state that models materials not only gases, but also compressible liquids & solids is stiffened gas

p(ρ, e) = (γ − 1) ρe +(ρ − ρ0) B

In this case, with basic conservation laws, diffuse interface model includes additional transport equations as

∂

∂t

# 1

γ − 1

$

+ ⃗u · ∇

# 1

γ − 1

$

= 0

∂

∂t

# ρ0B γ − 1

$

+ ⃗u · ∇# ρ0B γ − 1

$

= 0

∂

∂t

# ρB γ − 1

$ + ∇ ·

# ρB γ − 1⃗u

$

= 0

Transport equations for interfaces

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 · ∇ρ

$ +4 ∂

∂t(ρe) + ⃗u · ∇ (ρe) 5

= 0 This is,

∂ρ

∂t + ⃗u · ∇ρ = 0

∂

∂t(ρe) + ⃗u · ∇ (ρe) = 0

Effective transport equations: Stiffened gas

With stiffened gas EOS, tranport equation for ρe of interface

∂

∂t(ρe) + ⃗u · ∇ (ρe) = 0 takes form

∂

∂t

# p

γ − 1− ρ − ρ0

γ − 1B

$

+ ⃗u · ∇

# p

γ − 1− ρ − ρ0

γ − 1B

$

= 0 To ensure pressure in equilibrium, split equation into two parts

∂

∂t

# p

γ − 1

$

+ ⃗u · ∇

# p

γ − 1

$

= 0

∂

∂t

# ρ − ρ0

γ − 1B

$

+ ⃗u · ∇# ρ − ρ0

γ − 1B

$

= 0

That is, to retain pressure equilibrium across interface, we must have

∂

∂t

# 1

γ − 1

$

+ ⃗u · ∇

# 1

γ − 1

$

= 0

∂

∂t

# ρ − ρ0

γ − 1B

$

+ ⃗u · ∇# ρ − ρ0

γ − 1B

$

= 0

However, to ensure conservation of mass in region away from interface, latter equation needs to split further & rewrite as

∂

∂t

# ρ0B γ − 1

$

+ ⃗u · ∇# ρ0B γ − 1

$

= 0

∂

∂t

# ρB γ − 1

$ + ∇ ·

# ρB γ − 1⃗u

$

= 0

Note that from

∂

∂t

# ρB γ − 1

$

+ ⃗u · ∇

# ρB γ − 1

$

= 0

& by employing chain rule with respect to ρ, we have

∂

∂t

# ρB γ − 1

$

+ ⃗u · ∇

# ρB γ − 1

$

= B

γ − 1

∂ρ

∂t + B

γ − 1⃗u · ∇ρ

= B

γ − 1

# ∂ρ

∂t + ⃗u · ∇ρ

$

= B

γ − 1

# ∂ρ

∂t + ⃗u · ∇ρ

$

= − ρB γ − 1∇⃗u yielding

∂

∂t

# ρB γ − 1

$ + ∇ ·

# ρB γ − 1⃗u

$

= 0

Fluid-mixture model: Stiffened gas

Saurel, Abgrall (SISC 1999), Shyue (JCP 1998, 1999)

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·!

ρ⃗u ⊗ ⃗u + pI"

= 0

∂

∂t(ρE) + ∇ ·!

ρE⃗u + pI⃗u"

= 0

∂

∂t

# 1

γ − 1

$

+ ⃗u · ∇

# 1

γ − 1

$

= 0

∂

∂t

# ρ0B γ − 1

$

+ ⃗u · ∇# ρ0B γ − 1

$

= 0

∂

∂t

# ρB γ − 1

$ + ∇ ·

# ρB γ − 1⃗u

$

= 0 p(ρ, e, γ, B) = (γ − 1) ρe + (ρ − ρ0) B

Fluid-mixture model: Volume-fraction version I

Suppose there exists Mf ≥ 1different immiscible phases initially where each of them occupies distinct portion with volume fraction αk ∈ [0, 1], k = 1, 2, . . . , Mf,2Mf

k=1αk = 1 Define mixture states for ρ & ρe as

ρ =

Mf

%

k=1

αkρk, ρe =

Mf

%

k=1

αkρkek

Suppose pressure p is in equilibrium& constitutive law for each fluid phase is characterized by stiffened gas, we then have

ρe =

Mf

%

k=1

αk

# p

γk− 1− ρk− ρ0,k

γk− 1 Bk

$

= p

γ − 1 −ρ − ρ0

γ − 1B

It follows

1 γ − 1 =

Mf

%

k=1

αk

# 1

γk− 1

$

ρ0B γ − 1 =

M_f

%

k=1

αk

# ρ0,kBk

γk− 1

$

ρB γ − 1 =

Mf

%

k=1

αk

# ρkBk

γk− 1

$

Note that rather than using 1/(γ − 1) & ρ0B/(γ − 1) as state variables in our model system, instead we may use αk as

∂αk

∂t + ⃗u · ∇αk= 0, k = 1, 2, . . . , Mf

This leads to volume-fraction based model of form

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·!

ρ⃗u ⊗ ⃗u + pI"

= 0

∂

∂t(ρE) + ∇ ·!

ρE⃗u + pI⃗u"

= 0

∂

∂t

# ρB γ − 1

$ + ∇ ·

# ρB γ − 1⃗u

$

= 0

∂αk

∂t + ⃗u · ∇αk = 0, k = 1, 2, . . . , Mf

If partial density αkρk is used as state variables, rather than mixture density ρ, this leads to so called 5-equation transport model (Allaire et al. , JCP 2002), if Mf = 2, that is,

Fluid-mixture model: Volume-fraction version II

∂

∂t (αkρk) + ∇ · (αkρk⃗u) = 0, k = 1, 2, . . . , Mf

∂

∂t(ρ⃗u) + ∇ ·!

ρ⃗u ⊗ ⃗u + pI"

= 0

∂

∂t(ρE) + ∇ ·!

ρE⃗u + pI⃗u"

= 0

∂αk

∂t + ⃗u · ∇αk= 0, k = 1, 2, . . . , Mf

This gives 2Mf + 2 equations in total for Mf-fluid problem Pressure is computed directly by solving

ρe =

Mf

%

k=1

αkρkek(ρk, p)

which is easy to do if EOS can be written in Mie-Gr¨uneisen form (see next lecture)

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Liquid-falling problem

Underwater explosion with two circular obstacles

Underwater explosion with two circular obstacles

Underwater explosion with two circular obstacles

Shock in water over dispersed gas/solid in

cylindrical nozzle

Shock in water over dispersed gas/solid in

cylindrical nozzle

Shock in water over dispersed gas/solid in

cylindrical nozzle

Shock in water over dispersed gas/solid in

cylindrical nozzle

Shock in water over dispersed gas/solid in

cylindrical nozzle

Shock in water over dispersed gas/solid in

cylindrical nozzle

Fluid-mixture type models & methods for

compressible multicomponent flow

II: Nonlinear equation of state

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

May 8-12, 2017, CSRC, Beijing, China

Outline

1. Model problems 2. Model nonlinear EOS

3. Single-phase flow model: Nonlinear EOS 4. Multi-fluid flow model: Nonlinear EOS 5. Slip line problem

Motivation: Flying aluminum-plate problem

vacuum

Al target

Al flyer

time = 0.5µs time = 1µs time = 2µs time = 4µs

Motivation: 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 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)

Nonlinear EOS for real materials

To deal with problems with real materials, consider Mie-Gr¨uneisen type equation of state of form

p(V, e) = pref(V ) + Γ(V )

V !e − eref(V )" , Typical examples are

1. van der WaalsEOS for real gas Γ(V ) = γ − 1

1 − b/V , (V = 1/ρ) pref(V ) = aV⁻²

eref(V ) = −aV⁻¹

a: molecular cohesive forces, b: finite size of molecules, γ: ratio of specific heats

Mie-Gr¨uneisen EOS: Jones-Wilkins-Lee case

2. Jones-Wilkins-Lee(JWL) EOS for gaseous explosives

Γ(V ) = γ − 1 eref(V ) = e⁰+A V0

R1

exp# −R1V V0

$

+B V0

R2

exp# −R2V V0

$

pref(V ) = p⁰+ A exp# −R1V V0

$

+ B exp# −R2V V0

$

Reference state lies along an isentrope where deref/dV = pref

Material-dependent parameters: γ, A, B, R1,R2 V0

Mie-Gr¨uneisen EOS: Cochran-Chan case

3. Cochran-Chan (CC) EOS for solid explosives Γ(V ) = γ − 1

eref(V ) = e⁰+−AV0

1 − E1

%# V V0

$^1−E1

−1

&

+ BV0

1 − E2

%# V V0

$^1−E₂

−1

&

pref(V ) = p⁰+ A# V V0

$^−E₁

− B# V V0

$^−E₂

Reference state lies along an isentrope where deref/dV = pref

Material-dependent parameters: γ, A, B, E1,E2 V0

Mie-Gr¨uneisen EOS: Shock-Hugoniot case

4. Hugoniot-based EOS for solid metals based on linear shock-particle velocity relation

Γ(V ) = Γ0

# V V0

$α

pref(V ) = p⁰+ c²0(V0−V ) [V0−s(V0−V )]² eref(V ) = e⁰+1

2!pref(V ) + p⁰" (V0−V ) Reference state lies along a Hugoniot locus where Rankine-Hugoniot jump conditions are satisfied Material-dependent parameters: Γ0,V0,c0,s, α

Non-Mie-Gr¨uneisen EOS

1. Two molecular (oxygen-nitrogen) vibrating gas p(ρ, e) = ρRT (e)

ρe = C_V^trT + ρ RΘvib

exp(Θvib/T ) − 1 e.g., R = 287.086J/(Kg · K), C_V^tr = R/(γtr−1), γtr = 1.4, Θvib = 10³K

Non-Mie-Gr¨uneisen EOS

2. Osborne model p(ρ, e) = 1

ρ0e + φ0

[ζ(a0+ a1ζ) + ρ0e(b0+ ζ(b1+ b2ζ)+

ρ0e(c0+ c1ζ))]

ζ = V0

V −1

e.g., for water ρ0 = 10⁻², φ0 = 2 × 10⁻²,

⎡

⎣

a0 a1 a2

b0 b1 b2

c0 c1 c2

⎤

⎦=

⎡

⎣

3.84 × 10⁻⁴ 1.756 × 10⁻³ 0 1.312 × 10⁻² 6.265 × 10⁻² 2.133 × 10⁻¹ 5.132 × 10⁻¹ 6.761 × 10⁻¹ 0

⎤

⎦

Mie-Gr¨uneisen EOS: Material parameters

JWL EOS ρ0(kg/m³) A(GPa) B(GPa) R1 R2 Γ

TNT1 1630 371.2 3.23 4.15 0.95 0.30

TNT2 1630 548.4 9.375 4.94 1.21 1.28

Water 1004 1582 −4.67 8.94 1.45 1.17

CC EOS ρ0(kg/m³) A(GPa) B(GPa) E1 E2 Γ

TNT 1840 12.87 13.42 4.1 3.1 0.93

Copper 8900 145.67 147.75 2.99 1.99 2

Shock EOS ρ0(kg/m³) c0(m/s) s Γ0 α

Aluminum 2785 5328 1.338 2.0 1

Copper 8924 3910 1.51 1.96 1

Thermodynamic stability

Define fundamental derivative of gas dynamics G = −V

2

(∂²p/∂V²)s

(∂p/∂V )s

, s : specific entropy Assume fluid state satisfy G> 0 forthermodynamic stability, i.e.,

(∂²p/∂V²)s> 0 & (∂p/∂V )s< 0 (∂²p/∂V²)s> 0 meansconvex EOS

(∂p/∂V )s < 0 means real speed of sound, for c² =# ∂p

∂ρ

$

s

= −V²# ∂p

∂V

$

s

> 0

In practice, G= (K_s^′ + 1)/2with Ks = ρc² &

K_s^′ =# ∂K_s

∂p

$

s

= +

ρ∂p

∂ρ + ρ²∂²p

∂ρ² + ∂p

∂ρ

∂p

∂e + 2p ∂²p

∂ρ∂e− p

ρ

∂p

∂e + p ρ²

# ∂p

∂e

$² +p²

ρ²

∂²p

∂e²

&

, Ks

G ≤0 case, i.e., with existence of nonclassical nonlinear wave& phase transition (see next page), will not be discussed here

Phase diagram: van der Waals EOS

0.5 1 1.5 2 2.5 3

0.6 0.7 0.8 0.9 1 1.1 1.2

G = 0 G < 0

Saturation Curve Liquid−Vapor Coexistence Region

G > 0

C T = 1

T > 1

Vapor Region Liquid

Region

p/pc

V /Vc

Nonlinear EOS: Fundamental numerical challenge

Consider interface only problem for single-phase flow, &

assume fluid is gas governed by van der Waals EOS ρe =# 1 − bρ

γ − 1

$

-p + aρ². − aρ²

To maintain pressure equilibrium numerically, scheme should be consistent approximation of

∂

∂t

+# 1 − bρ γ − 1

$

-p + aρ². − aρ² /

+ u ∂

∂x

+# 1 − bρ γ − 1

$

-p + aρ². − aρ² /

= 0 based on numerical approximation of ρ from

∂ρ

∂t + u∂ρ

∂x = 0

Most state-of-the-art shock-capturing scheme for hyperbolic conservation laws fails, yielding spurious oscillations in mechanical-equilibrium states near diffused interfaces e.g., Results extracted from Pantano et al. (JCP 2017)

To correct this, one viable approach is to consider splitting

∂

∂t

+# 1 − bρ γ − 1

$ p

/ + u ∂

∂x

+# 1 − bρ γ − 1

$ p

/

= 0

∂

∂t

+# 1 − bρ γ − 1

$

aρ²−aρ² /

+ u ∂

∂x

+# 1 − bρ γ − 1

$

aρ²−aρ² /

= 0 To retain pressure in equilibrium, we find

∂

∂t

# 1 − bρ γ − 1

$ + u ∂

∂x

# 1 − bρ γ − 1

$

= 0

∂

∂t

+# 1 − bρ γ − 1

$

aρ²−aρ² /

+ u ∂

∂x

+# 1 − bρ γ − 1

$

aρ²−aρ² /

= 0

For single-phase flow where γ & b, are constants

∂

∂t

# 1 − bρ γ − 1

$ + u ∂

∂x

# 1 − bρ γ − 1

$

= 0 reduce to

∂ρ

∂t + u∂ρ

∂x = 0 hence it can be ignored

Thus to retain pressure equilibrium across interfaces, & also for states across shock/rarefaction waves, one introduces

∂

∂t

# 2 − γ − bρ γ − 1 aρ²

$ + ∂

∂x

# 2 − γ − bρ γ − 1 aρ²u

$ +

# 2 − γ − 2bρ

γ − 1 aρ²$ ∂u

∂x = 0

Single-phase flow model: van der Waals EOS

1. Mathematical model

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·0

ρ⃗u ⊗ ⃗u + pI1

= 0

∂

∂t(ρE) + ∇ ·0

ρE⃗u + pI⃗u1

= 0

∂

∂t

# 2 − γ − bρ γ − 1 aρ²

$

+ ∇ ·# 2 − γ − bρ γ − 1 aρ²⃗u

$ +

# 2 − γ − 2bρ γ − 1 aρ²

$

∇ ·⃗u = 0 2. Compute pressure from EOS

p =# γ − 1 1 − bρ

$ + ρe−

#2 − γ − bρ γ − 1 aρ²

$/

Extension to multi-fluid flow: van der Waals EOS

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·0

ρ⃗u ⊗ ⃗u + pI1

= 0

∂

∂t(ρE) + ∇ ·0

ρE⃗u + pI⃗u1

= 0

∂

∂t

# 1 γ − 1

$

+ ⃗u · ∇

# 1 γ − 1

$

= 0

∂

∂t

# bρ γ − 1

$ + ∇ ·

# bρ γ − 1⃗u

$

= 0

∂

∂t

# 2 − γ − bρ γ − 1 aρ²

$

+ ∇ ·# 2 − γ − bρ γ − 1 aρ²⃗u

$ +

# 2 − γ − 2bρ γ − 1 aρ²

$

∇ ·⃗u = 0

Extension to multi-fluid flow: Mie-Gr¨uneisen EOS

∂ρ

∂t + ∇ · (ρ⃗u) = 0

∂

∂t(ρ⃗u) + ∇ ·0

ρ⃗u ⊗ ⃗u + pI1

= 0

∂

∂t(ρE) + ∇ ·0

ρE⃗u + pI⃗u1

= 0

∂

∂t

# 1 Γ

$

+ ⃗u · ∇# 1 Γ

$

−ρ# 1 Γ

$^′

∇ ·⃗u = 0

∂

∂t 0pref

Γ

1+ ⃗u · ∇0pref Γ

1

−ρ0pref Γ

1^′

∇ ·⃗u = 0

∂

∂t-ρeref. + ⃗u · ∇ -ρeref. + ρ -ρeref.^′∇ ·⃗u = 0

Fluid-mixture model: Hyperbolicity

Write our model system in quasi-linear form as

∂w

∂t + A(w) ∂w

∂x = 0 with state vector w and the matrix A defined by

w = +

ρ, ρu, ρE, 1 Γ, pref

Γ , ρeref /T

A =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

0 1 0 0 0 0

K − u² u(2 − Γ) Γ −pΓ −Γ Γ

u(K − H) H − u²Γ u(Γ + 1) −upΓ −uΓ uΓ

−ϕu ϕ 0 u 0 0

−χu χ 0 0 u 0

−ψu ψ 0 0 0 u

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

Eigen-structure of matrix A is Eigenvalues

Λ = diag(λ1, λ2, · · · , λ6) = diag(u−c, u, u+c, · · · , u) Eigenvectors

R = (r1, r2, · · · , r6) =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

1 1 1 0 0 0

u − c u u + c 0 0 0

H − uc u²/2 H + uc p 1 −1

ϕ 0 ϕ 1 0 0

χ 0 χ 0 1 0

ψ 0 ψ 0 0 1

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦ with Ar_k = λ_kr_k. Here K = Γu²/2, H = E + (p/ρ),

ϕ = −Γ^′/Γ², χ = (Γp^′ref−Γ^′pref)/Γ², & ψ = eref + ρe^′ref

Shock molybdenum over MORB (Mid-Ocean Ridge Basalt) liquid

density pressure

Molybdenum

MORB time = 50µs

time = 100µs

Shock over MORB: AMR result

density pressure

Flying aluminum in water over copper half space

density pressure

Water

Aluminum Copper

time = 50µs

time = 100µs

time = 150µs

Slip line (shear flow) problem

Consider plane interface moving vertically from the left to right in x1-direction as shown inN = 2 below. Interface conditions for this problem are

Dynamic condition: pR= pL

Kinematic condition: u1,R = u1,L & (u2,R−u2,L) ̸= 0

slip line

¯ u1

¯

u2 −¯u2

Slip line problem: Model equations

In this case, transport equations for motion of ρ, ρu2, &

ρe + ρu²2/2 are, respectively,

∂ρ

∂t + ¯u1

∂ρ

∂x1

= 0 (from mass conservation law)

∂

∂t (ρu2) + ¯u1

∂

∂x1

(ρu2) = 0 (from momentum in x2-dir)

∂

∂t (ρe) + ¯u1

∂

∂x1

(ρe) +

∂

∂t

# 1 2ρu²2

$ + ¯u1

∂

∂x1

# 1 2ρu²2

$

= 0 (from total energy)

Assume single phase ideal gasflow with EOS p(ρ, e) = (γ − 1)ρe heats

In this instance, to ensure pressure equilibrium, as it should be for this slip line problem, from

∂

∂t

# p

γ − 1

$ + ¯u1

∂

∂x1

# p

γ − 1

$ +

∂

∂t

# 1 2ρu²2

$ + ¯u1

∂

∂x1

# 1 2ρu²2

$

= 0, we should have

∂

∂t

# 1 2ρu²2

$ + ¯u1

∂

∂x1

# 1 2ρu²2

$

= 0

Note that when the problem is solved numerically by an interface-capturing method, it is common to compute pressure, p, based on conservative variables as

p = (γ − 1) 4

E− 52

i=1(ρui)² 2ρ

6 ,

while generally (ρu2)²/2ρ̸=ρu²2/2when a slip line is smeared out

Immediate consequence of this is loss of pressure equilibrium, and so incorrect solution of other state variables (see below)

Slip line problem: Godunov method

Suppose first order Godunov method of form Qⁿ⁺¹_ij = Qⁿ_ij − ∆t

∆x1

!F1-Qⁿ_i,j, Qⁿ_i+1,j. − F1-Qⁿ_i−1,j, Qⁿ_ij."

is used for numerical discretization of this slip line problem Here Qⁿ_ij denotes approximate value of cell average of solution q over cell C_ij at time t_n, i.e.,

Qⁿ_ij = 1

∆x1∆x2

7 7

Cij

q(x1, x2, tn)dx1dx2,

and F1(Qi±1,j, Qij) denotes numerical flux, say, defined by F1(Qi±1,j, Qij) = f1-q^∗_i±1/2,j.

(q_i±1/2,j^∗ Riemann solution)

Slip line problem: Godunov results

Errors depend strongly on transverse velocity jump

0 0.2 0.4 0.6 0.8 1

1 2 3 4 5 6

0 0.2 0.4 0.6 0.8 1

1 1.05 1.1 1.15 1.2 1.25 1.3 1.35

0 0.2 0.4 0.6 0.8 1

9.95 10 10.05 10.1 10.15

0 0.2 0.4 0.6 0.8 1

−1

−0.5 0 0.5 1

x1

x1

u1 u2

ρ p

Slip line problem: Remarks

To devise a more accurate method for numerical resolution of slip lines, we may use

Diffuse interface approach

Include transverse kinetic energy equation in the model

& use its solution for pressure update p = (γ − 1)

#

E− (ρu1)² 2ρ +ρu²₂

2

$

This transverse kinetic equation should be modified so that there is no difficulty to work with shock waves Sharp interface approach

Front tracking or Lagrangian moving grid method

Slip line problem: quasi Lagrangian results

Errors are on the order of machine epsilon

0 0.5 1 1.5

0 1 2 3 4 5 6

0 0.5 1 1.5

1 1 1 1 1

0 0.5 1 1.5

10 10 10 10 10 10 10 10

0 0.5 1 1.5

−1

−0.5 0 0.5 1

x1

x1

u1 u2

ρ p

Slip line problem: Lagrangian grid

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.05 0.1 0.15

Initial grid

0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2

0 0.05 0.1 0.15 0.2

Final grid

References: Fluid-mixture type algorithms

K.-M. Shyue: An efficient shock-capturing algorithm for compressible multicomponent problems, JCP 1998 K.-M. Shyue: A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state, JCP 1999

K.-M. Shyue: A fluid-mixture type algorithm for

compressible multicomponent llow with Mie-Grueneisen equation of state, JCP 2001

K.-M. Shyue: A fluid-mixture type algorithm for barotropic two-fluid flow problems, JCP 2004

K.-M. Shyue: A volume-fraction based algorithm for hybrid barotropic and non-barotropic two-fluid flow problems, Shock Waves 2006

Fluid-mixture type algorithm: Limitations

1. Model is not effectiveto deal problems with large difference in pressure such as cavitating flow

2. Nonequilibrium interface dynamics such asdrag, phase change, & heat conduction are difficult to incorporate into model

Need to consider general 2-phase flow model for solutions (One nature choice is ???)

5 -equation transport model: 2-phase case

Allaire et al. (JCP 2002) 1. Mathematical model

∂

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

∂

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

∂

∂t(ρ⃗u) + ∇ ·0

ρ⃗u ⊗ ⃗u + pI1

= 0

∂

∂t(ρE) + ∇ ·0

ρE⃗u + pI⃗u1

= 0

∂α1

∂t + ⃗u · ∇α1 = 0 2. Pressure is computed by

p = +

(ρE) −(ρu)²

2ρ +0pref Γ

1−- ρeref.

/ , # 1 Γ

$

Homogeneous relaxation models &

methods for

compressible two-phase flow

I: Reduced models

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

May 8-12, 2017, CSRC, Beijing, China

Outline

Main theme: Compressible 2-phase (liquid-gas) solver for metastable fluids: application to cavitation & flashing flows

1. Motivation

2. Constitutive law for metastable fluid

3. Homogeneous relaxation model (HRM) for compressible 2-phase flow with & without heat & mass transfer 4. Reduced5-, 4-, & 3-equation model of HRM

Flashing flow means a flow with dramatic evaporationof liquid due to pressure drop

Outline

Main theme: Compressible 2-phase (liquid-gas) solver for metastable fluids: application to cavitation & flashing flows

1. Motivation

2. Constitutive law for metastable fluid

3. Homogeneous relaxation model (HRM) for compressible 2-phase flow with & without heat & mass transfer 4. Reduced5-, 4-, & 3-equation model of HRM

Flashing flow means a flow with dramatic evaporationof liquid due to pressure drop

Numerical solver of relaxation type will be discussed in next lecture

Phase transition: Thermodynamic path I

Sample thermodynamic path for phase transition with non-convex EOS (require phase boundary modelling)

0.5 1 1.5 2 2.5 3

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3

Liquid Vapor

s = s0

Phase boundary

V

p

Phase transition: Thermodynamic path II

Figure extracted from Saurel et al. (JFM 2008)

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

Initial condition 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

Dodecane 2-phase problem: Phase diagram

Thermodynamic 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