• 沒有找到結果。

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

N/A
N/A
Protected

Share "A mixture-energy-consistent numerical method for compressible two-phase ﬂow with interfaces, cavitation, and evaporation waves"

Copied!
63
0
0

(1)

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

(2)

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

(3)

### Cavitating Richtmyer-Meshkov problem

Gas volume fraction Mixture pressure

t=0ms

t=2ms

t=3.1ms

t=6.4ms

t=8.6ms

(4)

### High-speed underwater projectile

With thermo-chemical relaxation No thermo-chemical relaxation

(5)

### High-pressure fuel injector

With thermo-chemical relaxation No thermo-chemical relaxation

(6)

### 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)γk1 + ηk Helmholtz free energy ak = ek−Tksk

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

(7)

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

(8)

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

(9)

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

(10)

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

(11)

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

(12)

### Homogeneous relaxation model (Cont.)

µ, θ, ν → ∞: 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, &

isnonzero only at 2-phase mixture

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

(13)

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

(14)

### Homogeneous relaxation model (Cont.)

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

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(15)

### Homogeneous relaxation model (Cont.)

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

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

µθ → ∞

µν → ∞ θ → ∞

ν → ∞

(16)

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

Assume frozen thermal & chemical relaxation,HRMreduces to

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

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

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

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

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

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

(17)

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

Assume frozen thermal & chemical relaxation,HRMreduces to

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

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

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

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

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

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

Take formal asymptotic expansion ansatz q = q0+ εq1+ · · ·

Find equilibrium equation for q0 as µ = 1/ε → ∞ (ε → 0+)

(18)

### Reduced model: Asymptotic analysis as µ → ∞

Define material derivative d dt = ∂

∂t + ~u · ∇

Rewrite energy & volume-fraction equations in form dp1

dt + ρ1c21∇ ·~u = ρ1c2I1 α1

1

ε(p2−p1) dp2

dt + ρ2c22∇ ·~u = ρ2c2I2 α2

1

ε(p1−p2) dα1

dt = 1

ε(p1−p2)

Assume formal asymptotic expansion as α1 = α01+ εα11+ · · ·

pk = p0k+ εp1k+ · · · for k = 1, 2

(19)

### Reduced model: Asymptotic analysis (Cont.)

We get d

dt p01+ εp11+ · · ·

+ ρ1c21∇ ·~u = ρ1c2I1

α1

1 ε

 p02−p01

+ ε p12−p11

+ · · · d

dt p02+ εp12+ · · ·

+ ρ2c22∇ ·~u = ρ2c2I2

α2

1 ε

 p01−p02

+ ε p11−p12

+ · · · d

dt α10+ εα11+ · · ·

= 1 ε

 p01−p02

+ ε p11−p12

+ · · ·

(20)

### Reduced model: Asymptotic analysis (Cont.)

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

p01 = p02 =p0 =⇒ p0I = p0, c0Ik2 = c0k2 O(1)

dp0

dt + ρ01c012∇ ·~u = ρ01c012

α01 p12−p11 dp0

dt + ρ02c022∇ ·~u = ρ02c022

α02 p11−p12 dα01

dt = p11−p12

(21)

### Reduced model: Asymptotic analysis (Cont.)

Subtracting former two equations, we find

01c012 −ρ02c022

∇ ·~u = ρ01c012

α1002c022 α02

!

p12−p11

i.e.,

01

dt = p11−p12 = ρ02c022 −ρ01c012 ρ01c01201+ ρ02c02202

!

∇ ·~u

(22)

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

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

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

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

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

tα1+ ~u · ∇α1 =

 ρ2c22−ρ1c21 ρ1c211 + ρ2c222



∇ ·~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)

(23)

### p relaxation: Subcharacteristic condition

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

ρc2p = X2 k=1

αk

ρkc2k & ρc2f = X2 k=1

αkρkc2k

(24)

### p relaxation: Subcharacteristic condition

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

ρc2p = X2 k=1

αk

ρkc2k & ρc2f = X2 k=1

αkρkc2k

0 0.2 0.4 0.6 0.8 1

101 102 103 104

frozen p relax

αwater

cp&cf

Non-monotonic cp

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

(25)

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

ρc2pT = 1 ρc2p +T

 Γ2

ρ2c22 − Γ1 ρ1c21

2  1 α1ρ1cp1

+ 1

α2ρ2cp2



(26)

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

ρc2pT g = 1 ρc2p +T

"

α1ρ1

Cp1

ds1

dp

2

2ρ2

Cp2

ds2

dp

2#

(27)

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

(28)

### 5 -equation model: liquid-vapor phase transition

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

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+ 1 qI

Q+ 1 ρI

˙ m K¯s=

1

Ks1 + α2

Ks2

1

, Ksι = ριc2ι qI =

Ks1 α1

+ Ks2 α2

   Γ1

α1

2

α2



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

Ks1 α1

+Ks2 α2

   c21 α1

+ c22 α2



, m = ν(g˙ 2−g1)

(29)

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

(30)

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

(31)

### HRM: Phasic-internal-energy-based

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

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

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

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

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

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

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

(32)

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

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)

(33)

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

(34)

### Definition (mixture-energy consistent)

(i) Mixture total energy conservation consistency E0 = E0,C

where E0 = (α1E1)0+ (α2E2)0 (ii) Relaxed pressure consistency

e0,C = α1e1



p,(α1ρ1)0 α1



+ α2e2



p,(α2ρ2)0 α2

 ,

where e0,C = E0,C(ρ~u)0·(ρ~0u)0

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

(35)

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

(36)

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

(37)

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

(38)

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

(39)

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

(40)

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

(41)

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

(42)

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

(43)

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

(44)

### p relaxation: Algebraic approach

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

(45)

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

(46)

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

(47)

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

(48)

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

(49)

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

(50)

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

(51)

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

(52)

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

+ Y2

2−1)Cv2

p+ p2

(53)

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

(54)

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

(55)

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

(56)

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

(57)

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

(58)

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

(59)

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

(60)

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

(61)

### Vapor-bubble compression: pT g relaxation results

Vapor mass fraction Mixture pressure

t=0.4ms

t=0.8ms

t=1.2ms

(62)

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

(63)

## Thank you

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

Our model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates Our grid system is a time-varying grid. Extension of the model to

To date we had used PSO and successfully found optimal designs for experiments up to 8 factors for a mixture model, nonlinear models up to 6 parameters and also for more involved

In Section 3, the shift and scale argument from [2] is applied to show how each quantitative Landis theorem follows from the corresponding order-of-vanishing estimate.. A number

In this paper, we propose a practical numerical method based on the LSM and the truncated SVD to reconstruct the support of the inhomogeneity in the acoustic equation with

Then, we recast the signal recovery problem as a smoothing penalized least squares optimization problem, and apply the nonlinear conjugate gradient method to solve the smoothing

Then, we recast the signal recovery problem as a smoothing penalized least squares optimization problem, and apply the nonlinear conjugate gradient method to solve the smoothing

This kind of algorithm has also been a powerful tool for solving many other optimization problems, including symmetric cone complementarity problems [15, 16, 20–22], symmetric