Skip to content

Introduction to Thermo-fluid Dynamics of Two-phase Flow ​

  • Fundamentals of two-phase flow
    • Two-phase flow definition
  • Fundamentals of computational simulation code for two-phase flow
  • Fundamentals of mass, momentum, and energy conservation equations for two-phase flow
    • Homogeneous flow model, drift-flux model, two-fluid model
  • Fundamentals of several constitutive equations
    • Flow regime map
    • Void fraction
    • Interfacial area concentration
    • Frictional pressure drop
    • Interfacial force
  • Two-phase flow modeling

Fundamentals of two-phase flow ​

What is two-phase flow ​

Two-phase flow is composed of two immiscible mediums.

  • Gas-liquid two-phase flow
  • Solid-gas two-phase flow
  • Solid-liquid two-phase flow
  • Liquid-liquid two-phase flow
gas-liquid flow
ClassTypical regimesGeometryConfigurationExamples
Separated flowsFilm flowtable of two-phase flow1Liquid film in gas; Gas film in liquidFilm condensation; Film boiling
Annular flowtable of two-phase flow2Liquid core and gas film; Gas core and liquid filmFilm boiling; Boilers
Jet flowtable of two-phase flow3Liquid jet in gas; Gas jet in liquidAtomization; Jet condenser
Mix or Transactional flowsCap, Slug or Churn-turbulent flowtable2 of two-phase flow1Gas pocket in liquidSodium boiling in forced convection
Bubbly annular flowtable2 of two-phase flow2Gas bubbles in liquid film with gas coreEvaporators with wall nucleation
Droplet annular flowtable2 of two-phase flow3Gas core with droplets and liquid filmSteam generator
Bubbly droplet annular flowtable2 of two-phase flow4Gas core with droplets and liquid film with gas bubblesBoiling nuclear reactor channel
Dispersed flowsBubbly flowtable3 of two-phase flow1Gas bubbles in liquidChemical reactors
Droplet flowtable3 of two-phase flow2Liquid droplets in gasSpray cooling
Particulate flowtable3 of two-phase flow3Solid particles in gas or liquidTransportaion of powder
flow regimeMishima-Ishii regime map
  • Bubbly-to-Slug: ⟨α⟩→0.3

Why is two-phase flow so difficult ​

Different Density:

density of two-phase flow

Peculiar characteristics of two-phase flow dynamics ​

Thermal and Kinematic Non-Equilibrium

  • Driving Mechanism: Natural circulation flow, forced convective flow (laminar or turbulent flow)
  • Flow Regime: Dispersed (bubbly, slug, churn flow), separated (annular or stratified flow)
  • Flow Orientation: Upward, downward, inclined, counter current, elbow, sudden expansion & contraction
  • Flow Channel Scale: Micro, mini, conventional (or medium-size), large
  • Phase Change: Wall nucleation & condensation, interfacial heat transfer, non-condensable gas effect
  • Wall Heat Transfer: One-component with phase change, two-component without phase change
  • Wall Effect: CHF, bubble nucleation, hydrophilic or hydrophobic
  • Characteristic Length Scale: Kolmogorov scale, bubble size, sub-channel size, flow channel size
  • Pressure and Temperature: Room temperature-to-nuclear reactor core temperature
  • Gravity: Microgravity-to-elevated gravity
  • Flow Instability: Geysering, chugging, density-wave oscillation
  • Pressure Wave and Shock Wave: Water hammer, steam explosion
  • Critical Flow: Choking flow
  • Flow Structure: Internal flow, external flow
  • Coupling between Two-Phase Flow & Structure: Flow induced vibration

Fundamentals of computational simulation code for two-phase flow ​

Key parameters in two-phase flow analysis ​

  • Mass conservation eq.
  • Momentum conservation eq.
  • Energy conservation eq.

Simplified computational code structure ​

procedure of solving two-phase flow

No. of balance eq. Single-phase flow: 3 Two-phase flow: 3(mixture model), 6 (=3×2), or others? Assumptions

  • Relative velocity b/w gas & liquid phases
  • Temperature difference b/w gas & liquid phases Two-phase flow model
  • Homogeneous flow model
  • Mixture model
  • Two-fluid model Pros. & Cons.
  • Numerical instability
  • Available constitutive eqs.

Averaging methods ​

  • Time domain
    • Instantaneous
    • Time average
  • Space domain
    • Local
    • Control volume (or surface) average
    • 1D (Nuclear system analysis)
    • Porous (Steam generator analysis)
    • Sub-channel (Sub-channel analysis)
    • 3D (Detailed analysis of flow field)

Fundamentals of mass, momentum, and energy conservation equations for two-phase flow ​

Instantaneous single-phase flow ​

Reynolds number: (inertia force / viscous force)

N-S equation (Instantaneous formulation):

ρDvDt=−∇p−∇⋅Cμ+ρgma=F

F = (Surface force) + (Volumetric force)

  • (Surface force) = (Pressure)+(Shear stress)
  • (Volumetric force) = (Gravity)

a = (Local acceleration) + (Convective acceleration)

Time-averaged single-phase flow ​

A×B―≠A¯×B¯A×B―=A¯×B¯+α

N-S equation (Time-averaged formulation):

ρDv¯Dt=−∇p¯−∇⋅(Cμ―+CT―⏟Reynolds stress)+ρg.

Reynolds stress: Momentum transfer due to turbulence

Turbulence models

  • Mixing length model
  • k−ε model

Time-averaged two-phase flow ​

Two-phase flow

  • Instantaneous formulation of gas and liquid phases
ρDvDt=−∇p−∇⋅Cμ+ρg
  • Time-averaged formulation
    • Kinematic condition (Relative velocity)
      • Homogeneous model
      • Mixture model (1 mixture momentum eq.)
        • Slip model
        • Drift-flux model
    • Two-fluid model (2 momentum eq.)
  • Thermal condition (Temperature difference)
    • Thermal equilibrium (1 mixture energy eq.)
    • Thermal non-equilibrium (2 energy eq.)
  • Evaporation (Tg>Tf)
  • Condensation (Tg<Tf)

Drift-flux model ​

4 equations:

  • Mixture continuitv equation
∂ρm∂t+∇⋅(ρmvm)=0
  • Continuity equation for phase 2
∂α2ρ2――∂t+∇⋅(α2ρ2――vm)=Γ2−∇⋅(α2ρ2――V2m)
  • Mixture momentum equation
∂ρmvm∂t+∇⋅(ρmvmvm)=−∇pm+∇⋅[C―+C~T+C~D]+ρmgm+Mm

where

CD≡−∑k=12αkρkVkmVkm,Vkm≡vk−vmvm≡∑k=12αkρkvkρm,ρm≡∑k=12αkρk
  • Mixture Energy equation
Vkj≡vk−jj=∑k=12αkvk

Total number of unknown variables: 27: Flow regime

ρk――,p¯¯k,T¯¯k,i^k,Ti――,σ――,ρm,vm,pm,Vkm,C―,CT,Mm,αk,im,q―,qT,Φmμ,Φmσ,Φmi,Γ2

Two-fluid model [core] ​

Totally 6 equations:

  • Continuity equation
∂αkρk――∂t+∇⋅(αkρk――vk^)=Γk
  • Momentum equation
∂αkρk――vk^∂t+∇⋅(αkρk――vk^vk^)=−∇(αkpk――)+∇⋅[αk(C――k+CkT)]+αkρk――g^k+Mk
  • Energy equation
  • 3 jump conditions
∑k=12Γk=0. gas-liquid interface

Total number of unknown variables: 36: Flow regime

αk,ρk――,vk^,Γk,pk――Ck――,TkT,Cki――,Mik,MmH,pki――,ik^,qk――,qkT,aiqk′′――Tk――,Ti――,H21――,σ¯¯,ai,psat

With HEM, it becomes 3 equations

Summary of two-phase flow formulation ​

2-field diagram

Fundamentals of several constitutive equations ​

One-dimensional computational code structure ​

1D diagram

Important parameters in two-phase flow analysis [core] ​

  • Flow regime map
    • Basic idea of flow regime transition criteria
  • Void fraction
    • Basic idea of drift-flux model
  • Interfacial area concentration
    • Significance of interfacial transfer terms
  • Frictional pressure drop
  • Interfacial force
    • Drag force
    • Non-drag force (Lift force, wall lubrication force etc.)
  • scaling design
  • mechanistic modeling
  • instrumentation development
  • experiment
  • computational calculation

Flow regime map ​

1D diagram

The flow regime depends on the void fraction

Void fraction ​

Drift velocity, which is generated by the force balance b/w buoyancy force and drag force:

vgj=vg−j=(1−α)vr
vgj=vg−j=vg−(jg+jf)=vg−[αvg+(1−α)vf]=(1−α)(vg−vf)=(1−α)vr
  • ⟨⋅⟩: × (void fraction)
  • ⟨⟨⋅⟩⟩: Area-averaging

One-dimensional (or area-averaged) drift-flux model[core] :

⟨⟨vg⟩⟩≡⟨jg⟩⟨α⟩=C0⟨j⟩+⟨⟨vgj⟩⟩

where C0 is distribution parameter and ⟨⟨vgj⟩⟩ is the drift velocity:

C0≡⟨αj⟩⟨α⟩⟨j⟩,⟨⟨vgj⟩⟩≡⟨αvgj⟩⟨α⟩
vgj=vg−j⇒αvgj=αvg−αj⇒⟨αvgj⟩=⟨αvg⟩−⟨αj⟩⇒⟨αvg⟩=⟨αj⟩+⟨αvgj⟩⇒⟨αvg⟩⟨α⟩=⟨αj⟩⟨α⟩+⟨αvgj⟩⟨α⟩⇒⟨αvg⟩⟨α⟩=⟨αj⟩⟨α⟩⟨j⟩⟨j⟩+⟨αvgj⟩⟨α⟩

In Homogeneous flow (Uniform gas distribution, C0=1, & no relative velocity, ⟨⟨vgj⟩⟩=0 m/s ) Gas velocity:

⟨⟨vg⟩⟩≡⟨jg⟩⟨α⟩=⟨j⟩

Mixture volumetric flux Center of volumetric velocity

⟨α⟩=⟨jg⟩C0⟨j⟩+⟨⟨vgj⟩⟩=⟨jg⟩/⟨j⟩C0+⟨⟨vgj⟩⟩/⟨j⟩

In HEM: ⟨⟨vgj⟩⟩=0,C0=1

C0=⟨αj⟩⟨α⟩⟨j⟩=1A∫αj dA1A∫α dA1A∫j dA=1(If α is constant)

Distribution parameter C0:

{=1flat surface / uniform α>1with surface tension<1wall peaking

Dix model for C0:

C0=β[1+(1β−1)b]

where

b=(ρgρl)0.1
β=(Gxρg)/(Gxρg+G(1−x)ρf)
  • For regular pipe:
C0=1.2−0.2(ρgρf)0.5.
  • For rectangular channel / duct:
C0=1.35−0.35(ρgρf)0.5.
  • For rod bundle:
C0=1.10−0.10(ρgρf)0.5.
  • For Subcooled boiling flow:
C0=1.2−0.2(ρgρf)0.5(1−e−18⟨α⟩)⟹C0=1.2 If ⟨α⟩=0.

C0: Channel geometry dependent

Drift velocity ⟨⟨vgj⟩⟩:

Suggested by Lahey & Moody:

vgj=2.9[(ρf−ρgρf2)σg]0.25
  • Bubbly:
⟨⟨vgj⟩⟩=2(Δρgσρf2)1/4(1−⟨α⟩)1.75
  • Slug:
⟨⟨vgj⟩⟩=0.35(Δρgσρf2)0.5
  • Churn:
⟨⟨vgj⟩⟩=2(Δρgσρf2)1/4

Development of void fraction model

  1. HEM, and 20% overprediction
  2. Slip model: Slip ratio s=
s=vgvf=jg/αjf/(1−α)=(1−α)jgαjf→sjfjg=1α−1→α=11+sjfjg
  1. Drift-flux model

Interfacial area concentration ​

  • Mass equation
∂αkρk∂t+∇⋅(αkρkvk)=Γ¯k=aimk
  • Momentum equation
∂αkρkvk∂t+∇⋅(αkρkvkvk)=−αk∇pk+∇⋅αk(Ckμ+CkT)+αkρkg+vkiΓk+Mjk−∇αk⋅Ci

Energy equation

∂αkρkHk∂t+∇⋅(αkρkHkvk)=−∇⋅αk(q――k+qkT)+αkDkDtpk+HkiΓk+aiqik−1+Φk

(Interfacial transfer term)=( Interfacial area concentration ) x ( Driving potential )

ai=6αDsm

Frictional pressure drop ​

Lockhart-Martinelli method[core] :

(−ΔpΔz)2ϕ,fric =Φf2(−ΔpΔz)1ϕ,fric,f Darcy friction factor 

where (single phase)

(−ΔpΔz)1ϕ,fric,f=λ1Dρfjf22,λ={64/Ref for laminar flow 0.316/Ref0.25 for turbulent flow 

and λ is Dancy friction factor; and Reynolds number:

Ref≡ρfjfDμf

For Two-phase multiplier, here's the Chiskolm equation:

Φf2=1+CX+1X2,C=20for gas & liquid turbulent flows

and

X≡(−ΔpΔz)1ϕ,fric,f(−ΔpΔz)1ϕ,fric,gis Martinelli parameter

[Final-term]

pressure drop for cylinderπD24p1=πD24p2+πDΔzτw⇒−p2−p1Δz=4Dτw⇒−ΔpΔz=4Dτw
GasLiquidC: Chisholm coefficient
TurbulentTurbulent20
TurbulentLaminar12
LaminarTurbulent10
LaminarLaminar5
Rek{<1000,Laminar≥2000,Turbulent,k=g or f

Quiz - Void fraction ​

Superficial gas velocity =0.50 m/s, Superficial liquid velocity =0.30 m/s : What is the void fraction for air-water flow in a vertical 7 cm pipe under atmospheric pressure? Assume churn flow regime.

  • Flow regime
  • Flow direction
  • Flow channel geometry
  • Channel diameter
⟨jg⟩⟨α⟩=C0⟨j⟩+⟨⟨vgj⟩⟩

where

C0=1.2−0.2ρgρf=1.2

and

⟨⟨vgj⟩⟩=2(Δρgσρf2)0.25=0.23 m/s

Thus

⟨α⟩=⟨jg⟩C0⟨j⟩+⟨⟨vgj⟩⟩=0.501.2(0.50+0.30)+0.231=0.42⟨⟨vg⟩⟩=⟨jg⟩⟨α⟩=0.500.42=1.2 m/s

Quiz - Flow regime transition ​

Consider air-water flow in a vertical 7 cm pipe under atmospheric pressure. Derive bubbly-to-slug flow transition condition with critical void fraction of 0.30 .

  • Flow direction
  • Flow channel geometry
  • Channel diameter
⟨jg⟩⟨α⟩=C0⟨j⟩+⟨⟨vgj⟩⟩

where

C0=1.2−0.2ρgρf=1.2⟨⟨vgj⟩⟩=2(Δρgσρf2)0.25(1−⟨α⟩)1.75

and

⟨jf⟩=1−C0⟨α⟩crit C0⟨α⟩crit ⟨jg⟩−⟨⟨vgj⟩⟩C0⟨jf⟩=1−(1.2−0.2ρg/ρf)⟨α⟩crit (1.2−0.2ρg/ρf)⟨α⟩crit ⟨jg⟩−2(Δρgσρf2)0.25(1−⟨α⟩)1.75(1.2−0.2ρg/ρf)

where

⟨α⟩crit =0.3.

Quiz - Interfacial area concentration ​

Consider air-water flow in a vertical 7 cm pipe under atmospheric pressure. Assume the bubble size given by twice the size of Laplace length scale. Assume uniformly distributed spherical bubbles. Calculate the interfacial area concentration for bubbly flow with the void fraction of 0.20.

La=σΔρg=0.0027 m⟨Db⟩=2La=2×0.0027=0.0055 m⟨α⟩=16π⟨Db⟩3n

Thus

⟨ai⟩=π⟨Db⟩2n=6⟨α⟩⟨Db⟩=6×0.200.0055=220 m−1

Quiz2 ​

A saturated steam-water mixture at 290∘C vertically flows in a rectangular channel with a width of 40 mm and a gap length of 2.4 mm . The volumetric gas and liquid flow rates are 5.0×10−5 m3/s and 3.0×10−5 m3/s, respectively. The distribution parameter is calculated by C0=1.35−0.35(ρg/ρf)0.5, whereas the drift velocity is calculated by ⟨⟨vgj⟩⟩=0.35[(ρf−ρg)gDhρf]0.5. The bubble Sauter mean diameter is calculated by ⟨DSm⟩=2La. Laplace length La is given by [σ(ρf−ρg)g]0.5.

Here, ρg(=39.2 kg/m3) and ρf(=732 kg/m3) are the densities of gas and liquid, respectively. g(=9.8 m/s2) is the gravitational acceleration. Dh is the hydraulic equivalent diameter. σ(0.0167 N/m) is the surface tension. Use three digits for significant digits.

  1. What is the value of the distribution parameter?
C0=1.35−0.35(ρgρf)0.5=1.27.

1.21−1.33

  1. What is the value of the drift velocity in  m/s?
Dh=4Ap=2×0.04×0.00240.04+0.0024 m≈0.00453 m⟨⟨vgj⟩⟩=0.35[(ρf−ρg)gDhρf]0.5≈0.0717 m/s.

0.0681−0.0753

  1. What is the value of superficial gas velocity in  m/s?
⟨jg⟩=QgA=5.0×10−50.04×0.0024 m/s=0.521 m/s.

0.495−0.547

  1. What is the value of gas velocity in  m/s?
⟨jf⟩=QfA=3.0×10−50.04×0.0024 m/s=0.3125 m/s⟨⟨vg⟩⟩=C0⟨j⟩+⟨⟨vgj⟩⟩=C0(⟨jg⟩+⟨jf⟩)+⟨⟨vgj⟩⟩≈1.27×(0.521+0.3125)+0.0717≈1.13 m/s.

1.07−1.19

  1. What is the value of the void fraction?
⟨α⟩=⟨jg⟩⟨⟨vg⟩⟩≈0.5211.13≈0.461.

0.438−0.484

  1. What is the value of the bubble Sauter mean diameter in  m?
⟨DSm⟩=2La=2[σ(ρf−ρg)g]0.5=2×0.0167(732−39.2)×9.8 m≈0.00314 m.

0.003−0.0033

  1. What is the value of the interfacial area concentration in 1/m?
ai=6⟨α⟩DSm≈882 m−1.

838−926

  1. What is the value of the two-phase mixture density in  kg/m3?
ρm=⟨α⟩ρg+(1−⟨α⟩)ρf=0.461×39.2+(1−0.461)×732 kg/m3≈412 kg/m3.

392−433

  1. What is the value of the two-phase mass flux in  kg/m2s?
G=ρg⟨jg⟩+ρf⟨jf⟩≈39.2×0.521+732×0.3125≈249 kg/m2s.

237−262

  1. What is the value of the flow quality?
⟨x⟩=GgG=ρg⟨jg⟩G≈39.2×0.521249≈0.0819.

0.0778−0.086

Released under the MIT License.