Relaxation methods for compressible multiphase flow

Relaxation methods for

compressible multiphase flow

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

Joint work: Marica Pelanti (Paris Tech)

Outline

1. Model problems

Shock diffraction in 2-phase mixture Cavitating Richtmyer-Meshkov instability High pressure fuel injector

2. Homogeneous relaxation model (HRM)

Numerical modelling of wave dynamics in multiphase mixtures of compressible flow

3. Numerical scheme

Finite volume method Stiff relaxation solver 4. Numerical examples

Shock wave diffraction down a step: Van Dyke

Shock wave diffraction in 2-phase mixture

Shock wave induced by high pressure liquid injection Inject liquid pressure p = 10⁸ Pa

Ambientgas pressure p = 10⁵ Pa

Liquid & gas inthermal equilibrium with T = 640 K α-dependent homogeneous fluid (liquid & gas)mixture

liquid →

gas

Shock diffraction in 2-phase mixture: α = 10

Mixture density Mixture pressure

Shock diffraction in 2-phase mixture: α = 10

Mixture density Mixture pressure

Shock diffraction in 2-phase mixture (Cont.)

α = 10⁻⁴ (left) & α = 10⁻² (right), vapor temperature, liquid temperature, & volume fraction (top to bottom),

Cavitating Richtmyer-Meshkov problem

Gas volume fraction Mixture pressure

t=0ms

t=2ms

t=3.1ms

t=6.4ms

t=8.6ms

High-pressure fuel injector

With thermo-chemical relaxation No thermo-chemical relaxation

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(α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) B(q, ∇q) is non-conservative product (q: state vector)

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

Homogeneous relaxation model (Cont.)

1. µ (p1−p2): Volume transfer via pressure relaxation µexpresses rate toward mechanical equilibrium p1 → p₂,

& isnonzero in all flow regimes of interest

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

isnonzero only at 2-phase mixture

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

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

If µ, θ, ν → ∞: stiff (instantaneous) exchanges

Homogeneous relaxation model (Cont.)

HRM model in compact form

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

q= [α1, α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, θ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 (stiff or non-stiff limit)

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

Homogeneous relaxation model (Cont.)

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

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

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

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∆ξ1

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) x2 = x2(ξ1, ξ2) 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−1 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

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

pT relaxation

Now 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) =µpI(p2−p1) + θTI(T2−T1)

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

∂tα1 =µ (p1−p2)

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, ρ = ρ0, ~u = ~u0, E = E0, e = e0

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, ρ = ρ0, ~u = ~u0, E = E0, e = e0

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(α1ρ1) = ν (g2−g1)

∂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= 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 ρ 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)

High-speed underwater projectile

With thermo-chemical relaxation No thermo-chemical relaxation

Thank you

Constitutive law

Stiffened gas equation of state (SG EOS) with Pressure

pk(ek, ρk)= (γk−1)ek−γkp∞k−(γk−1)ρkη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 (s2−s1): 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)