A mixture-energy-consistent numerical method for compressible two-phase flow with interfaces, cavitation, and evaporation waves

A mixture-energy-consistent numerical method for compressible two-phase flow

with

interfaces, cavitation, and evaporation waves

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

May 30, 2014

Joint work with Marica Pelanti at ENSTA, Paris Tech, France

Outline

1. Model problems 2. Constitutive law

Saturation curve for metastable liquids 3. Compressible 2-phase flow model

Homogeneous relaxation model (HRM) Reduced 5-equation model

4. Numerical scheme

Finite volume method Stiff relaxation solver 5. Numerical examples

Cavitating Richtmyer-Meshkov problem

Gas volume fraction Mixture pressure

t=0ms

t=2ms

t=3.1ms

t=6.4ms

t=8.6ms

High-speed underwater projectile

With thermo-chemical relaxation No thermo-chemical relaxation

High-pressure fuel injector

With thermo-chemical relaxation No thermo-chemical relaxation

Constitutive law

Stiffened gas equation of state (SG EOS) with Pressure

pk(ek, ρk)= (γk−1)ρk(ek−ηk) − γkp∞k

Temperature

Tk(pk, ρk)= pk+ p∞k

(γk−1)Cvkρk

Entropy

sk(pk, Tk)= Cvklog Tkγ_k

(pk+ p∞k)^γk−1 + η_k^′ Helmholtz free energy ak = ek−Tksk

Gibbs free energy gk = ak+ pkvk, vk = 1/ρk

Constitutive law: SG EOS parameters

Ref: Le Metayer et al. , Intl J. Therm. Sci. 2004

Fluid Water

Parameters/Phase Liquid Vapor

γ 2.35 1.43

p^∞ (Pa) 10⁹ 0

η (J/kg) −11.6 × 10³ 2030 × 10³

η^′ (J/(kg · K)) 0 −23.4 × 10³

Cv (J/(kg · K)) 1816 1040

Fluid Dodecane

Parameters/Phase Liquid Vapor

γ 2.35 1.025

p∞ (Pa) 4 × 10⁸ 0

η (J/kg) −775.269 × 10³ −237.547 × 10³

η^′ (J/(kg · K)) 0 −24.4 × 10³

Cv (J/(kg · K)) 1077.7 1956.45

Constitutive law: Saturation curves

Assume two phases in diffusive equilibriumwith equal Gibbs free energies (g1 = g2), saturation curveforphase transitionsis G(p, T )=A+B

T +ClogT+Dlog(p+p∞1)−log(p+p∞2) = 0 A= Cp1−Cp2+ η^′₂−η^′₁

Cp2−Cv2

, B= η1 −η2

Cp2−Cv2

C = Cp2−Cp1

Cp2−Cv2

, D= Cp1−Cv1

Cp2−Cv2

Constitutive law: Saturation curves

Assume two phases in diffusive equilibriumwith equal Gibbs free energies (g1 = g2), saturation curveforphase transitionsis G(p, T )=A+B

T +ClogT+Dlog(p+p∞1)−log(p+p∞2) = 0 A= Cp1−Cp2+ η^′₂−η^′₁

Cp2−Cv2

, B= η1 −η2

Cp2−Cv2

C = Cp2−Cp1

Cp2−Cv2

, D= Cp1−Cv1

Cp2−Cv2

or, from dg1 = dg2, we get Clausius-Clapeyron equation dp(T )

dT = Lh

T (v2−v1) Lh = T (s₂−s₁): latent heat of vaporization

Constitutive law: Saturation curves (Cont.)

Saturation curves for water & dodecane in T ∈ [298, 500]K

300 350 400 450 500

0 5 10 15 20 25 30

water dodecane Pressure(bar)

300 350 400 450 500

0 500 1000 1500 2000 2500

water dodecane Latent heat of vaporization (kJ/kg)

300 350 400 450 500

10⁻⁴ 10⁻³ 10⁻²

water dodecane Liquid volume v1(m³/kg)

300 350 400 450 500

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

water dodecane Vapor volume v2(m³/kg)

Homogeneous relaxation model (HRM)

Consider HRM for 2-phase flow of form

∂t(α1ρ1) + ∇ · (α1ρ1~u) =ν (g2−g1)

∂t(α2ρ2) + ∇ · (α2ρ2~u) =ν (g1−g2)

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

∂t(α₁E₁) + ∇ · (α₁E₁~u + α₁p₁~u)+B (q, ∇q)= µpI(p2−p1) + θTI(T2 −T1) + νgI(g2−g1)

∂t(α2E2) + ∇ · (α2E2~u + α2p2~u)−B(q, ∇q)= µpI(p1−p2) + θTI(T1 −T2) + νgI(g1−g2)

∂tα1+ ~u · ∇α1 =µ (p1−p2) + νvI(g1 −g2) B(q, ∇q) is non-conservative product (q: state vector)

B= ~u · [Y1∇(α2p2) − Y2∇(α1p1)]

Homogeneous relaxation model (Cont.)

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

µexpresses rate toward mechanical equilibrium p₁ → p₂,

& isnonzero in all flow regimes of interest

2. Heattransfer via temperature relaxation: θ (T2−T1) θexpresses rate towards thermal equilibrium T₁ → T₂, &

isnonzero only at 2-phase mixture

3. Mass transfer viathermo-chemical relaxation: ν (g2−g1) ν expresses rate towards diffusive equilibrium g₁ → g₂, &

isnonzero only at 2-phase mixture& metastable state T_liquid> T_sat

Homogeneous relaxation model (Cont.)

HRM model in compact form

∂tq + ∇ · f (q) + w (q, ∇q) = ψµ(q) + ψθ(q) + ψν(q) where

q= [α₁, α₁ρ₁, α₂ρ₂, ρ~u, α₁E₁, α₂E₂, α₁]^T f = [α1ρ1~u, α2ρ2~u, ρ~u ⊗ ~u + (α1p1+ α2p2)IN,

α1(E1+ p1) ~u, α2(E2+ p2) ~u, 0]^T w= [0, 0, 0, B (q, ∇q) , −B (q, ∇q) , ~u · ∇α1]^T

ψµ= [0, 0, 0, µpI(p2−p1) , µpI(p1−p2) , µ (p1−p2)]^T ψθ = [0, 0, 0, θTI (T2−T1) , θTI(T1−T2) , 0]^T

ψν = [ν (g2−g1) , ν (g1−g2) , 0, νgI(g2 −g1) , νgI(g1 −g2) , νvI(g1−g2)]^T

Homogeneous relaxation model (Cont.)

Flow hierarchy in HRM: H. Lund (SIAP 2012)

HRM µ → ∞ θν → ∞ µθν → ∞

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

Homogeneous relaxation model (Cont.)

Consider HRM limits as µ → ∞, µθ → ∞, & µθν → ∞

HRM µ → ∞ θν → ∞ µθν → ∞

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

HRM: Reduced model as θ = ν = 0 & µ → ∞

Assume frozen thermal & chemical relaxation,HRMreduces to

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

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

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

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

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

∂tα₁+ ~u · ∇α₁ =µ (p₁−p₂)

HRM: Reduced model as θ = ν = 0 & µ → ∞

Assume frozen thermal & chemical relaxation,HRMreduces to

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

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

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

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

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

∂tα₁+ ~u · ∇α₁ =µ (p₁−p₂)

Take formal asymptotic expansion ansatz q = q⁰+ εq¹+ · · ·

Find equilibrium equation for q⁰ as µ = 1/ε → ∞ (ε → 0⁺)

Reduced model: Asymptotic analysis as µ → ∞

Define material derivative d dt = ∂

∂t + ~u · ∇

Rewrite energy & volume-fraction equations in form dp₁

dt + ρ1c²₁∇ ·~u = ρ₁c²_I1 α1

1

ε(p2−p1) dp2

dt + ρ2c²₂∇ ·~u = ρ2c²_I2 α2

1

ε(p1−p2) dα₁

dt = 1

ε(p₁−p₂)

Assume formal asymptotic expansion as α1 = α⁰₁+ εα¹₁+ · · ·

pk = p⁰_k+ εp¹_k+ · · · for k = 1, 2

Reduced model: Asymptotic analysis (Cont.)

We get d

dt p⁰₁+ εp¹₁+ · · ·

+ ρ1c²₁∇ ·~u = ρ₁c²_I1

α1

1 ε

 p⁰₂−p⁰₁

+ ε p¹₂−p¹₁

+ · · · d

dt p⁰₂+ εp¹₂+ · · ·

+ ρ2c²₂∇ ·~u = ρ2c²_I2

α2

1 ε

 p⁰₁−p⁰₂

+ ε p¹₁−p¹₂

+ · · · d

dt α₁⁰+ εα¹₁+ · · ·

= 1 ε

 p⁰₁−p⁰₂

+ ε p¹₁−p¹₂

+ · · ·

Reduced model: Asymptotic analysis (Cont.)

Collecting equal power of ε, we have O(1/ε)

p⁰₁ = p⁰₂ =p⁰ =⇒ p⁰_I = p⁰, c⁰_Ik² = c⁰_k² O(1)

dp⁰

dt + ρ⁰₁c⁰₁²∇ ·~u = ρ⁰₁c⁰₁²

α⁰₁ p¹₂−p¹₁
dp⁰

dt + ρ⁰₂c⁰₂²∇ ·~u = ρ⁰₂c⁰₂²

α⁰₂ p¹₁−p¹₂
dα⁰₁

dt = p¹₁−p¹₂

Reduced model: Asymptotic analysis (Cont.)

Subtracting former two equations, we find

ρ⁰₁c⁰₁² −ρ⁰₂c⁰₂²

∇ ·~u = ρ⁰₁c⁰₁²

α₁⁰ +ρ⁰₂c⁰₂² α⁰₂

!

p¹₂−p¹₁

i.e.,

dα⁰₁

dt = p¹₁−p¹₂ = ρ⁰₂c⁰₂² −ρ⁰₁c⁰₁² ρ⁰₁c⁰₁²/α⁰₁+ ρ⁰₂c⁰₂²/α⁰₂

!

∇ ·~u

Reduced model as θ = ν = 0 & µ → ∞ (Cont.)

Thus, as θ = 0, ν = 0 & µ → ∞ leading order approximation of RHM model becomes so-called reduced 5-equation model (e.g., Kapila et al. 2001, Murrone et al. 2005)

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

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

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

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

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

 ρ2c²₂−ρ1c²₁ ρ₁c²₁/α₁ + ρ₂c²₂/α₂



∇ ·~u

Mixture pressure p = p(ρe, ρ1, ρ2, α1)determined from algebraic equation (linear one with SG EOS)

ρe = α1ρ1e1(p, ρ1) + α2ρ2e2(p, ρ2)

p relaxation: Subcharacteristic condition

Mechanical equilibrium sound speed cp ≤ cf (frozen speed) 1

ρc²_p = X2 k=1

αk

ρkc²_k & ρc²_f = X2 k=1

αkρkc²_k

p relaxation: Subcharacteristic condition

Mechanical equilibrium sound speed cp ≤ cf (frozen speed) 1

ρc²_p = X2 k=1

αk

ρkc²_k & ρc²_f = X2 k=1

αkρkc²_k

0 0.2 0.4 0.6 0.8 1

10¹ 10² 10³ 10⁴

frozen p relax

αwater

cp&cf

Non-monotonic cp

leads tostiffness in equations &

difficulties in numerical solver, e.g., positivity- preserving in volume fraction

HRM: Model as ν = 0, µ → ∞ & θ → ∞

Assume frozen chemical relaxation ν = 0, HRM in mechanical-thermal limit as µ → ∞ & θ → ∞ reads (Saurel et al. 2008, Fl˚atten et al. 2010)

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

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

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

∂t(ρY1) + ∇ · (ρY1~u) = 0

Mechanical-thermal equilibrium speed of sound satisfies 1

ρc²_pT = 1 ρc²_p +T

 Γ₂

ρ2c²₂ − Γ₁ ρ1c²₁

2  1 α1ρ1cp1

+ 1

α2ρ2cp2



HRM: Limit model as z → ∞, z = µ, θ, & ν

As all relaxation parameters go to infinity; z → ∞, z = µ, θ,

& ν, limit system of HRM is homogeneous equilibrium model (HEM) that follows standard mixture Euler equation

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

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

∂tE + ∇ · (E~u + p~u) = 0 This gives local resolution at interface only Speed of sound satisfies

1

ρc²_{pT g} = 1 ρc²_p +T

"

α1ρ1

Cp1

ds1

dp

2

+α2ρ2

Cp2

ds2

dp

2#

Equilibrium speed of sound: Comparison

Sound speeds follow subcharacteristic condition cpT g ≤cpT ≤cp ≤cf

Limit of sound speed

αlimk→1cf = lim

αk→1cp = lim

αk→1cpT = ck, lim

αk→1cpT g 6= ck

0 0.2 0.4 0.6 0.8 1

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

frozen p relax pT relax pTg relax

α_water cpTg,cpT,cp&cf

5 -equation model: liquid-vapor phase transition

Modelling phase transitionin metastable liquids Saurel et al. (JFM 2008) proposed

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

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

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

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

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

K¯s

K_s¹∇ ·~u+ 1 qI

Q+ 1 ρI

˙ m K¯s=

α1

K_s¹ + α2

K_s²

⁻1

, K_s^ι = ριc²_ι qI =

K_s¹ α1

+ K_s² α2

   Γ1

α1

+Γ2

α2



, Q = θ(T2−T1) ρI =

K_s¹ α1

+K_s² α2

   c²₁ α1

+ c²₂ α2



, m = ν(g˙ 2−g1)

5 -equation model (Cont.)

Mathematically, 5-equation model approaches tosame relaxation limits as HRM, but is difficult to solve numerically to ensure solution to be feasible

5 -equation model (Cont.)

Mathematically, 5-equation model approaches tosame relaxation limits as HRM, but is difficult to solve numerically to ensure solution to be feasible

Saurel et al. (JCP 2009) & Zein et al. (JCP 2010) proposedHRM based on phasic internal energy as alternative way to solve 5-equation model

HRM: Phasic-internal-energy-based

HRM based on phasic internal energy formulation of Saurel et al. (JCP 2009) & Zein et al. (JCP 2010) is

∂t(α1ρ1) + ∇ · (α1ρ1~u) = ν (g2−g1)

∂t(α2ρ2) + ∇ · (α2ρ2~u) = ν (g1−g2)

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

∂t(α1ρ1e1) + ∇ · (α1ρ1e1~u) + α1p1∇ ·~u=

µpI(p2−p1) + θTI(T2 −T1) + νgI(g2−g1)

∂t(α2ρ2e2) + ∇ · (α2ρ2e2~u) + α2p2∇ ·~u=

µpI(p1−p2) + θTI(T1 −T2) + νgI(g1−g2)

∂tα1+ ~u · ∇α1 = µ (p1−p2) + νvI(g2 −g1) To ensure conservation of mixture total energy include

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

HRM: Phasic-total-energy-based

Numerically more advantageousto use HRM based on phasic-total-energy formulation than phasic-internal-energy one; for ease of devise mixture-energy-consistent discretization Pelanti & Shyue (JCP 2014), i.e.,

∂t(α1ρ1) + ∇ · (α1ρ1~u) = ν (g2−g1)

∂t(α2ρ2) + ∇ · (α2ρ2~u) = ν (g1−g2)

∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α₁p₁+ α₂p₂) = 0

∂t(α1E1) + ∇ · (α1E1~u + α1p1~u) + B (q, ∇q) = µpI(p2−p1) + θTI(T2 −T1) + νgI(g2−g1)

∂t(α2E2) + ∇ · (α2E2~u + α2p2~u) − B (q, ∇q) = µpI(p1−p2) + θTI(T1 −T2) + νgI(g1−g2)

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

Relaxation scheme

To find approximate solution of HRM, in each time step, fractional-step method is employed:

1. Non-stiff hyperbolic step

Solve hyperbolic system without relaxation sources

∂tq + ∇ · f (q) + w (q, ∇q) = 0 using state-of-the-art solver over time interval ∆t 2. Stiff relaxation step

Solve system of ordinary differential equations

∂tq = ψµ(q) + ψθ(q) + ψν(q) in various flow regimes under relaxation limits

Relaxation scheme (Cont.)

Definition (mixture-energy consistent)

(i) Mixture total energy conservation consistency E⁰ = E^0,C

where E⁰ = (α1E1)⁰+ (α2E2)⁰ (ii) Relaxed pressure consistency

e^0,C = α₁^∗e1



p^∗,(α1ρ1)⁰ α^∗₁



+ α^∗₂e2



p^∗,(α2ρ2)⁰ α₂^∗

 ,

where e^0,C = E^0,C− ^(ρ~u)_2ρ^0·(ρ~0^u)0

Method proposed here with phasic-total-energy formulation is mixture-energy consistent

Non-stiff hyperbolic step: Mapped grid method

Consider solution of model system

∂tq + ∇ · f (q) + w (q, ∇q) = 0 in 2D general non-rectangular geometry

Model in integral form over any control volume C is d

dt Z

C

q dΩ = − Z

∂C

f (q) · ~n ds − Z

C

w (q, ∇q) dΩ

where ~n is outward-pointing normal vector at boundary ∂C

C

∂C →

Hyperbolic step: Mapped grid (Cont.)

Then finite volume method on control volume C reads Qⁿ⁺¹ = Qⁿ− ∆t

M(C)

Ns

X

j=1

hjF˘j −∆tW^∗M(C)

Qⁿ :=R

Cq(z, tn) dz/M(C)

M(C)measure (area in 2D or volume in 3D) of C Ns number of sides

hj length of j-th side (in 2D) or area of cell edge (in 3D) measured in physical space

F˘j numerical approximation tonormal flux in average across j-th side of grid cell

W^∗ cell average of w in cell C

Hyperbolic step: Mapped grid (Cont.)

Assume mapped (i.e., logically rectangular) grid in 2D, we get Qⁿ⁺¹_ij = Qⁿ_ij − ∆t

κij∆ξ₁

F_i+^{1 1}

2,j−F_i−^{1 1}

2,j



−

∆t κij∆ξ2

F_i,j+^{2 1}

2

−F_i,j−^{2 1}

2

−∆tW_ij^∗∆ξ1∆ξ2

i − 1 i − 1

i

i j

j

j + 1 j + 1

Cˆij

Cij

ξ1

ξ2

mapping

∆ξ1

∆ξ2

logical grid physical grid

←−

x1 = x1(ξ1, ξ2) x₂ = x₂(ξ₁, ξ₂) x1

x2

κij = M(Cij)/∆ξ1∆ξ2

Mapped grid method: Wave propagation (Cont.)

Godunov-type in wave propagation form is Qⁿ⁺¹_ij = Qⁿ_ij− ∆t

κij∆ξ1

A⁺₁∆Q_i−¹

2,j+ A⁻₁∆Q_i+¹

2,j

−

∆t κij∆ξ2



A⁺₂∆Q_i,j−¹

2 + A⁻₂∆Q_i,j+¹

2



A⁺₁∆Q

A⁻₁∆Q

fluctuations^A+1∆Q_i−¹

2,j, A⁻

1∆Q_i+¹

2,j, A⁺₂∆Q_i,j−¹

2, &

A⁻₂∆Q_i,j+¹

2 : Solve

one-dimensional Riemann problemsin direction normal to cell edges

W_ij^∗ may be included in fluctuations

Mapped grid mtehod: Wave propagation (Cont.)

Speeds & limited of waves are used to calculatesecond order correction:

Qⁿ⁺¹_ij := Qⁿ⁺¹_ij − ∆t κij∆ξ1

Fe_i+^{1 1}

2,j− eF_i−^{1 1}

2,j

−

∆t κij∆ξ2

Fe_i,j+^{2 1}

2

− eF_i,j−^{2 1}

2



For example, at cell edge (i − ¹₂, j) correction flux takes

Fe_i−^{1 1}

2,j = 1 2

Nw

X

k=1

 λ^1,ki−¹₂,j

 1 − ∆t κ_i−¹

2,j∆ξ1

 λ^1,ki−¹₂,j

! Wf_i−^1,k1

2,j

κ_i−¹

2,j = (κi−1,j+ κij)/2, Wf_i−^1,k1

2,j is limited waves to avoid oscillations near discontinuities

Mapped grid method: Wave propagation (Cont.)

Transverse wave propagation is included to ensure second order accuracy & also improve stability

A⁻₂A⁺₁∆Q

A⁺₂A⁻₁∆Q

A⁺₂A⁺₁∆Q

A⁻₂A⁻₁∆Q

Mapped grid method: Wave propagation (Cont.)

Method can be shown to be quasi conservative & stable under a variant of CFL (Courant-Friedrichs-Lewy) condition

∆t max

i,j,k



  λ^1,k_i−¹

2,j

  κip,j∆ξ1

,  λ^2,k_i,j−¹

2

  κi,jp∆ξ2



 ≤ 1,

ip = i if λ^1,k_i−1

2,j > 0 & i − 1 if λ^1,k_i−1 2,j < 0

Mapped grid method: Wave propagation (Cont.)

Semi-discrete wave propagation method takes form

∂tQ(t) = L(Q(t)) where in 2D

L(Qij(t) = − 1 κij∆ξ1



A⁺₁∆Q_i−¹

2,j+ A⁻₁∆Q_i+¹

2,j+ A1∆Qij



− 1

κij∆ξ2



A⁺₂∆Q_i,j−¹

2 + A⁻₂∆Q_i,j+¹

2 + A2∆Qij



ODEs are integrated in time using strong stability-preserving (SSP) multistage Runge-Kutta, e.g., 3-stage 3rd-order

Q^∗ = Qⁿ+ ∆tL (Qⁿ) Q^∗∗ = 3

4Qⁿ+1

4Q^∗+1

4∆tL (Q^∗) Qⁿ⁺¹ = 1

3Qⁿ+2

3Q^∗+2

3∆tL (Q^∗∗)

Relaxation scheme: Stiff solvers

1. Algebraic-basedapproach

Saurel et al. (JFM 2008), Zein et al. (JCP 2010), LeMartelot et al. (JFM 2013), Pelanti-Shyue (JCP 2014) Imposeequilibrium conditions directly, without making explicit of interface states pI, gI,. . .

2. Differential-based approach

Saurel et al. (JFM 2008), Zein et al. (JCP 2010) Imposedifferential of equilibrium conditions, require explicit of interface states pI, gI,. . .

3. Optimization-based approach (for mass transfer only) Helluy & Seguin (ESAIM: M2AN 2006), Faccanoni et al.(ESAIM: M2AN 2012)

p relaxation: Algebraic approach

Assume frozen thermal & thermo-chemical relaxation, i.e., θ = 0 & ν = 0, look for solution of ODEs in limitµ → ∞

∂t(α1ρ1) = 0

∂t(α2ρ2) = 0

∂t(ρ~u) = 0

∂t(α1E1) = µpI(p2−p1)

∂t(α2E2) = µpI(p1−p2)

∂tα1 =µ (p1−p2) under mechanical equilibrium with equal pressure

p1 = p2 = p

p relaxation (Cont.)

We find easily

αkρk= αk0ρk0, ρ = ρ0, ~u = ~u0, E = E0, e = e0

∂t(αE)_k =∂t(αρe)_k= −pI∂tαk, k = 1, 2

p relaxation (Cont.)

We find easily

αkρk= αk0ρk0, ρ = ρ0, ~u = ~u0, E = E0, e = e0

∂t(αE)_k =∂t(αρe)_k= −pI∂tαk, k = 1, 2

Integrating latter equation & using αkρk = αk0ρk0 leads to ek(pk, ρk) − ek0+ ¯pI

 1 ρk

− 1 ρk0



= 0

This gives condition for ρk in p, k = 1, 2, if assume e.g.,

¯

pI = (p⁰_I + p)/2, & imposemechanical equilibrium in EOS

p relaxation (Cont.)

Combining that with saturation condition for volume fraction α₁ρ₁

ρ1(p) + α₂ρ₂ ρ2(p) = 1

leads to algebraic equation (quadratic one with SG EOS) for relaxed pressure p

With that, ρk, αk can be determined & state vector q is updated from current time to next

pT relaxation

Now assume frozen thermo-chemical relaxation ν = 0, look for solution of ODEs in limits µ & θ → ∞

∂t(α₁ρ₁) = 0

∂t(α2ρ2) = 0

∂t(ρ~u) = 0

∂t(α1E1) =µpI(p2−p1) + θTI(T2−T1)

∂t(α2E2) =µpI(p1−p2) + θTI(T1−T2)

∂tα₁ =µ (p₁−p₂)

under mechanical-thermal equilibrium conditons p1 = p2 = p

T1 = T2 = T

pT relaxation (Cont.)

As before, for k = 1, 2, states remain in equilibrium are αkρk = α_k0ρ_k0, ρ = ρ₀, ~u = ~u₀, E = E₀, e = e₀ Lead to equilibrium in mass fraction Yk = αkρk/ρ = Yk0

pT relaxation (Cont.)

As before, for k = 1, 2, states remain in equilibrium are αkρk = α_k0ρ_k0, ρ = ρ₀, ~u = ~u₀, E = E₀, e = e₀ Lead to equilibrium in mass fraction Yk = αkρk/ρ = Yk0

Impose mechanical-thermal equilibrium to 1. Saturation condition

α1ρ1

ρ1(p,T) + α2ρ2

ρ2(p,T) = 1

or Y 1

ρ1(p,T) + Y2

ρ2(p,T) = 1 ρ

pT relaxation (Cont.)

Impose mechanical-thermal equilibrium to 1. Saturation condition

Y 1

ρ1(p,T) + Y2

ρ2(p,T) = 1 ρ 2. Equilibrium of internal energy

Y1 e1(p,T) + Y2 e2(p,T) = e Give 2 algebraic equations for2 unknowns p & T

pT relaxation (Cont.)

Impose mechanical-thermal equilibrium to 1. Saturation condition

Y 1

ρ1(p,T) + Y2

ρ2(p,T) = 1 ρ 2. Equilibrium of internal energy

Y1 e1(p,T) + Y2 e2(p,T) = e Give 2 algebraic equations for2 unknowns p & T

For SG EOS, it reduces to single quadraticequation for p&

explicit computation of T: 1

ρT = Y1

(γ1−1)Cv1

p+ p^∞1

+ Y2

(γ2−1)Cv2

p+ p^∞2

pT g relaxation

Look for solution of ODEs in limits µ, θ, & ν → ∞

∂t(α₁ρ₁) = ν (g₂−g₁)

∂t(α2ρ2) = ν (g1−g2)

∂t(ρ~u) = 0

∂t(α1E1) = µpI(p2−p1) + θTI(T2 −T1) + ν (g2 −g1)

∂t(α2E2) = µpI(p1−p2) + θTI(T1 −T2) + ν (g1 −g2)

∂tα1 =µ (p1−p2) + νvI(g2−g1)

under mechanical-thermal-chemical equilibrium conditons p1 = p2 = p

T1 = T2 = T g1 = g2

pT g relaxation (Cont.)

In this case, states remain in equilibrium are

ρ = ρ0, ρ~u = ρ0~u0, E = E0, e = e0

but αkρk6= α_k0ρ_k0 & Yk 6= Y_k0, k = 1, 2

Impose mechanical-thermal-chemical equilibrium to 1. Saturation condition for temperature

G(p,T) = 0 2. Saturation condition for volume fraction

Y1

ρ1(p,T) + Y2

ρ2(p,T) = 1 ρ 3. Equilibrium of internal energy

Y1 e1(p,T) +Y2 e2(p,T) = e

pT g relaxation (Cont.)

From saturation condition for temperature G(p,T) = 0 we get T in terms of p, while from

Y1

ρ1(p,T)+ Y2

ρ2(p,T) = 1 ρ

&

Y1 e1(p,T) +Y2 e2(p,T) = e we obtain algebraic equation for p

Y1 = 1/ρ2(p) − 1/ρ

1/ρ2(p) − 1/ρ1(p) = e − e2(p) e1(p) − e2(p) which is solved by iterative method

pT g relaxation (Cont.)

With that, T can be solved from either condition for volume fraction or equilibrium of internal energy (quadratic equation for SG EOS), yielding ρk & αk update

Expansion wave problem: Cavitation test

Liquid-vapor mixture (αvapor = 1/5) with initial states pliquid= pvapor = 1bar

Tliquid = Tvapor = 354.7284K< T^sat

ρvapor = 0.63kg/m³> ρ^sat_vapor, ρliquid= 1150kg/m³> ρ^sat_liquid g^sat> gvapor > gliquid

← −~u ~u →

← Membrane

Cavitation test: ~u = 2m/s

Snap shot of computed solution at time t = 3.2ms

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10x 10⁴

p relaxation p−pT p−pT−pTg p−pTg

x(m)

p(Pa)

0 0.2 0.4 0.6 0.8 1

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

p relaxation p−pT p−pT−pTg p−pTg

x(m)

u(m/s)

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35

p relaxation p−pT p−pT−pTg p−pTg

x(m) αv

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4x 10⁻⁴

p relaxation p−pT p−pT−pTg p−pTg

x(m) Yv

Cavitation test: ~u = 500m/s

Snap shot of computed solution at time t = 0.58ms

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10x 10⁴

p relaxation p−pT p−pT−pTg p−pTg

x(m)

p(Pa)

0 0.2 0.4 0.6 0.8 1

−500 0 500

p relaxation p−pT p−pT−pTg p−pTg

x(m)

u(m/s)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

p relaxation p−pT p−pT−pTg p−pTg

x(m) αv

0 0.2 0.4 0.6 0.8 1

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

p relaxation p−pT p−pT−pTg p−pTg

x(m) Yv

Dodecane liquid-vapor shock tube problem

Snap shot of computed solution at time t = 473µs

0 0.2 0.4 0.6 0.8 1

10⁰ 10¹ 10² 10³

p relaxation p−pT−pTg p−pTg

ρ(kg/m3 )

0 0.2 0.4 0.6 0.8 1

0 50 100 150 200 250 300 350

p relaxation p−pT−pTg p−pTg

u(m/s)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

p relaxation p−pT−pTg p−pTg

αv

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

p relaxation p−pT−pTg p−pTg

Yv

0 0.2 0.4 0.6 0.8 1

10⁰ 10¹ 10² 10³

p relaxation p−pT−pTg p−pTg

p(bar)

0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5

2x 10⁶

p relaxation p−pT−pTg p−pTg

(Tv−Tl)(K)

Vapor-bubble compression: pT g relaxation results

Vapor mass fraction Mixture pressure

t=0.4ms

t=0.8ms

t=1.2ms

Future perspective

Low Mach solver for weakly compressible flow

Model with granular pressure for gas-liquid-soild flow Robust & stable stiff solver for real materials

Thank you