• 沒有找到結果。

# Relaxation methods for compressible multiphase ﬂow

N/A
N/A
Protected

Share "Relaxation methods for compressible multiphase ﬂow"

Copied!
48
0
0

(1)

### compressible multiphase flow

Keh-Ming Shyue

Institute of Applied Mathematical Sciences National Taiwan University

Joint work: Marica Pelanti (Paris Tech)

(2)

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

(3)

(4)

### Shock wave diffraction in 2-phase mixture

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

Ambientgas pressure p = 105 Pa

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

liquid →

gas

(5)

### Shock diffraction in 2-phase mixture: α = 10

−4

Mixture density Mixture pressure

(6)

### Shock diffraction in 2-phase mixture: α = 10

−2

Mixture density Mixture pressure

(7)

### Shock diffraction in 2-phase mixture (Cont.)

α = 10−4 (left) & α = 10−2 (right), vapor temperature, liquid temperature, & volume fraction (top to bottom),

(8)

### Cavitating Richtmyer-Meshkov problem

Gas volume fraction Mixture pressure

t=0ms

t=2ms

t=3.1ms

t=6.4ms

t=8.6ms

(9)

### High-pressure fuel injector

With thermo-chemical relaxation No thermo-chemical relaxation

(10)

### Homogeneous relaxation model (HRM)

Consider HRM for 2-phase flow of form

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

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

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

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

t2E2) + ∇ · (α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)]

(11)

### Homogeneous relaxation model (Cont.)

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

& isnonzero in all flow regimes of interest

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

isnonzero only at 2-phase mixture

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

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

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

(12)

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

(13)

### Homogeneous relaxation model (Cont.)

Flow hierarchy in HRM (stiff or non-stiff limit)

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(14)

### Homogeneous relaxation model (Cont.)

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

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(15)

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

(16)

### 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 →

(17)

### Hyperbolic step: Mapped grid (Cont.)

Then finite volume method on control volume C reads Qn+1 = Qn− ∆t

M(C)

Ns

X

j=1

hjj −∆tWM(C)

Qn :=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

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

W cell average of w in cell C

(18)

### Hyperbolic step: Mapped grid (Cont.)

Assume mapped (i.e., logically rectangular) grid in 2D, we get Qn+1ij = Qnij − ∆t

κij∆ξ1

Fi+1 1

2,j−Fi−1 1

2,j



∆t κij∆ξ2

Fi,j+2 1

2

−Fi,j−2 1

2

−∆tWij∆ξ1∆ξ2

i − 1 i − 1

i

i j

j

j + 1 j + 1

ij

Cij

ξ1

ξ2

mapping

∆ξ1

∆ξ2

logical grid physical grid

←−

x1 = x11, ξ2) x2 = x21, ξ2) x1

x2

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

(19)

### Mapped grid method: Wave propagation (Cont.)

Godunov-type in wave propagation form is Qn+1ij = Qnij− ∆t

κij∆ξ1

A+1∆Qi−1

2,j+ A1∆Qi+1

2,j

−

∆t κij∆ξ2



A+2∆Qi,j−1

2 + A2∆Qi,j+1

2



A+1∆Q

A1∆Q

fluctuationsA+1∆Qi−1 2,j, A

1∆Qi+1

2,j, A+2∆Qi,j−1

2, &

A2∆Qi,j+1

2 : Solve

one-dimensional Riemann problemsin direction normal to cell edges

Wij may be included in fluctuations

(20)

### Mapped grid mtehod: Wave propagation (Cont.)

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

Qn+1ij := Qn+1ij − ∆t κij∆ξ1

Fei+1 1

2,j− eFi−1 1

2,j

−

∆t κij∆ξ2

Fei,j+2 1

2

− eFi,j−2 1

2



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

Fei−1 1

2,j = 1 2

Nw

X

k=1

λ1,ki−12,j

1 − ∆t κi−1

2,j∆ξ1

λ1,ki−12,j

! Wfi−1,k1

2,j

κi−1

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

2,j is limited waves to avoid oscillations near discontinuities

(21)

### Mapped grid method: Wave propagation (Cont.)

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

A2A+1∆Q

A+2A1∆Q

A+2A+1∆Q

A2A1∆Q

(22)

### 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,ki−1

2,j

κip,j∆ξ1

, λ2,ki,j−1

2

κi,jp∆ξ2

 ≤ 1,

ip = i if λ1,ki−1

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

(23)

### 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+1∆Qi−1

2,j+ A1∆Qi+1

2,j+ A1∆Qij



− 1

κij∆ξ2



A+2∆Qi,j−1

2 + A2∆Qi,j+1

2 + A2∆Qij



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

Q = Qn+ ∆tL (Qn) Q∗∗ = 3

4Qn+1

4Q+1

4∆tL (Q) Qn+1 = 1

3Qn+2

3Q+2

3∆tL (Q∗∗)

(24)

### 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)

(25)

### p relaxation

Assume frozen thermal & thermo-chemical relaxation, i.e., θ = 0 & ν = 0, 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) under mechanical equilibrium with equal pressure

p1 = p2 = p

(26)

### p relaxation (Cont.)

We find easily

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

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

(27)

### p relaxation (Cont.)

We find easily

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

t(αE)k =∂t(αρe)k= −pItα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 = (p0I + p)/2, & imposemechanical equilibrium in EOS

(28)

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

(29)

### pT relaxation

Now 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) + θTI(T2−T1)

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

tα1 =µ (p1−p2)

under mechanical-thermal equilibrium conditons p1 = p2 = p

T1 = T2 = T

(30)

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

(31)

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

(32)

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

(33)

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

(34)

### pT g relaxation

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

t1ρ1) = ν (g2−g1)

t2ρ2) = ν (g1−g2)

t(ρ~u) = 0

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

t2E2) = µ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

(35)

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

(36)

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

(37)

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

(38)

### Expansion wave problem: Cavitation test

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

Tliquid = Tvapor = 354.7284K< Tsat

ρvapor = 0.63kg/m3> ρsatvapor, ρliquid= 1150kg/m3> ρsatliquid gsat> gvapor > gliquid

← −~u ~u →

← Membrane

(39)

### 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 104

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

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

x(m) Yv

(40)

### 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 104

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

(41)

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

100 101 102 103

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

100 101 102 103

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 106

p relaxation p−pT−pTg p−pTg

(TvTl)(K)

(42)

### High-speed underwater projectile

With thermo-chemical relaxation No thermo-chemical relaxation

(43)

(44)

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

(45)

### 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) 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

(46)

### 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+ η2−η1

Cp2−Cv2

, B= η1 −η2

Cp2−Cv2

C = Cp2−Cp1

Cp2−Cv2

, D= Cp1−Cv1

Cp2−Cv2

(47)

### 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+ η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

(48)

### 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−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)

Now given the volume fraction for the interface cell C i , we seek a reconstruction that mimics the sub-grid structure of the jump between 0 and 1 in the volume fraction

A mixture-energy-consistent 6-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves. Modelling phase transition in metastable

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

Weak solution for problems with shock &amp; rarefaction waves Interface indicator H I takes value zero away from interfacs, yielding standard compressible Euler equations

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

We have presented a numerical model for multiphase com- pressible ﬂows involving the liquid and vapor phases of one species and one or more inert gaseous phases, extending the

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