• 沒有找到結果。

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

N/A
N/A
Protected

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

Copied!
70
0
0

(1)

### 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/m3, 500kg/m3, 0, 108Pa, 108 Vapor phase: Right-hand side (0.75m < x ≤ 1m)

v, ρl, u, p, αv)R = 2kg/m3, 500kg/m3, 0, 105Pa, 1 − 108 ,

Liquid Vapor

← Membrane

(6)

### Dodecane 2-phase problem: Phase diagram

10−4 10−3 10−2 10−1 100 101 102 103

101 102 103 104 105 106 107 108

Liquid Vapor

2-phase mixture

Saturation curve →

v

p

(7)

### Dodecane 2-phase problem: Phase diagram

Wave path in p-v phase diagram

10−4 10−3 10−2 10−1 100 101 102 103

101 102 103 104 105 106 107 108

Liquid Vapor

2-phase mixture

Saturation curve →

v

p

Isentrope

Hugoniot locus

(8)

### Dodecane 2-phase problem: Sample solution

0 0.2 0.4 0.6 0.8 1

100 101 102 103

t = 0 t=473µs

Density (kg/m3)

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

104 105 106 107 108

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

(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 = 102) for water with pliquid = pvapor = 1bar

Tliquid= Tvapor = 354.7284K< Tsat

ρvapor = 0.63kg/m3> ρsatvapor, ρliquid= 1150kg/m3> ρsatliquid gsat> 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

102 103 104

t = 0 t=3.2ms

Density (kg/m3)

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

104 105

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−4

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 104

t = 0 t=3.2ms

colorblue gvgl(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−3 105

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−4 10−3 10−2 10−1 100 101 102

103 104 105 106 107

Liquid

Vapor 2-phase mixture

Saturation curve →

v

p

(15)

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

0 0.2 0.4 0.6 0.8 1

10−1 100 101 102 103 104

t = 0 t=0.58ms

Density (kg/m3)

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

103 104 105

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 104

t = 0 t=0.58ms

colorblue gvgl(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)ek− γkpk− (γk− 1)ρkηk

Temperature

Tk(pk, ρk)= pk+ pk

k− 1)Cvkρk

Entropy

sk(pk, Tk)= Cvklog Tkγk

(pk+ pk)γk1 + η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) 109 0

η (J/kg) −11.6 × 103 2030 × 103

η (J/(kg · K)) 0 −23.4 × 103

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

Fluid Dodecane

Parameters/Phase Liquid Vapor

γ 2.35 1.025

p (Pa) 4 × 108 0

η (J/kg) −775.269 × 103 −237.547 × 103

η (J/(kg · K)) 0 −24.4 × 103

Cv (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+p1)−log(p+p2) = 0 A= Cp1− Cp2+ η2− η1

Cp2− Cv2

, B= η1 − η2 Cp2− Cv2

C = Cp2− Cp1

Cp2− Cv2, D= Cp1− Cv1

Cp2− Cv2

(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+p1)−log(p+p2) = 0 A= Cp1− Cp2+ η2− η1

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 (s2− s1): 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−4 10−3 10−2

water dodecane Liquid volume v1(m3/kg)

300 350 400 450 500

10−2 10−1 100 101 102 103

water dodecane Vapor volume v2(m3/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(αρ)1+ ∇ · (αρ~u)1 = 0

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

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

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

t(αE)1+ ∇ · (αE~u + αp~u)1 =pI~uI · ∇α1+ µpI(p2− p1) +λ~uI · (~u2− ~u1)

t(αE)2+ ∇ · (αE~u + αp~u)2 = −pI~uI · ∇α1− µpI(p2− p1) −λ~uI · (~u2 − ~u1)

tα1+~uI · ∇α1 =µ(p1− p2) (α1+ α2 = 1) αk: volume fraction, ρk: density, ~uk: velocity pkk, 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

, SI(Interfacial area)

(24)

### 7 -equation model: With phase transition

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

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

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

t(αρ~u)1+ ∇ · (αρ~u ⊗ ~u)1+ ∇(αp)1 = pI∇α1+ λ (~u2− ~u1) +~uI

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

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

µpI(p2− p1) + λ~uI · (~u2− ~u1) +Q+(eI + ~uI · ~uI/2) ˙m

t(αE)2+ ∇ · (αE~u + αp~u)2 = −pI~uI · ∇α1

µpI(p2− p1) − λ~uI · (~u2− ~u1) −Q−(eI + ~uI · ~uI/2) ˙m

tα1+ ~uI · ∇α1 = µ (p1− p2) + 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 − α1) max(p − pv, 0) tρ1U2 /2

˙

m= Cliqα1ρ1min(p − pv, 0) ρvtρ1U2 /2 Kunz et al. (2000)

˙

m+ = Cprodα12(1 − α1) ρ1t

, m˙= Cliqα1ρvmin(p − pv, 0) ρ1tρ1U2/2

(26)

### Mass transfer modelling

Singhal et al. (2002)

˙

m+= Cprod√ κ σ ρ1ρv

 2 3

max(p − pv, 0) ρ1

1/2

˙

m= Cliq√ κ σ ρ1ρv

 2 3

min(p − pv, 0) ρ1

1/2

Senocak & Shyy (2004)

˙

m+= max(p − pv, 0)

1− ρc)(Vvn− V1n)2t, ˙m= ρ1min(p − pv, 0) ρv1− ρc)(Vvn− V1n)2t Hosangadi & Ahuja (JFE 2005)

˙

m+= Cprod

ρv

ρl(1 − α1)min(p − pv, 0) ρU2 /2

˙

m= Cliq

ρv

ρl

α1max(p − pv, 0) ρU2/2

(27)

### Phase transition model: 7-equation

We assume

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

˙

m = ν (g2− g1)

for mass transfer

θ ≥ 0 expresses rate towardsthermal equilibrium T1 → T2

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

metastable state Tliquid > Tsat

(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

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

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

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

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

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

s

Ks1∇ · ~u+ Q qI

+ m˙ ρI

s= α1

Ks1 + α2

Ks2

1

, qI = Ks1 α1

+Ks2 α2

   Γ1

α1

2

α2



ρI = Ks1 α1

+Ks2 α2

   c21 α1

+ c22 α2



, Ksι = ριc2ι

(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

ρc2p = α1

ρ1c21 + α2

ρ2c22 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

ρc2pT g = 1 ρc2p +T

"

α1ρ1

Cp1

 ds1

dp

2

+ α2ρ2

Cp2

 ds2

dp

2#

(32)

### Equilibrium speed of sound: Comparison

Sound speeds followsubcharacteristic condition cpT g≤cp

Sound speed limits follow

αlimk1cp = ck, lim

αk1cpT g 6= ck

0 0.2 0.4 0.6 0.8 1

10−1 100 101 102 103 104

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

Solve system of ordinary differential equations

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

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

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

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

t1E1) + ∇ · (α1E1~u + α1p1~u) + B (q, ∇q) = µpI(p2− p1) + Q + eI

t2E2) + ∇ · (α2E2~u + α2p2~u) − B (q, ∇q) = µpI(p1− p2) − Q − eI

tα1+ ~u · ∇α1 =µ (p1− p2) + Q qI

+ m˙ ρI

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

(38)

### Phase transition model: 6-equation

µ, θ, ν → ∞: instantaneous exchanges (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: θ (T2− T1) θexpresses rate towards thermal equilibrium T1→ T2,

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

isnonzero only at 2-phase mixture& metastable state Tliquid> Tsat

(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 + (α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, Q, −Q, Q/qI]T

ψν = [ ˙m, − ˙m, 0, eIm, −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 cpT g ≤ cpT ≤ cp ≤ cf

Limit of sound speed

αlimk1cf = lim

αk1cp = lim

αk1cpT = ck, lim

αk1cpT g 6= ck

0 0.2 0.4 0.6 0.8 1

10−1 100 101 102 103 104

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 (µ & θ → ∞)

Solve system of ordinary differential equations

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 qI, ρI, & eI

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 qI, ρI, & eI

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 µ → ∞

t1ρ1) = 0

t2ρ2) = 0

t(ρ~u) = 0

t1E1) = µpI(p2− p1)

t2E2) = µpI(p1− p2)

tα1 =µ (p1− p2)

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

p1 = p2 = p

(46)

### Stiff mechanical relaxation step (Cont.)

We find easily

t1ρ1) = 0 =⇒ α1ρ1 = α01ρ01

t2ρ2) = 0 =⇒ α2ρ2 = α02ρ02

t(ρ~u) = 0 =⇒ ρ~u = ρ0~u0

t1E1) = µpI(p2− p1) =⇒ ∂t(αρe)1 = −pItα1

t2E2) = µpI(p1− p2) =⇒ ∂t(αρe)2 = −pItα2

Integrating latter two equations with respect to time Z

t(αρe)k dt = − Z

pItαk dt

=⇒ αkρkek− α0kρ0ke0k = −p¯I αk− α0k

 or

=⇒ ek− e0k = −p¯I 1/ρk− 1/ρ0k

 (use αkρk = α0kρ0k) Take p¯I = (p0I + 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

100 101 102 103

t = 0 t=473µs

Density (kg/m3)

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

104 105 106 107 108

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

(50)

### Dodecane 2-phase problem: Phase diagram

Wave path after p-relaxationin p-v phase diagram

10−4 10−3 10−2 10−1 100 101 102 103

101 102 103 104 105 106 107 108

Liquid Vapor

2-phase mixture

Saturation curve →

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−4 10−3 10−2 10−1 100 101 102 103

101 102 103 104 105 106 107 108

Liquid Vapor

2-phase mixture

Saturation curve →

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

103.02 103.03 103.04 103.05

t = 0 t=3.2ms

Density (kg/m3)

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

103 104 105

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

Vapor volume fraction

(53)

### Expansion wave problem: Phase diagram

Wave path comparison between solutions after p- &

pT g-relaxation in p-v phase diagram

10−3 103

104 105 106 107

Liquid 2-phase mixture

v

p

(54)

### Stiff thermal relaxation step

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

t1ρ1) = 0

t2ρ2) = 0

t(ρ~u) = 0

t1E1) =µpI(p2− p1) + θ (T2− T1)

t2E2) =µpI(p1− p2) + θ (T1− T2)

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

t1ρ1) = 0 =⇒ α1ρ1 = α01ρ01

t2ρ2) = 0 =⇒ α2ρ2 = α02ρ02

t(ρ~u) = 0 =⇒ ρ~u = ρ0~u0

tkEk) = θ qI

(T2− T1) =⇒ ∂t(αρe)k= qItαk

Integrating latter two equations with respect to time Z

t(αρe)k dt = Z

qItαk dt

=⇒ αkρkek− α0kρ0ke0k = −q¯I αk− α0k

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

T2 e2, α02ρ02/(1 − α1) − T1 e1, α01ρ011 = 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 ρ0 2. Equilibrium of internal energy

Y1 e1(p,T) + Y2 e2(p,T) = e0 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+ p1

+ Y2

2− 1)Cv2

p+ p2

(57)

### Stiff thermo-chemical relaxation step

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

t1ρ1) =ν (g2 − g1)

t2ρ2) =ν (g1 − g2)

t(ρ~u) = 0

t1E1) =µpI(p2− p1) + θ (T2− T1) + ν (g2− g1)

t2E2) =µpI(p1− p2) + θ (T1− T2) + ν (g1− g2)

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

qI (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

ρ = ρ0, ρ~u = ρ0~u0, E = E0, e = e0 but αkρk6= α0kρ0k & Yk6= Yk0, 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 ρ0 3. Equilibrium of internal energy

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

(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 ρ0

&

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

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

1/ρ2(p) − 1/ρ1(p) = e0− e2(p) e1(p) − e2(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) = e0 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/m3)

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

104 105 106 107 108

p relax pT relax pTg relax

Pressure(bar)

0 0.2 0.4 0.6 0.8 1

−6

−4

−2 0 2 4 6x 108

p relax pT relax pTg relax

gvgl(103K)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

p relax pT relax pT relax

Vapor volume fraction

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

p relax pT relax pTg relax

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 104

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

p[Pa]

0 0.2 0.4 0.6 0.8 1

−500 0 500

p relaxation p−pT 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 p−pT−pTg p−pTg

α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

Yv

(63)

### High-pressure fuel injector

Inject fluid: Liquid dodecane containing small amountαvapor

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

Ambient fluid: Vapor dodecanecontaining small amount αliquid

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

liquid →

vapor

(64)

v,l

−4

### ): 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 = 104 (left) vs. 102 (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)

### Thank you

Key words: Compressible two-phase flow, Five-equation model, Interface sharpening, THINC reconstruction, Mie-Gr¨ uneisen equation of state, Semi-discrete wave propagation method..

Key words: Compressible two-phase flow, Five-equation model, Interface sharpening, THINC reconstruction, Mie-Gr¨ uneisen equation of state, Semi-discrete wave propagation method..

Standard interface capturing results for toy problem, observing poor interface resolution.. Passive advection

Comments on problems with complex equation of state Even for single phase flow case, standard method will fail to attain pressure equilibrium, say, for example, van der Waals gas..

• Many statistical procedures are based on sta- tistical models which specify under which conditions the data are generated.... – Consider a new model of automobile which is

The main tool in our reconstruction method is the complex geometri- cal optics (CGO) solutions with polynomial-type phase functions for the Helmholtz equation.. This type of

FIGURE 5. Item fit p-values based on equivalence classes when the 2LC model is fit to mixed-number data... Item fit plots when the 2LC model is fitted to the mixed-number

Then, based on these systematically generated smoothing functions, a unified neural network model is pro- posed for solving absolute value equationB. The issues regarding