Multiphase CFD Modeling of Flow Boiling & Phase Change Dynamics

Numerical prediction of subcooled & saturated flow boiling heat transfer rates, bubble dynamics, and void fraction distribution in internal pipe flows.

Overview of Internal Flow Boiling Regimes

Flow boiling occurs when a liquid flows through a heated tube and undergoes a continuous thermodynamic phase change into vapor along its length. The flow pattern transitions through distinct visual and hydrodynamic regimes—driven by mass flux, heat flux, saturation pressure, and surface orientation.

Orientation & Gravity Effects: Flow regimes differ fundamentally by pipe angle. In vertical pipes, forces are axisymmetric. In horizontal pipes (such as the validation cases presented below), gravity stratifies the phases, creating distinct Stratified (Smooth & Wavy) flow patterns where vapor accumulates along the top crown while liquid wets the bottom invert.

Sequential Flow Regimes (Inlet to Outlet)

Single-Phase Liquid Flow

Both Horizontal & Vertical

Subcooled liquid enters the pipe; heat is transferred strictly via single-phase forced convection prior to local thermodynamic saturation.

Bubbly Flow (Subcooled & Saturated)

Both (Asymmetric in Horizontal)

Discrete vapor bubbles nucleate along the heated wall at the Onset of Nucleate Boiling (ONB) and detach into the fluid core.

Plug / Slug Flow

Intermittent / Transitional

As vapor fraction increases, discrete bubbles coalesce into large vapor plugs. In horizontal tubes, buoyancy forces these plugs to sweep along the upper wall section.

Stratified Flow (Smooth & Wavy)

Horizontal Pipe Dominant

Gravity forces full phase separation: liquid flows along the tube bottom while vapor flows above. As vapor velocity increases, interfacial shear creates stratified waves that roll along the liquid-gas interface, washing up the side walls.

Churn Flow

Highly Transient • Vertical Dominant

A highly oscillatory, turbulent transitional regime where large vapor slugs break down due to shear forces. Occurs primarily in vertical flows as liquid repeatedly falls back and gets swept upward; in horizontal pipes, it manifests as unsteady wavy-annular flow.

Annular Flow

Both (Thin Film Upper-Wall in Horizontal)

A continuous high-velocity vapor core forms in the center of the tube, forcing a liquid film to coat the internal pipe circumference.

Mist Flow / Dryout

Post-Dryout Transition

The liquid wall film evaporates completely (Dryout), leaving entrained liquid droplets suspended in a continuous vapor stream transitioning to single-phase vapor.

Validation & Benchmark Summary for Flow Boiling Regime

To demonstrate the capability of multiphase CFD in predicting phase-change phenomena, two distinct internal flow boiling benchmark cases were evaluated using R134a and water across subcooled and saturated regimes in a uniformly heated horizontal tube.

The numerical models were validated against the established empirical correlation of Liu & Winterton (1991) and the fundamental heat transfer formulations in Incropera et al.

Key Performance Metrics Evaluated:

Heat Transfer Coefficients: Agreement with published empirical correlations across varying mass fluxes and heat fluxes.
Phase Change Dynamics: Accurate prediction of onset of nucleate boiling (ONB) and spatial void fraction evolution along the tube length.
Thermal Performance: Capturing subcooled fluid heating and saturated boiling transitions under fixed wall heat flux conditions.

Flow Boiling Validation & Analysis

CFD Performance Benchmarking & Empirical Model Formulation

Geometry and Operating Conditions

The validation analysis is performed on straight horizontal circular tubes subjected to uniform outer wall heat fluxes. Fluid enters each domain as a subcooled liquid and undergoes saturated flow boiling along the length. Two distinct working fluids (Refrigerant R134a and Water) are modeled across specific geometric scales and operating parameters:

Case 1: Refrigerant R134a

Geometry: ID = 10.88 mm, OD = 12.7 mm, L = 1.0 m
Inlet Condition: Liquid, \(u_{\text{in}} = 0.235 \text{ m/s}\)
Saturation State: \(5^\circ\text{C}\) at \(349.6 \text{ kPa}\)
Heat Flux: \(q_{s,o}'' = 8,567 \text{ W/m}^2\) (outer wall)

Case 2: Water

Geometry: ID = 10.88 mm, OD = 12.7 mm, L = 1.0 m
Inlet Condition: Liquid, \(u_{\text{in}} = 0.5 \text{ m/s}\)
Saturation State: \(425 \text{ K}\) at \(500 \text{ kPa}\)
Heat Flux: \(q_{s,o}'' = 42,835 \text{ W/m}^2\) (outer wall)

Empirical Formulation Framework

To benchmark CFD predictions, the comprehensive Kandlikar correlation (from Incropera) and the Liu & Winterton formulation are paired with core energy balances and slip ratio models to compute two-phase heat transfer coefficients, local wall temperatures, exit vapor quality, and phase distribution dynamics.

Kandlikar Correlation

Calculates saturated flow boiling by taking the maximum value between the convective dominant (\(h_a\)) and nucleate dominant (\(h_b\)) mechanisms.

\[ h = \max \left( h_a, \, h_b \right) \]

\(h_a = h_{sp} \left[ 0.6683 \left(\frac{\rho_l}{\rho_v}\right)^{0.1} x^{0.16}(1-x)^{0.64}f(Fr) + 1058 \, Bo^{0.7}(1-x)^{0.8} G_{s,f} \right]\)
\(h_b = h_{sp} \left[ 1.136 \left(\frac{\rho_l}{\rho_v}\right)^{0.45} x^{0.72}(1-x)^{0.08}f(Fr) + 667.2 \, Bo^{0.7}(1-x)^{0.8} G_{s,f} \right]\)

Liu & Winterton Model

Asymptotic combination model integrating forced liquid convection and nucleate pool boiling mechanisms.

\[ h_{\text{TP}} = \sqrt{(F \cdot h_L)^2 + (S \cdot h_{\text{pool}})^2} \]

Energy Balance & Zivi Slip Model

Global energy conservation defines exit vapor quality (\(\bar{X}_{\text{exit}}\)), while Zivi's slip ratio accounts for phase velocity differences to compute accurate exit void fractions (\(\alpha\)).

\[ \bar{X}_{\text{exit}} = \frac{Q_{\text{total}}}{\dot{m} \, h_{fg}} \]

\[ S = \left(\frac{\rho_l}{\rho_v}\right)^{1/3}, \quad \alpha = \left[ 1 + \left(\frac{1 - \bar{X}}{\bar{X}}\right) \left(\frac{\rho_v}{\rho_l}\right) S \right]^{-1} \]

Validation and Results

The performance metrics below benchmark numerical CFD results against established empirical models, including the Liu & Winterton model, the Kandlikar correlation, and global energy balances coupled with the Zivi Slip Ratio model for both R134a and Water flow boiling test cases.

Parameter / Metric	Liu & Winterton Model	Kandlikar Model	Global Energy Balance & Zivi Slip Model	CFD Results	Deviations
HTC – R134a (\(\text{W}/\text{m}^2\cdot\text{K}\))	2,000	2,800	—	1,800	10%
Wall Superheat – R134a (\(\text{K}\))	5	8.5	—	5.5	—
HTC – Water (\(\text{W}/\text{m}^2\cdot\text{K}\))	12,900	9,400	—	8,600	9%
Wall Superheat – Water (\(\text{K}\))	3.9	5.7	—	7.6	—
Vapor Quality (\(x\)) – R134a	—	—	6.1%	3.5%	—
Vapor Quality (\(x\)) – Water	—	—	1.9%	1.5%	—
Vapor Void Fraction – R134a	—	—	54%	55%	2%
Vapor Void Fraction – Water	—	—	50%	65%	30%

Transient Visualization: Void Fraction and Heat Transfer Coefficient Evolution

Transient Evolution of Void Fraction and Heat Transfer Coefficient Contours

The following animations illustrate the time-dependent phase behavior and thermal performance inside the heated tube. The top contours demonstrate the spatial development of vapor volume fraction (α), capturing the transition from single-phase liquid to slug and annular regimes over time. Below, the corresponding heat transfer coefficient (h) contours show local thermal spikes driven by micro-layer evaporation and bubble agitation. Comparisons between R134a and Water highlight significant differences in transient phase distribution and local cooling capacity due to their distinct latent heat, density ratios, and surface wetting characteristics.

Conclusion: Due to the complex physics of flow boiling and the absence of a single, unified CFD model, every project requires case-by-case validation against established empirical correlations before applying the methodology to complex target geometries.