A mixture-energy-consistent numerical method for compressible two-phase flow
with
interfaces, cavitation, and evaporation waves
Keh-Ming Shyue
Institute of Applied Mathematical Sciences National Taiwan University
Joint work with Marica Pelanti at ENSTA, Paris Tech, France
Outline
Compressible 2-phase solver for matastable fluids: application to cavitation & flashing flows
1. Motivation
2. Constitutive law for metastable fluid
3. 6-equation single-velocity 2-phase flow model Model with & without heat & mass transfer 4. Stiff relaxation solver
Mixture-energy-consistent method means total energy conservation& pressure consistency
Flashing flowmeans a flow with dramatic evaporationof liquid due to pressure drop
Motivation: Dodecane 2-phase Riemann problem
Riemann data for metastable dodecane modelled by SG EOS Liquid phase: Left-hand side (0 ≤ x ≤ 0.75m)
(ρv, ρl, u, p, αv)L = 2kg/m3, 500kg/m3, 0, 108Pa, 10−8 Vapor phase: Right-hand side (0.75m < x ≤ 1m)
(ρv, ρl, u, p, αv)R = 2kg/m3, 500kg/m3, 0, 105Pa, 1 − 10−8 ,
Liquid Vapor
← Membrane
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, contact, phase, shock
Dodecane 2-phase problem: Sample solution
Allphysical quantities arediscon- tinuous across phase boundary
Expansion wave problem: Cavitation test
Liquid-vapor mixture(αvapor = 10−2) for water with pliquid = pvapor = 1bar
Tliquid= Tvapor = 354.7284K< Tsat
ρvapor = 0.63kg/m3> ρsatvapor, ρliquid= 1150kg/m3> ρsatliquid gsat> gvapor > gliquid
Outgoing velocity u = 2m/s
← −~u ~u →
← Membrane
Expansion wave problem: Sample solution
Cavitation pocket formation &
mass transfer
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 gv− gl(J/kg)
Equilibrium Gibbs free energy inside cavitation pocket
High-speed underwater projectile
Constitutive law: Metastable fluid
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
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
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−Cp2+ η′2−η′1
Cp2−Cv2
, B= η1 −η2
Cp2−Cv2
C = Cp2−Cp1
Cp2−Cv2
, D= Cp1−Cv1
Cp2−Cv2
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−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
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)
Compressible 2-phase flow: 6-equation model
6-equation single-velocity 2-phase model with stiff mechanical relaxation of Saurel et al. (JCP 2009) reads
∂t(α1ρ1) + ∇ · (α1ρ1~u) = 0
∂t(α2ρ2) + ∇ · (α2ρ2~u) = 0
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0
∂t(α1ρ1e1) + ∇ · (α1ρ1e1~u) + α1p1∇ ·~u=µpI(p2−p1)
∂t(α2ρ2e2) + ∇ · (α2ρ2e2~u) + α2p2∇ ·~u=µpI(p1−p2)
∂tα1 + ~u · ∇α1 =µ (p1−p2)
µ: relaxation parameter for volume-transfer rate as p1 → p2; assume stiff µ → ∞ limit
pI: interfacial pressure,
pI = p1/Z1+ p2/Z2
1/Z1+ 1/Z2
, Zi = ρici
6 -equation model (Cont.)
Include conservation of mixture total energy also
∂tE + ∇ · (E~u + p~u) = 0
for purpose of maintaining numerical conservation of total energy
Phasic entropy equations for sk are αkρkTk
dsk
dt = µ (p1−p2)2 Zk
Z1+ Z2
≥0 for k = 1, 2, yielding nonnegative variationof mixture entropy
s = Y1s1+ Y2s2
Model is hyperbolicwith monotone speed of sound cf as c2f = Y1c21+ Y2c22, Yk= αkρk
ρ
6 -equation model: Reduced model
6-equation model approaches to reduced 5-equation model asymptotically as µ → ∞
∂t(α1ρ1) + ∇ · (α1ρ1~u) = 0
∂t(α2ρ2) + ∇ · (α2ρ2~u) = 0
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p= 0
∂tE + ∇ · (E~u +p~u) = 0
∂tα1+ ~u · ∇α1 =
ρ2c22−ρ1c21 ρ1c21/α1 + ρ2c22/α2
∇ ·~u Mixture pressure p determined from total internal energy
ρe = α1ρ1e1(p, ρ1) + α2ρ2e2(p, ρ2)
Model is hyperbolicwith non-monotone sound speed cp: 1
ρc2p = α1
ρ1c21 + α2
ρ2c22
6 -equation model: Reduced model (Cont.)
Volume-fraction equationis differential form ofpressure equilibriumcondition
p1(ρ1, s1) = p2(ρ2, s2)
Assume K = (ρ2c22−ρ1c21) / (ρ1c21/α1+ ρ2c22/α2) < 0, i.e., ρ2c22 < ρ1c21 (phase 1 less compressible)
Compaction effect (K ∇ · ~u > 0)
α1 increases when ∇ · ~u < 0 (compression or shock waves) Relaxation effect (K ∇ · ~u < 0)
α1 decreases when ∇ · ~u > 0 (expansion waves) No effect
α1 remains unchanged when ∇ · ~u = 0 (contacts)
p relaxation: Subcharacteristic condition
Mechanical equilibrium sound speed cp ≤ cf (frozen speed) 1
ρc2p =
2
X
k=1
αk
ρkc2k & ρc2f =
2
X
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
leads tostiffness in equations &
difficulties in numerical solver, e.g., positivity- preserving in volume fraction
6 -equation model: Phasic-total-energy-based
Alternative 6-equation model based on phasic total energy is
∂t(α1ρ1) + ∇ · (α1ρ1~u) = 0
∂t(α2ρ2) + ∇ · (α2ρ2~u) = 0
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1+ α2p2) = 0
∂t(α1E1) + ∇ · (α1E1~u + α1p1~u) + B (q, ∇q) = µpI(p2−p1)
∂t(α2E2) + ∇ · (α2E2~u + α2p2~u) − B (q, ∇q) = µpI(p1−p2)
∂tα1+ ~u · ∇α1 = µ (p1−p2)
B(q, ∇q) is non-conservative product (q: state vector) B= ~u · [Y1∇(α2p2) − Y2∇(α1p1)]
Model is hyperbolicwith monotonic speed of sound cf as well
& is basis for mixture-energy-consistent method
6 -equation model (Cont.)
6-equation model in compact form
∂tq + ∇ · f (q) + w (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
Relaxation scheme
Fractional-step method is employed to solve 6-equation model That is,
1. Non-stiff hyperbolic step
Solve hyperbolic system without relaxation sources
∂tq + ∇ · f (q) + w (q, ∇q) = 0 usingstate-of-the-art solver over time interval ∆t 2. Stiff mechanical relaxation step
Solve system of ordinary differential equations
∂tq = ψµ(q)
with initial solution from step 1 as µ → ∞
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−p2)
∂tα1 =µ (p1−p2)
with initial condition q0 (solution after non-stiff hyperbolic step) & under mechanical equilibrium condition
p1 = p2 = p
Stiff mechanical relaxation step (Cont.)
We find easily
∂t(α1ρ1) = 0 =⇒ α1ρ1 = α01ρ01
∂t(α2ρ2) = 0 =⇒ α2ρ2 = α02ρ02
∂t(ρ~u) = 0 =⇒ ρ~u = ρ0~u0
∂t(α1E1) = µpI(p2−p1) =⇒ ∂t(αρe)1 = −pI∂tα1
∂t(α2E2) = µpI(p1−p2) =⇒ ∂t(αρe)2 = −pI∂tα2
Integrating latter two equations with respect to time Z
∂t(αρe)k dt = − Z
pI∂tαk dt
=⇒ αkρkek−α0kρ0ke0k = −¯pI αk−α0k
or
=⇒ ek−e0k = −¯pI 1/ρk−1/ρ0k
(use αkρk = α0kρ0k) Take p¯I = (p0I + p)/2or p, for example
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
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
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
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
6 -equation model: With heat & mass transfer
6-equation single-velocity 2-phase model with stiff mechanical, thermal, & chemical relaxationsreads
∂t(α1ρ1) + ∇ · (α1ρ1~u) =m˙
∂t(α2ρ2) + ∇ · (α2ρ2~u) =−m˙
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1 + α2p2) = 0
∂t(α1E1) + ∇ · (α1E1~u + α1p1~u) + B (q, ∇q) = µpI(p2−p1) + Q + gIm˙
∂t(α2E2) + ∇ · (α2E2~u + α2p2~u) − B (q, ∇q) = µpI(p1−p2) − Q − gIm˙
∂tα1+ ~u · ∇α1 =µ (p1−p2) + Q qI
+ m˙ ρI
Assume Q = θ (T2−T1), m = ν (g˙ 2−g1)
6 -equation with heat & mass transfer (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,
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
6 -equation with heat & mass transfer (Cont.)
Modified 6-equation 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, Q, −Q, Q/qI]T
ψν = [ ˙m, − ˙m, 0, gIm, −g˙ Im, ˙˙ m/ρI]T
6 -equation with heat & mass transfer (Cont.)
Flow hierarchy in 6-equation model: H. Lund (SIAP 2012)
6-eqns µ → ∞ θν → ∞ µθν → ∞
µθ → ∞
µν → ∞ θ → ∞
ν → ∞
6 -equation with heat & mass transfer (Cont.)
Stiff limits as µ → ∞, µθ → ∞, & µθν → ∞ sequentially
6-eqns µ → ∞ θν → ∞ µθν → ∞
µθ → ∞
µν → ∞ θ → ∞
ν → ∞
Relaxation scheme: With heat & mass transfer
Continue from previous algorithm for 6-equation model with stiff mechanical relaxation, 2 sub-steps are included
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
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)
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) =µpI(p2−p1) + θ (T2−T1)
∂t(α2E2) =µpI(p1−p2) + θ (T1−T2)
∂tα1 =µ (p1−p2) + θ qI
(T2−T1) under mechanical-thermal equilibrium conditons
p1 = p2 = p T1 = T2 = T
Stiff thermal relaxation step (Cont.)
We find easily
∂t(α1ρ1) = 0 =⇒ α1ρ1 = α01ρ01
∂t(α2ρ2) = 0 =⇒ α2ρ2 = α02ρ02
∂t(ρ~u) = 0 =⇒ ρ~u = ρ0~u0
∂t(αkEk) = θ qI
(T2−T1) =⇒ ∂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ρ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ρ01/α1 = 0
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+ p∞1
+ Y2
(γ2−1)Cv2
p+ p∞2
Reduced model: Thermal relaxation step
Reduced model after thermal relaxation step is (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 Mixture entropy ρs satisfy
∂t(ρs) + ∇ · (ρs~u) = θ
1 + peq
qI
(T2−T1)2 T1T2
≥0 Mechanical-thermal equilibriumsound speed cpT satisfies
1
ρc2pT = 1 ρc2p+T
Γ2
ρ2c22
− Γ1
ρ1c21
2 1 α1ρ1cp1
+ 1
α2ρ2cp2
Stiff thermo-chemical relaxation step
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) + θ (T2−T1) + ν (g2−g1)
∂t(α2E2) =µ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
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
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
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
Reduced model: Thermo-chemical relaxation step
Reduced model after thermo-chemical relaxation 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
Mixture entropy ρs satisfies
∂t(ρs) + ∇ · (ρs~u) = ν(g2 −g1)2 Teq
≥0 Speed of sound cpT g satisfies
1
ρc2pT g = 1 ρc2p +T
"
α1ρ1
Cp1
ds1
dp
2
+ α2ρ2
Cp2
ds2
dp
2#
Equilibrium speed of sound: Comparison
Sound speeds follow subcharacteristic condition cpT g ≤cpT ≤cp ≤cf
Limit of sound speed
αlimk→1cf = lim
αk→1cp = lim
αk→1cpT = ck, lim
αk→1cpT 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
5 -equation model: liquid-vapor phase transition
Modelling phase transitionin metastable liquids Saurel et al. (JFM 2008) proposed
∂t(α1ρ1) + ∇ · (α1ρ1~u) =m˙
∂t(α2ρ2) + ∇ · (α2ρ2~u) =−m˙
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0
∂tE + ∇ · (E~u + p~u) = 0
∂tα1+ ∇ · (α1~u) =α1
K¯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)
Expansion wave problem: p-pT relaxation
Mechanical-thermal-equilibrium solution at t = 3.2ms
0 0.2 0.4 0.6 0.8 1
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
104 105
t = 0 t=3.2ms
Pressure(bar)
0 0.2 0.4 0.6 0.8 1
−1
−0.5 0 0.5 1
t = 0 t=3.2ms
Tv− Tl(K)
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
t = 0 t=3.2µs
Vapor volume fraction
0 0.2 0.4 0.6 0.8 1
−1
−0.5 0 0.5 1 1.5
t = 0 t=3.2ms
Vapor mass fraction
Expansion wave problem (Cont.)
Comparison p-, pT - & p-pT g-relaxation solution at t = 3.2ms
0 0.2 0.4 0.6 0.8 1
750 800 850 900 950 1000 1050 1100 1150
p relax pT relax pTg relax
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
p relax pT relax pT relax
Velocity (m/s)
0 0.2 0.4 0.6 0.8 1
103 104 105
p relax pT relax pTg relax
Pressure(bar)
0 0.2 0.4 0.6 0.8 1
−5
−4
−3
−2
−1 0 1x 106
p relax pT relax pTg relax
gv− gl(103K)
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 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 1.2 1.4x 10−4
p relax pT relax pTg relax
Vapor mass fraction
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
Dodecane 2-phase problem: p-pT relaxation
Mechanical-thermal-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
−2 0 2 4 6 8 10 12x 105
t = 0 t=473µs
Tv− Tl(K)
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
Dodecane 2-phase Riemann problem (Cont.)
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
gv− gl(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
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
High-pressure fuel injector: α
v,l= 10
−4Mixture density Mixture pressure
High-pressure fuel injector: α
v,l= 10
−2Mixture density Mixture pressure
High-pressure fuel injector (Cont.)
Vapor volume fraction: αv,l = 10−4 (left) vs. 10−2 (right)
High-pressure fuel injector (Cont.)
Vapor mass fraction: αv,l = 10−4 (left) vs. 10−2 (right)
High-pressure fuel injector: Remark
Numerical solver (to be discussed) uses 400 × 200 uniform Cartesian grid
6-equation single-velocity two-phase model with stiff mechanical relaxation
Consitutive law is stiffened gasequation of state
Model is solved instiff limit towards equilibrium pressure pliquid = pvapor, while admitting different temperatures Tliquid6= Tvapor & entropies sliquid6= svapor
Observation:
Higher αv,l in fluid mixture, higer volume transfer expansion (Is this correct statement ? If so, to what extent)
Fuel injector: p-pT relaxation
Vapor mass fraction: αv,l = 10−4 (left) vs. 10−2 (right)
Fuel injector: p-pT -pT g relaxation
Vapor volume fraction: αv,l = 10−4 (left) vs. 10−2 (right)
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.
