Recent advances in numerical methods for compressible two-phase ﬂow with heat & mass transfers

(1)

Recent advances in

numerical methods for compressible two-phase flow with heat & mass transfers

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

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

(2)

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. Mathematical model with & without heat & mass transfer 4. Stiff relaxation solver

(3)

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. Mathematical model with & without heat & mass transfer 4. Stiff relaxation solver

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

Solver preserves total energy conservation & employ convex pressure law

(4)

Phase transition with non-convex EOS

Sample wave path for phase transition problem 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

(5)

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

(6)

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

(7)

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

v

p

Isentrope

Hugoniot locus

(8)

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

t = 0 t=473µs

Vapor volume fraction

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

(9)

Dodecane 2-phase problem: Sample solution

Allphysical quantities arediscon- tinuous across phase boundary

(10)

Expansion wave problem: Cavitation test

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

Liquid-vapor mixture(αvapor = 10⁻²) for water with pliquid = pvapor = 1bar

Tliquid= Tvapor = 354.7284K< T^sat

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

Outgoing velocityu = 2m/s

← −~u ~u →

← Membrane

(11)

Expansion wave problem: Sample solution

Cavitation pocket formation &

mass transfer

(12)

Expansion wave problem: Sample solution

0 0.2 0.4 0.6 0.8 1

10² 10³ 10⁴

t = 0 t=3.2ms

Density (kg/m³)

0 0.2 0.4 0.6 0.8 1

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

t = 0 t=3.2ms

Velocity (m/s)

0 0.2 0.4 0.6 0.8 1

10⁴ 10⁵

t = 0 t=3.2ms

Pressure(bar)

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⁻⁴

t = 0 t=3.2ms

Vapor mass fraction

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10 12x 10⁴

t = 0 t=3.2ms

colorblue gv−gl(J/kg)

Equilibrium Gibbs free energy inside cavitation pocket

(13)

Expansion wave problem: Phase diagram

Solution remains in 2-phase mixture; phase separation has not reached

10⁻³ 10⁵

Liquid 2-phase mixture

v

p

(14)

Expansion wave ~u = 500m/s: Phase diagram

With faster ~u = 500m/s, phase separation becomes more evident

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

10³ 10⁴ 10⁵ 10⁶ 10⁷

Liquid

Vapor 2-phase mixture

v

p

(15)

Expansion wave ~u = 500m/s: Sample solution

0 0.2 0.4 0.6 0.8 1

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

t = 0 t=0.58ms

Density (kg/m³)

0 0.2 0.4 0.6 0.8 1

−500 0 500

t = 0 t=0.58ms

Velocity (m/s)

0 0.2 0.4 0.6 0.8 1

10³ 10⁴ 10⁵

t = 0 t=0.58ms

Pressure(bar)

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

t = 0 t=0.58ms

Vapor mass fraction

0 0.2 0.4 0.6 0.8 1

−2 0 2 4 6 8 10 12x 10⁴

t = 0 t=0.58ms

colorblue gv−gl(J/kg)

Equilibrium Gibbs free energy inside cavitation pocket

(16)

Constitutive law: Metastable fluid

Stiffened gas equation of state (SG EOS) with Pressure

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

Temperature

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

(γk− 1)C^vkρk

Entropy

sk(pk, Tk)= Cvklog Tkγk

(p_k+ p∞k)^γ^k⁻¹ + η_k^′ Helmholtz free energy ak = ek− Tksk

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

(17)

Metastable fluid: 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³

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

(18)

Metastable fluid: Saturation curves

Assume two phases in chemicalequilibrium with 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− C^p2+ η^′₂− η^′1

Cp2− Cv2

, B= η1 − η² Cp2− Cv2

C = Cp2− Cp1

Cp2− C^v2, D= Cp1− Cv1

Cp2− C^v2

(19)

Metastable fluid: Saturation curves

Assume two phases in chemicalequilibrium with 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− C^p2+ η^′₂− η^′1

Cp2− Cv2

, B= η1 − η² Cp2− Cv2

C = Cp2− Cp1

Cp2− C^v2, D= Cp1− Cv1

Cp2− C^v2 or, from dg1 = dg2, we get Clausius-Clapeyron equation

dp(T )

dT = Lh

T (v2− v1) L_h = T (s2− s¹): latent heat of vaporization

(20)

Metastable fluid: 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)

(21)

Mathematical Models

Phase tranition models for compressible 2-phase flow include 1. 7-equation model (Baer-Nunziato type)

Zein, Hantke, Warnecke (JCP 2010) 2. Reduced 5-equation model (Kapila type)

Saurel, Petitpas, Berry (JFM 2008) 3. Homogeneous 6-equation model

Zein et al. , Saurel et al. , Pelanti & Shyue (JCP 2014) 4. Homogeneous equilibrium model

Dumbser, Iben, & Munz (CAF 2013), Hantke, Dreyer, &

Warnecke (QAM 2013) 5. Navier-Stokes-Korteweg model

Prof. Kr¨oner’s talk tomorrow

(22)

7 -equation model: Without phase transition

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

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

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

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

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

∂t(αE)₁+ ∇ · (αE~u + αp~u)1 =pI~uI · ∇α1+ µpI(p2− p¹) +λ~uI · (~u²− ~u¹)

∂t(αE)₂+ ∇ · (αE~u + αp~u)2 = −pI~uI · ∇α¹− µpI(p2− p¹) −λ~uI · (~u² − ~u¹)

∂_tα1+~u_I · ∇α¹ =µ(p1− p²) (α1+ α2 = 1) αk: volume fraction, ρk: density, ~uk: velocity pk(ρk, ek): pressure, ek: specific internal energy

Ek = ρkek+ ρk~uk· ~uk/2: specific total energy, k = 1, 2

(23)

7 -equation model: Closure relations

pI & ~uI: interfacial pressure & velocity, e.g., 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

µ & λ: non-negative relaxation parameters that express rates pressure & velocity toward equilibrium, respectively

µ = SI

Z1+ Z2

, λ = SIZ1Z2

Z1+ Z2

, S_I(Interfacial area)

(24)

7 -equation model: With phase transition

7-equation model with heat& mass transfers (Zein et al. ):

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

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

∂t(αρ~u)₁+ ∇ · (αρ~u ⊗ ~u)1+ ∇(αp)¹ = pI∇α¹+ λ (~u2− ~u¹) +~uIm˙

∂t(αρ~u)₂+ ∇ · (αρ~u ⊗ ~u)2+ ∇(αp)² = −p^I∇α1− λ (~u2− ~u1)−~uIm˙

∂t(αE)₁+ ∇ · (αE~u + αp~u)1 = pI~uI · ∇α¹+

µpI(p2− p¹) + λ~uI · (~u²− ~u¹) +Q+(eI + ~uI · ~u^I/2) ˙m

∂t(αE)₂+ ∇ · (αE~u + αp~u)2 = −p^I~uI · ∇α¹−

µp_I(p2− p¹) − λ~uI · (~u²− ~u¹) −Q−(e_I + ~u_I · ~uI/2) ˙m

∂tα1+ ~uI · ∇α¹ = µ (p1− p²) + Q qI

+ m˙ ρI

(25)

Mass transfer modelling

Typical apporach to mass transfer modelling assumes

˙

m = ˙m⁺+ ˙m⁻

Singhal et al. (1997) & Merkel et al. (1998)

˙

m⁺= Cprod(1 − α¹) max(p − p^v, 0) t^∞ρ1U∞² /2

˙

m⁻= Cliqα1ρ1min(p − p^v, 0) ρvt∞ρ1U∞² /2 Kunz et al. (2000)

˙

m⁺ = Cprodα₁²(1 − α¹) ρ1t∞

, m˙⁻= Cliqα1ρvmin(p − p^v, 0) ρ1t∞ρ1U∞²/2

(26)

Mass transfer modelling

Singhal et al. (2002)

˙

m⁺= Cprod√ κ σ ρ1ρv

2 3

max(p − p^v, 0) ρ1

1/2

˙

m⁻= Cliq√ κ σ ρ1ρv

2 3

min(p − p^v, 0) ρ1

1/2

Senocak & Shyy (2004)

˙

m⁺= max(p − p^v, 0)

(ρ1− ρc)(Vvn− V1n)²t^∞, ˙m⁻= ρ1min(p − p^v, 0) ρv(ρ1− ρc)(Vvn− V1n)²t^∞ Hosangadi & Ahuja (JFE 2005)

˙

m⁺= Cprod

ρv

ρl(1 − α¹)min(p − p^v, 0) ρ∞U∞² /2

˙

m⁻= Cliq

ρv

ρl

α1max(p − p^v, 0) ρ^∞U∞²/2

(27)

Phase transition model: 7-equation

We assume

Q = θ (T2− T1) for heat transfer &

˙

m = ν (g2− g¹)

for mass transfer

θ ≥ 0 expresses rate towardsthermal equilibrium T1 → T2

ν ≥ 0expresses rate towards diffusive equilibrium g1 → g², & isnonzero only at2-phase mixture &

metastable state T_liquid > T_sat

(28)

7 -equation model: Numerical approximation

Write 7-equation model in compact form

∂tq + ∇ · f(q) + w (q, ∇q) = ψ^µ(q) + ψλ(q) + ψθ(q) + ψν(q) Solve by fractional-step method

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) + ψ_ν(q) in various flow regimes underrelaxation limits

(29)

Reduced 5-equation model: With phase transition

Saurel et al. considered 7-equation model in asymptotic limits λ & µ → ∞, i.e., flow towards mechanical equilibrium:

~u1 = ~u2 = ~u&p1 = p2 = p, i.e., reduced 5-equationmodel

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

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

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

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

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

K¯_s

K_s¹∇ · ~u+ Q qI

+ m˙ ρI

K¯s= α1

K_s¹ + α2

K_s²

⁻1

, qI = K_s¹ α1

+K_s² α2

Γ1

α1

+Γ2

α2

ρI = K_s¹ α1

+K_s² α2

c²₁ α1

+ c²₂ α2

, K_s^ι = ριc²_ι

(30)

Phase transition model: 5-equation

Mixture entropys = Y1s1+ Y2s2 admits nonnegative variation

∂t(ρs) + ∇ · (ρs~u) ≥ 0

Mixture pressurep determined from total internal energy ρe = α1ρ1e1(p, ρ1) + α2ρ2e2(p, ρ2)

Model is hyperbolicwith non-monotonic sound speed cp: 1

ρc²_p = α1

ρ1c²₁ + α2

ρ2c²₂ Limit interface model, i.e., as θ & ν → ∞

(thermo-chemical relaxation), is homogeneous equilibrium model

(31)

Homogeneous equilibrium model

Homogeneous equilibrium model (HEM) 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

System is closed by

p1 = p2 = p, T1 = T2 = T, & g1 = g2 = g Speed of sound cpT g satisfies

1

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

"

α1ρ1

C_p1

ds1

dp

2

+ α2ρ2

C_p2

ds2

dp

2#

(32)

Equilibrium speed of sound: Comparison

Sound speeds followsubcharacteristic condition cpT g≤cp

Sound speed limits follow

αlimk→1cp = ck, lim

αk→1cpT g 6= c^k

0 0.2 0.4 0.6 0.8 1

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

p−relax pTg−relax

αwater

cp&cpTg

(33)

5 -equation model: Numerical approximation

Write 5-equation model in compact form

∂tq + ∇ · f(q) + w (q, ∇q) = ψ^θ(q) + ψν(q) Solve by fractional-step method

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

∂_tq = ψ_θ(q) + ψ_ν(q) in various flow regimes underrelaxation limits

(34)

HEM: Numerical approximation

Write HEM in compact form

∂tq + ∇ · f(q) = 0

Compute solution numerically, e.g., Godunov-type method, requires Riemann solverfor elementary waves to fulfil

1. Jump conditions across discontinuities 2. Kinetic condition

3. Entropy condition

(35)

Numerical approximation: summary

1. Solver based on7-equation model is viableone for wide variety of problems, but is expensive to use

2. Solver based on reduced 5-equation model isrobust one for sample problems, but is difficult to achieve admissible solutions under extreme flow conditions

3. Solver based onHEM is mathematically attracive one

(36)

Numerical approximation: summary

1. Solver based on7-equation model is viableone for wide variety of problems, but is expensive to use

2. Solver based on reduced 5-equation model isrobust one for sample problems, but is difficult to achieve admissible solutions under extreme flow conditions

3. Solver based onHEM is mathematically attracive one Numerically advantageous to use 6-equationmodel as opposed to 5-equation model (Saurel et al. , Pelanti & Shyue)

(37)

6 -equation model: With phase transition

6-equation single-velocity 2-phase model with stiff mechanical, thermal, & chemical relaxationsreads

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

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

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

∂t(α1E1) + ∇ · (α¹E1~u + α1p1~u) + B (q, ∇q) = µpI(p2− p¹) + Q + e^Im˙

∂_t(α2E2) + ∇ · (α²E2~u + α2p2~u) − B (q, ∇q) = µp_I(p1− p²) − Q − eIm˙

∂_tα1+ ~u · ∇α¹ =µ (p1− p²) + Q qI

+ m˙ ρI

B (q, ∇q) is non-conservative product (q: state vector) B = ~u · [Y1∇ (α2p2) − Y²∇ (α1p1)]

(38)

Phase transition model: 6-equation

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

µexpresses rate toward mechanical equilibrium p1 → p2,

& isnonzero in all flow regimes of interest

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

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

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

(39)

Phase transition model: 6-equation

6-equation model in compact form

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

q= [α1ρ1, α2ρ2, ρ~u, α1E1, α2E2, α1]^T

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

ψµ= [0, 0, 0, µpI(p2− p¹) , µpI(p1− p²) , µ (p1− p²)]^T ψθ = [0, 0, 0, Q, −Q, Q/q^I]^T

ψν = [ ˙m, − ˙m, 0, e^Im, −e˙ ^Im, ˙˙ m/ρI]^T

(40)

Phase transition model: 6-equation

Flow hierarchy in 6-equation model: H. Lund (SIAP 2012)

6-eqns µ → ∞ θν → ∞ µθν → ∞

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(41)

Phase transition model: 6-equation

Stiff limits as µ → ∞, µθ → ∞, & µθν → ∞ sequentially

6-eqns µ → ∞ θν → ∞ µθν → ∞

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(42)

Equilibrium speed of sound: Comparison

Sound speeds follow subcharacteristic condition c_{pT g} ≤ cpT ≤ cp ≤ cf

Limit of sound speed

αlimk→1cf = lim

αk→1cp = lim

αk→1cpT = ck, lim

αk→1cpT g 6= c^k

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

(43)

6 -equation model: Numerical approximation

As before, we begin by solving non-stiff hyperbolic equations in step 1, & continue by applying 3 sub-steps as

2. Stiff mechanical relaxation step

Solve system of ordinary differential equations (µ → ∞)

∂tq = ψµ(q)

with initial solution from step 1 as µ → ∞ 3. Stiff thermal relaxation step (µ & θ → ∞)

∂_tq = ψ_µ(q) + ψ_θ(q)

4. Stiff thermo-chemical relaxation step (µ, θ, & ν → ∞) Solve system of ordinary differential equations

∂tq = ψµ(q) + ψθ(q) + ψν(q) Take solution from previous step as initial condition

(44)

6 -equation model: Stiff relaxation solvers

1. Algebraic-basedapproach

Imposeequilibrium conditions directly, without making explicit of interface states q_I, ρ_I, & e_I

Saurel et al. (JFM 2008), Zein et al. (JCP 2010), LeMartelot et al. (JFM 2013), Pelanti-Shyue (JCP 2014) 2. Differential-based approach

Imposedifferential of equilibrium conditions, require explicit of interface states q_I, ρ_I, & e_I

Saurel et al. (JFM 2008), Zein et al. (JCP 2010) 3. Optimization-based approach (for mass transfer only)

Helluy & Seguin (ESAIM: M2AN 2006), Faccanoni et al.(ESAIM: M2AN 2012)

(45)

Stiff mechanical relaxation step

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− p²)

∂tα1 =µ (p1− p²)

with initial condition q⁰ (solution after non-stiff hyperbolic step) & under mechanical equilibrium condition

p1 = p2 = p

(46)

Stiff mechanical relaxation step (Cont.)

We find easily

∂t(α1ρ1) = 0 =⇒ α1ρ1 = α⁰₁ρ⁰₁

∂_t(α2ρ2) = 0 =⇒ α2ρ2 = α⁰₂ρ⁰₂

∂_t(ρ~u) = 0 =⇒ ρ~u = ρ⁰~u⁰

∂t(α1E1) = µpI(p2− p1) =⇒ ∂t(αρe)₁ = −p^I∂tα1

∂t(α2E2) = µpI(p1− p2) =⇒ ∂t(αρe)₂ = −p^I∂tα2

Integrating latter two equations with respect to time Z

∂t(αρe)_k dt = − Z

pI∂tαk dt

=⇒ αkρkek− α⁰kρ⁰_ke⁰_k = −p¯I αk− α⁰k

or

=⇒ ek− e⁰k = −p¯I 1/ρk− 1/ρ⁰k

(use αkρk = α⁰_kρ⁰_k) Take p¯I = (p⁰_I + p)/2or p, for example

(47)

Stiff mechanical relaxation step (Cont.)

We find condition for ρ_k in p, k = 1, 2

Combining that with saturation condition for volume fraction α1+ α2 = α1ρ1

ρ1(p) + α2ρ2

ρ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

(48)

Stiff mechanical relaxation step (Cont.)

We find condition for ρ_k in p, k = 1, 2

Combining that with saturation condition for volume fraction α1+ α2 = α1ρ1

ρ1(p) + α2ρ2

ρ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

Relaxed solution depends strongly on initial condition from non-stiff hyperbolic step

(49)

Dodecane 2-phase Riemann problem: p relaxation

Mechanical-equilibrium solution at t = 473µs

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

−20 0 20 40 60 80 100 120 140 160

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

t = 0 t=473µs

(50)

Dodecane 2-phase problem: Phase diagram

Wave path after p-relaxationin 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

v

p

Isentrope

Hugoniot locus

(51)

Dodecane 2-phase problem: Phase diagram

Wave path comparison between solutions after p- &

pT g-relaxation 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

v

p

Isentrope

Hugoniot locus

(52)

Expansion wave problem: p relaxation

Mechanical-equilibrium solution at t = 3.2ms

0 0.2 0.4 0.6 0.8 1

10^3.02 10^3.03 10^3.04 10^3.05

t = 0 t=3.2ms

Density (kg/m³)

0 0.2 0.4 0.6 0.8 1

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

t = 0 t=3.2ms

Velocity (m/s)

0 0.2 0.4 0.6 0.8 1

10³ 10⁴ 10⁵

t = 0 t=3.2ms

Pressure(bar)

0 0.2 0.4 0.6 0.8 1

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

t = 0 t=3.2ms

(53)

Expansion wave problem: Phase diagram

Wave path comparison between solutions after p- &

pT g-relaxation in p-v phase diagram

10⁻³ 10³

10⁴ 10⁵ 10⁶ 10⁷

Liquid 2-phase mixture

v

p

(54)

Stiff thermal relaxation step

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

∂t(α1ρ1) = 0

∂t(α2ρ2) = 0

∂t(ρ~u) = 0

∂_t(α1E1) =µp_I(p2− p¹) + θ (T2− T¹)

∂t(α2E2) =µpI(p1− p²) + θ (T1− T²)

∂tα1 =µ (p1− p2) + θ qI

(T2− T1) under mechanical-thermal equilibrium conditons

p1 = p2 = p T1 = T2 = T

(55)

Stiff thermal relaxation step (Cont.)

We find easily

∂t(α1ρ1) = 0 =⇒ α1ρ1 = α⁰₁ρ⁰₁

∂_t(α2ρ2) = 0 =⇒ α2ρ2 = α⁰₂ρ⁰₂

∂t(ρ~u) = 0 =⇒ ρ~u = ρ⁰~u⁰

∂t(αkEk) = θ qI

(T2− T¹) =⇒ ∂t(αρe)_k= qI∂tαk

Integrating latter two equations with respect to time Z

∂t(αρe)_k dt = Z

qI∂tαk dt

=⇒ α_kρ_ke_k− α⁰kρ⁰_ke⁰_k = −q¯_I α_k− α⁰k

Take q¯I = (q_I⁰+ qI)/2 or qI, for example, & find algebraic equation for α1, by imposing

T2 e2, α⁰₂ρ⁰₂/(1 − α¹) − T1 e1, α⁰₁ρ⁰₁/α1 = 0

(56)

Stiff thermal relaxation step: Algebraic approach

Impose mechanical-thermal equilibrium directly 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

(57)

Stiff thermo-chemical relaxation step

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

∂_t(α1ρ1) =ν (g2 − g¹)

∂_t(α2ρ2) =ν (g1 − g²)

∂t(ρ~u) = 0

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

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

∂tα1 =µ (p1− p2) + θ

q_I (T2− T1) + ν

ρ_I (g2− g1) under mechanical-thermal-chemical equilibrium conditons

p1 = p2 = p T1 = T2 = T

g1 = g2

(58)

Stiff thermal-chemical relaxation step (Cont.)

In this case, states remain in equilibrium are

ρ = ρ⁰, ρ~u = ρ⁰~u⁰, E = E⁰, e = e⁰ but α_kρ_k6= α⁰kρ⁰_k & Y_k6= Yk⁰, 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⁰

(59)

Stiff thermal-chemical relaxation step (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/ρ¹(p) = e⁰− e²(p) e1(p) − e²(p) which is solved by iterative method

(60)

Stiff thermal-chemical relaxation step (Cont.)

Having known Yk & p, T can be solved from, e.g., Y1 e1(p,T) +Y2 e2(p,T) = e⁰ yielding updateρk & αk

Feasibility of solutions, i.e., positivity of physical quantities ρk, αk, p, & T , for example

Employ hybridmethod i.e., combination of above method with differential-based approach (not discuss here), when it becomes necessary

(61)

Dodecane 2-phase Riemann problem

Comparison p-,pT -& p-pT g-relaxation solution at t = 473µs

0 0.2 0.4 0.6 0.8 1

0 100 200 300 400 500

p relax pT relax pTg relax

Density (kg/m³)

0 0.2 0.4 0.6 0.8 1

−50 0 50 100 150 200 250 300 350

p relax pT relax pT relax

Velocity (m/s)

0 0.2 0.4 0.6 0.8 1

10⁴ 10⁵ 10⁶ 10⁷ 10⁸

Pressure(bar)

0 0.2 0.4 0.6 0.8 1

−6

−4

−2 0 2 4 6x 10⁸

gv−gl(10³K)

0 0.2 0.4 0.6 0.8 1

p relax pT relax pT relax

0 0.2 0.4 0.6 0.8 1

Vapor mass fraction

(62)

Expansion wave problem: ~u = 500m/s

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

p[Pa]

0 0.2 0.4 0.6 0.8 1

−500 0 500

u[m/s

0 0.2 0.4 0.6 0.8 1

α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

Yv

(63)

High-pressure fuel injector

Inject fluid: Liquid dodecane containing small amountαvapor

Pressure & temperature are in equilibrium with p = 10⁸ Pa & T = 640K

Ambient fluid: Vapor dodecanecontaining small amount αliquid

Pressure & temperature are in equilibrium with p = 10⁵ Pa & T = 1022K

liquid →

vapor

(64)

High-pressure fuel injector (α

_v,l

= 10

⁻⁴

): p-relax

Mixture density Mixture pressure

(65)

High-pressure fuel injector: p-relax

Vapor volume fraction Vapor mass fraction

(66)

Fuel injector: p-pT -pT g relaxation

Vapor mass fraction: αv,l = 10⁻⁴ (left) vs. 10⁻² (right)

(67)

High-pressure fuel injector

With thermo-chemical relaxation No thermo-chemical relaxation

(68)

High-speed underwater projectile

With thermo-chemical relaxation No thermo-chemical relaxation

(69)

References

M. Pelanti and K.-M. Shyue. A mixture-energy-consistent 6-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves. J. Comput.

Phys., 259:331–357, 2014.

R. Saurel, F. Petitpas, and R. Abgrall. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. J. Fluid. Mech., 607:313–350, 2008.

R. Saurel, F. Petitpas, and R. A. Berry. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. k J. Comput. Phys.,

228:1678–1712, 2009.

A. Zein, M. Hantke, and G. Warnecke. Modeling phase transition for compressible two-phase flows applied to metastable liquids. J. Comput. Phys., 229:2964–2998, 2010.

(70)

Recent advances in numerical methods for compressible two-phase ﬂow with heat & mass transfers