Computational fluid dynamics and population balance modelling of nucleate boiling of cryogenic liquids: Theoretical developments

Abstract

The main focus in the analysis of pool or flow boiling in saturated or subcooled conditions is the basic understanding of the phase change process through the heat transfer and wall heat flux partitioning at the heated wall and the two-phase bubble behaviours in the bulk liquid as they migrate away from the heated wall. This paper reviews the work in this rapid developing area with special reference to modelling nucleate boiling of cryogenic liquids in the context of computational fluid dynamics and associated theoretical developments. The partitioning of the wall heat flux at the heated wall into three components – single-phase convection, transient conduction and evaporation – remains the most popular mechanistic approach in predicting the heat transfer process during boiling. Nevertheless, the respective wall heat flux components generally require the determination of the active nucleation site density, bubble departure diameter and nucleation frequency, which are crucial to the proper prediction of the heat transfer process. Numerous empirical correlations presented in this paper have been developed to ascertain these three important parameters with some degree of success. Albeit the simplicity of empirical correlations, they remain applicable to only a narrow range of flow conditions. In order to extend the wall heat flux partitioning approach to a wider range of flow conditions, the fractal model proposed for the active nucleation site density, force balance model for bubble departing from the cavity and bubble lifting off from the heated wall and evaluation of nucleation frequency based on fundamental theory depict the many enhancements that can improve the mechanistic model predictions. The macroscopic consideration of the two-phase boiling in the bulk liquid via the two-fluid model represents the most effective continuum approach in predicting the volume fraction and velocity distributions of each phase. Nevertheless, the interfacial mass, momentum and energy exchange terms that appear in the transport equations generally require the determination of the Sauter mean diameter or interfacial area concentration, which strongly governs the fluid flow and heat transfer in the bulk liquid. In order to accommodate the dynamically changing bubble sizes that are prevalent in the bulk liquid, the mechanistic approach based on the population balance model allows the appropriate prediction of local distributions of Sauter mean diameter or interfacial area concentration, which in turn can improve the predictions of the interfacial mass, momentum and energy exchanges that occur across the interface between the phases. Need for further developments are discussed.

Keywords

Nucleate boiling wall heat flux partitioning model computational fluid dynamics population balance model

Introduction

Cryogenic liquids are encountered in many industries. For the propulsion and thermal management for space systems, liquid hydrogen and liquid oxygen are usually used as fuel for RS-25 and J-2X rocket engines.¹ In the field of medical applications, liquid nitrogen and liquid helium are utilised as coolants for superconductors such as can be found in commercial MRI, MRT, NMRI and Large Hadron Collider.^2,3 Natural gas, which is considered to be one of the clean and largest energy sources, requires its liquefied form to be stored and transported economically and safely in tanks and through pipelines in the oil and gas industries. Nowadays, all of the aforementioned technological developments have predicated an increased interest in the better understanding of the boiling of cryogenic liquids. It is apparent that boiling of ordinary liquids is rather different from boiling of cryogenic liquids. For an ordinary liquid such as water, it can remain as a liquid within a range of about 100 K. Nevertheless, cryogenic liquids only remain in their liquid form within a much narrower atmospheric liquid range (from melting point to boiling point at atmospheric pressure) of about 20 K (e.g. He – 6 K, H₂ – 7 K, N₂ – 14 K, liquefied natural gas – 20 K). Handling cryogenic liquids thus presents many challenges during the operations of two-phase cryogenic systems. For the case of the transportation of cryogenic liquids, small bubbles that may be generated from the boiling of cryogenic liquids due to the presence of inevitable heat leaks along the pipe can merge or coalesce to form larger bubbles and possibly reaching the same diameter of the pipe during flow boiling. These so-called Taylor bubbles can subsequently cause the blockage of pipelines which may result in the destruction of the pipes and valves as well as possible cavitation erosion to the cryogenic systems. Also, heat absorption from the surroundings into the storage tanks of liquefied natural gas may initiate the boiling off of liquefied natural gas during pool boiling. Such situation is undesirable since additional gas generation can dramatically increase the pressure within these storage tanks being already at their pressurised state.

Pool boiling or flow boiling is essentially a phase change process in which vapour bubbles are formed on a heated wall. The main difference between pool boiling and flow boiling is that the former takes place from a stagnant liquid on top of a submerged heated surface in a reservoir of liquid while the latter is produced through heating the liquid through conduits such as tubes or pipes. There are essentially different modes and regimes of boiling that can be observed in pool boiling and flow boiling.

In pool boiling, the different modes and regimes of boiling can be characterised by a typical saturated boiling curve depicted in Figure 1, which illustrates the relationship between the heat flux and the wall superheat temperature. As the heat flux and wall superheat temperature increases, nucleate boiling occurs at point a, which denotes the onset of nucleate boiling (ONB). The detachment of vapour bubbles from the heated surface is observed to be of spherical shape bubbles between points a and b but of slug and column shape bubbles between points b and c in the nucleate boiling regime. The heat flux finally reaches point c, which denotes the critical heat flux (CHF). After point c, the boiling process proceeds to transition boiling. The heat flux decreases as the wall superheat increases to point d, the Leidenfrost point, which signifies the minimum wall temperature required for film boiling.

Figure 1.

Typical pool boiling curve and flow patterns.

In flow boiling, there are, however, two distinct regimes, known as subcooled and saturated nucleate boiling, depending on the bulk temperature of the liquid moving past the heating surface. Figure 2 shows the evolution of the two-phase flow in a vertical heated pipe. During subcooled boiling, the temperature of the heating surface is above the saturation temperature, ebullition occurs in the immediate vicinity of the heating surface. Bubbles are produced and attached at the heating surface. As the bulk liquid is continually heated downstream, saturated nucleate boiling occurs with bubbles beginning to detach from the heating surface into the liquid stream. Further downstream, complex topologies of the two-phase flows persist throughout the entire vertical heated pipe such as flow transition regimes from bubbly to slug/churn to annular and disperse flows.

Figure 2.

Evolution of the steam–water flow in a vertical heated pipe.

Based on the aforementioned observations prevalent in nucleate pool boiling or subcooled and saturated nucleate boiling flow of cryogenic liquids, it is clear that considerations for modelling requires the boiling phenomena to be predicted not only at the heating wall but also away from the heating wall.⁴ In the framework of computational fluid dynamics (CFD), modelling subcooled or saturated nucleate boiling can be broadly classified into two categories: (i) heat transfer mechanistic modelling at the heating wall and (ii) two-phase flow and bubble behaviours in the bulk liquid away from the heating wall.

In general, the first category encompasses the development of suitable models to predict the heat transfer at the heated wall via partitioning the wall heat flux. On the basis of the review by Warrier and Dhir,⁵ there are three approaches whereby the heat transfer rate can be determined. The first approach employs empirical correlations to predict the total wall heat flux given the wall superheat temperature in regions covering single-phase convection, partial nucleate boiling and fully developed boiling.^6–9 The second approach considered utilising appropriate correlations to predict the different heat transfer mechanisms being considered in the wall heat flux partitioning and the fraction of the total heat flux that is transferred through each of these mechanisms.^10,11 The third approach, which is the primary focus in this review, adopts the mechanistic modelling of the heat transfer mechanisms occurring during boiling, which can be readily employed for both the prediction of the wall heat flux as well as the partitioning of the wall heat flux between the liquid and gas phases. Almost all CFD investigations on flow boiling are mainly based on the use of these mechanistic models.^12–14 The model partitions the wall heat flux into three heat flux components: single-phase convection, transient conduction and evaporation.

The second category involves the treatment of the fluid flow in the bulked liquid phase. Based on the interpenetrating media framework, the two-fluid model is considered to be the most detailed and accurate macroscopic formulation to resolve the thermal and hydrodynamic characteristics of the two-phase systems. In the so-called Eulerian–Eulerian approach, the field equations in the two-fluid model comprise two sets of conservation equations of mass, momentum and energy. One set is dedicated to resolve the gas phase while the other is solved for the liquid phase. Through an ensemble average of the transport equations, interphase exchange terms appear in the averaged equations, which represent the mass, momentum and energy transfers through the interface between the phases. These interphase exchange terms determine the rate of phase changes as well as mechanical and thermal non-equilibrium between phases. Most important of all, these terms provide the necessary closure to the two-fluid model formulation.

The prediction of local bubble sizes in cryogenic liquids is also expected to be strongly influenced by the complex bubble behaviours near the heated wall. In essence, complex mechanistic behaviours of bubble coalescence and bubble collapse due to condensation in subcooled nucleate boiling exist during the boiling of cryogenic liquids.^15,16 Physical insights strongly indicate that the numerical solutions of two-phase cryogenic systems can be significantly improved through consideration of the bubble mechanistic behaviours in determining the local bubble size such as the Sauter mean bubble diameter, which appears prominently in the interphase exchange terms of mass, momentum and energy. Because of the importance of accurately determining the Sauter mean bubble diameter, which strongly governs the inter-phase heat and mass transfer through the interfacial area as well as the momentum drag, the potential to extend the modelling of boiling of cryogenic liquids in order to better predict the non-uniform bubble size distribution via the population balance approach^17,18 is of enormous significance.

Description of the various models or correlations developed for the prediction of heat transfer and wall heat flux partitioning along with the formulation of governing equations for boiling of cryogenic liquids and population balance model is treated in the following sections.

Wall heat flux partition model

According to Kurul and Podowski¹² and Podowski,¹³ the mechanistic modelling of local near-wall heat transfer applicable to saturated and subcooled boiling can be accounted based on partitioning the total wall heat flux into three components

Q w = Q e + Q q + Q c

(1)

where Q_e is the rate of heat transfer that is directly employed to convert the liquid into vapour at the heated surface, Q_q is the quenching component associated with heating the liquid, which comes into contact with the heated surface after the departure of bubbles, and Q_c is the single-phase heat transfer component corresponding to the non-boiling regions of the heated surface. In equation (1), the wall heat flux component Q_e is evaluated based on the consideration of microlayer evaporation occurring at the nucleate boiling region which can be determined from

Q e = N a f ρ g V b h fg

(2)

where N_a is the active nucleation site density, f is the nucleation frequency,

V b = π D_{b}^{3} / 6

is the volume of departing bubbles and

h fg

is the latent heat. The wall heat flux component contributed by surface quenching Q_q is modelled as transient conduction to a semi-infinite medium with a step change in temperature at the heated surface, i.e.

Q q = (\frac{2}{\sqrt{π t w}} \sqrt{k l ρ l C pl}) A q (T w - T l)

(3)

where t_w is the waiting time and A_q is the bubble area that denotes the fraction of the wall area that is in contact with the fresh liquid and is cooled down by transient conduction. Here, T_l refers to the temperature of the bulk liquid which may be saturated or subcooled. In equation (3), the area A_q is calculated from

A q = N a (K \frac{π D_{b}^{2}}{4})

(4)

where K is the proportional constant for bubble diameter of influence. The single-phase heat transfer component Q_c is generally taken due to natural or convection which can be defined as

Q c = h l (1 - A q) (T w - T l)

(5)

Alternatively, Benjamin and Balakrishnan¹⁹ have considered the heat to be removed from the heated surface by the boiling liquid according to

Q w = Q me t g f + Q q t w f + Q c

(6)

where t_g is the growth time of the bubble and f is given by

1 / (t g + t w)

. In equation (6), the weighted sum of the first two wall heat flux components is employed because the two modes are complimentary to each other. The heat absorbed by the evaporating microlayer

Q me

during the growth time is given by

Q me = N a V l ρ l h fg

(7)

The difference in partitioning the total wall heat flux through this particular approach is the consideration of the volume of microlayer evaporated V_l which is given by

V l = \frac{\int_{0}^{t g} A l (t) Δ δ (t) d t}{\int_{0}^{t g} d t}

(8)

where

A l (t)

and

δ (t)

denote the dry area under the bubble and microlayer thickness. The rest of the wall heat flux components in equation (6) follow the same formulations in the aforementioned mechanistic modelling; the quenching and single-phase heat transfer rates are depicted by equations (3) and (5).

Closure relations

The respective wall heat flux components that have been formulated in the mechanistic models require the determination of three unknown parameters: active nucleation site density (N_a), bubble departure diameter (D_b) and nucleation frequency (f ). In order that the multi-dimensional modelling framework can become an effective tool to predict the boiling heat transfer of cryogenics liquids, it has to be supplemented with accurate closure relations. Driven by practical applications and industrial needs, the prevailing approach has been thus far to adopt empirical correlations for the evaluation of N_a, D_b and f. Useful empirical correlations derived from experimental data for a range of conditions and fluids have been developed. Some of the most common correlations summarised in Tables 1 to 3 are further described in the preceding sections.

Table 1.

Empirical correlations for active nucleation site density.

Contributors	Correlations
Gaertner and Westwater²⁰	$N a ~ Q_{w}^{1 / m}$	m = 0.48
Hsu and Graham²¹		0.48 < m < 1
Westwater and Kirby²²		m = 0.73 (glass), 0.48 < m < 1 (metals)
Mikic and Rohsenow²⁴	$N a ~ \frac{C}{D_{c}^{m}}$ $D c = \frac{4 σ T sat}{ρ g h fg Δ T \sup}$	m = 6.5, C = $D_{c, \max}^{m}$ m = 3, C = 4 × 10^-12
Johov²³		m = 6.5, C = $D_{c, \max}^{m}$ m = 3, C = 4 × 10^-12
Bier et al.²⁵	$\ln N a = \ln (N \max) [1 -] (\frac{D c}{D c, \max}) m]$	m = 0.42 (Freon-115) m = 0.26 (Freon-11)
Cornwel and Brown²⁶	$N a ~ T_{\sup}^{m}$	m = 4.5
Lemmert and Chawla²⁷	$N a = [n Δ T \sup] m$	n = 185, m = 1.805
Koncar et al.²⁸		n = 185, m = 1.805
Kurul and Podowski¹²		n = 210, m = 1.805
Krepper et al.^29,30	$N a = n [\frac{Δ T \sup}{Δ T refN}] m$	n = 0.8 × 10⁶, ΔT_refN = 10 K (Bartolomej et al.³¹) n = 3 × 10⁷ (DEBORA 1–2³³) n = 5 × 10⁶ (DEBORA 3–7³³)
Kirichenko³³	$N a = n [\frac{ρ g h fg Δ T \sup}{σ T sat}] m$	n = $1 \times 10 - 7$ , m = 2, $P / P cr \geq 0.04$ n = $6.25 \times 10 - 6$ , m = 3, $P / P cr < 0.04$
Kocamustafaogullari and Ishii³⁴	$N a = \frac{1}{D_{b}^{2}} [(\frac{2 σ T sat}{ρ g h fg Δ T eff}) / (0.5 D b)] - 4.4 f (ρ *)$	$ρ * = Δ ρ / ρ g$ , $Δ ρ = ρ l - ρ g$ $f (ρ ) = 2.157 e - 7 ρ - 3.2 (1 + 0.0049 ρ ) 4.13$ $Δ T eff = S Δ T \sup$
Wang and Dhir³⁵	$N a = 5 \times 10 - 31 (1 - \cos ϕ) D_{c}^{- 6.0}$
Benjamin and Balakrishnan¹⁹	$N a = 218.8 \Pr 1.63 (\frac{1}{γ}) ω - 0.4 Δ T_{\sup}^{3}$	$γ \equiv (\frac{k w ρ w C p, w}{k l ρ l C p, l}) 1 / 2 ω = 14.5 - 4.5 (\frac{R a P}{σ}) + (\frac{R a P}{σ}) 0.4$
Hibiki and Ishii³⁶	$N a = 4.72 \times 105 {1 - \exp (- ϕ 2 / 4.17)} \times [\exp {2.5 \times 10 - 6 f (ρ +) (Δ T \sup ρ g h fg / 2 σ T sat)} - 1]$	$ρ + = \log 10 (Δ ρ / ρ g)$ $f (ρ +) = - 0.01064 + 0.48246 ρ +$ $- 0.22712 ρ + 2 + 0.05468 ρ + 3$
Basu et al.³⁷	$N a = 0.34 \times 104 (1 - \cos ϕ) Δ T_{\sup}^{2.0}$	$Δ T ONB$ < $Δ T \sup$ < 15 K
Basu et al.³⁷	$N a = 0.34 (1 - \cos ϕ) Δ T_{\sup}^{5.3}$	$Δ T \sup$ ≥ 15 K

Table 2.

Empirical correlations for bubble departure diameter.

Contributors	Correlations
Fritz [38]	$D b (D bF) = 0.0208 ϕ [\frac{σ}{g (ρ l - ρ g)}] 0.5$
Zuber³⁹	$D b = [\frac{6 σ}{g (ρ l - ρ g)} \frac{k Δ T \sup}{Q w}] 1 / 3$
Han and Griffith⁴⁰	$D b = 0.843 ϕ [\frac{2 σ}{g (ρ l - ρ g)}] 0.5$
Cole and Shulman⁴¹	$D b = [\frac{σ}{g (ρ l - ρ g)}] 0.5 P - 1$
Cole⁴²	$D b = ϕ [\frac{σ}{g (ρ l - ρ g)}] 0.5 \frac{ρ l C p, l Δ T \sup}{ρ g h fg}$
Cole and Rohsenow⁴³	$D b = 1.5 \times 10 - 4 [\frac{σ}{g (ρ l - ρ g)}] 0.5 (\frac{ρ l C p, l T sat}{ρ g h fg}) 1.25$
Tolubinskiy and Kostanchuk⁴⁶	$D b = D ref \exp (- \frac{Δ T sub}{Δ T refd})$
Ünal⁴⁸	$D b = \frac{2.42 \times 10 - 5 P 0.709 a}{\sqrt{b φ}}$	$a = \frac{(Q w - h l Δ T sub) 1 / 3 k l}{2 C 1 / 3 h fg \sqrt{π k l / ρ l C p, l ρ g}} \sqrt{\frac{k w ρ w C p, w}{k l ρ l C p, l},}$ $C = \frac{h fg μ l [C p, l / (0.013 h fg \Pr 1.7)] 3}{[σ / (ρ l - ρ g) g] 0.5},$ $b = Δ T sub / 2 (1 - ρ g / ρ l),$ $φ = {(\frac{u l}{0.61}) 0.47 for u l \geq 0.61 m / s 1.0 for u l < 0.61 m / s$
Kocamustafaogullari and Ishii³⁴	$D b = 0.0012 (\frac{ρ l - ρ g}{ρ l}) 0.9 D bF$
Farajisarir⁴⁶	$D b = 10.02 \times 109 Ja - 1.65 T - 1.65$	$T * = \frac{T w - T l}{Δ T \sup}$
Morel et al.⁵¹	$D b = 2.42 \times 10 - 5 P 0.709 \frac{A}{\sqrt{bk}}$	$A = \frac{Δ T \sup}{2 ρ g h fg} \sqrt{\frac{k w ρ w C p . w}{π},}$ $b = {\frac{Δ T sub}{2 (1 - ρ g / ρ l)} for St < 0.0065 \frac{1}{2 (1 - ρ g / ρ l)} for St > 0.0065$ $St = \frac{Q w}{0.0065 ρ l C p, l u l Δ T sub}$ $k = max (1, (\frac{u l}{0.61}) 0.47)$

Table 3.

Empirical correlations for nucleation frequency.

Contributors	Correlations
Cole⁵²	$f = [\frac{4 g (ρ l - ρ g)}{3 D b ρ l}] 0.5$
Zuber³⁹	$f = \frac{0.59}{D b} [\frac{σ g (ρ l - ρ g)}{ρ_{l}^{2}}] 0.25$
Kocamustafaogullari and Ishii³⁴	$f = \frac{1.18}{D b} [\frac{σ g (ρ l - ρ g)}{ρ_{l}^{2}}] 0.25$
Hatton and Hall⁵⁵	$f = \frac{3}{π η l} [\frac{16 k l σ T sat}{(h fg ρ g) 2 D b D c}] 2$
Ivey⁵⁶	$f = 0.9 (\frac{g}{D b}) 0.5 forhydrodynamicregion$
	$f = 0.44 (\frac{g}{D_{b}^{3}}) 0.25 fortransitionregion$
	$f = \frac{C}{D_{b}^{2}} forthermodynamicregion$
Stephan⁵⁷	$f = \frac{1}{π} \sqrt{\frac{g}{2 D b}} [1 + \frac{4 σ}{D_{b}^{2} ρ g g}] 0.5$

Active nucleation site density

At the heated surface, not all the nucleation sites will be activated as the temperature of the surface exceeds the saturation liquid temperature at the local pressure. For nucleation sites to be activated, they are governed by the distributions of cavities on the wall surface, heater and liquid properties and contact angle between liquid and the wall. Active nucleation site density plays an important role in the evaluation of the heat flux components of the wall heat flux partition model. A large amount of research has been carried out to correlate suitable relationships of the active nucleation site density such as shown in Table 1.

Gaertner and Westwater²⁰ have formulated the active nucleation site density based on dissolution of nickel salts in water and observed the deposition on the copper-heater surface during boiling. An estimate of the active nucleation site density was obtained by counting the number of holes in the deposited nickel layer. The correlation for the active nucleation site density is expressed as a function of the wall heat flux. In their investigation, the exponent m matched the experimental data with a value of 0.48. Hsu and Graham²¹ further reported that the range of the exponent m varies between 0.48 and 1. Westwater and Kirby²² concluded that the range of the exponent m is 0.73 for glass and between 0.66 and 0.5 for metal. This demonstrated that the rate of activation was fairly steady regardless of the finish of the heated a wide range of cavity sizes exists.

Johov²³ related, however, the active nucleation site density to the critical cavity diameter (D_c) of the activated sites on the heating surface. For a pressure of 1 bar, the exponent m and C are taken to be 3 and $4 \times 10 - 12$ , respectively. Similarly, Mikic and Rohsenow²⁴ related the active site density to the cavity size existing on the surface where the exponent m is an empirical constant taken to be 6.5 and C is a variable depending on the diameter of the largest cavity present on the heated surface ( $D c, \max$ ). Bier et al.²⁵ also correlated the active nucleation site density as a function of the critical cavity diameter where $N \max$ is the value corresponding to D_c = 0. The value of the exponent m was found to vary depending on the surface preparation procedure. It is taken to be 0.42 and 0.26 for Freon-115 and Freon-11 boiling on copper surfaces.

Cornwell and Brown²⁶ performed experiments for the pool boiling of saturated water on a copper surface at 1 bar. Using electron microscope, a systematic analysis was carried out to determine the active nucleation site density. They have correlated the active nucleation site density as a function of the wall superheat temperature (ΔT_sup) with a value of 4.5 for the exponent m. The functional dependence on the local wall superheat temperature was justified by assuming that only conical cavities were present on the heated surface and before any nucleation could occur, vapour are required to be trapped in the cavities.

In a similar dependence on the wall superheat temperature, Lemmert and Chawla²⁷ investigated the active nucleation site density for experiments on pool boiling of saturated water with an additional consideration of a constant n. Končar et al.²⁸ have proposed to adopt values of n and m as 185 and 1.805 while Kurul and Podowski¹² suggested values for n and m as 210 and 1.805.

Recently, Krepper et al.^29,30 proposed that the active nucleation site density has a strong influence on the wall superheat temperature, but has a small influence on the gas volume fraction and barely any influence on the liquid temperature. They suggested a similar correlation to the active nucleation site density correlation of Lemmert and Chawla²⁷ with additional consideration of the reference temperature ( $Δ T refN$ ). Here, the exponent m has a physical significance which is related to the pressure in MPa. Referring to Kolev,³¹ the correlation was derived to match different data sets with vastly different parameter values. They justified the high dependency of the active nucleation site density to the microscale topography of the boiling surface, and in turn strong dependency on the processes that have been adopted to produce the finishing of the heated surface, which can be very diverse and not specifically controlled in most boiling experiments. In their modelling framework, best results were achieved by taking values of $n = 0.8 \times 106 m - 2$ and $Δ T refN = 10 \circ K$ for the experiment of Bartolomej and Chanturiya,³² where high pressure water was used, and values ranging from $n = 3 \times 107 m - 2$ to $n = 5 \times 106 m - 2$ for the DEBORA²⁹ experiments carried out in France for the boiling of R12.

Kirichenko et al.³³ proposed a correlation for the active nucleation site density for the nucleate boiling of cryogenic liquids. He has shown that the correlation has a dependency on not only the wall superheat temperature such has those proposed by the investigations above but also include the surface tension. The values of n and m depend on the pressure of the system: $n = 1 \times 10 - 7$ and $m = 2$ when $P / P cr \geq 0.04$ $n = 6.25 \times 10 - 6$ and $m = 3$ when $P / P cr < 0.04$ where $P cr$ is the critical pressure of the fluid.

Kocamustafaogullari and Ishii³⁴ performed a parametric study to correlate the active nucleation site density for the existing data of pool boiling in the literature as a function of bubble departure diameter as well as the surface tension. They have also extended the correlation to forced convection nucleate boiling by the use of an effective superheat ( $Δ T eff$ ) where a suppression factor (S) is employed to reduce the actual wall superheat temperature.

Wang and Dhir³⁵ performed a systematic study of the effect of wettability on the density of active nucleation sites. Based on the pool boiling of water at 1 bar on vertical copper surfaces, the contact angle (ϕ) was varied from 18° to 90° by controlling the degree of oxidation of the surface. The shape of the cavities and the cumulative number density of them were obtained with an electron microscope, while the number of active sites was determined from still pictures. Hence, they correlated the active nucleation site density as a function of the contact angle as well as the critical cavity diameter. It should be noted that the correlation is valid for $D c < 5.8 \times 10 - 6 m$ which corresponds to $Δ T \sup > 11.2 \circ C$ for atmospheric water.

Benjamin and Balakrishnan¹⁹ investigated the nucleation site density during nucleate pool boiling through experiments on saturated pure liquids at low to moderate heat fluxes. The effect of surface–liquid interaction during the boiling phenomena was investigated on the nucleation sites. Different working mediums were used namely distilled water, carbon tetrachloride, n-hexane and acetone as liquid with the heated surface made of stainless steel and aluminium with different surface finishes. The active nucleation site density has been correlated as a function of the Prandtl number (⪻), surface–liquid interaction parameter (γ), dimensionless surface roughness parameter (ω) and the wall superheat temperature. For the dimensionless surface roughness parameter, R_a is the arithmetic average roughness which is defined as the average values of the peaks and valleys on the surface. The correlation is valid for $1.7 < ⪻ < 5$ , $4.7 < γ < 93$ , $0.02 mm < R a < 1.17 mm$ , $5 K < Δ T \sup < 25 K$ , $2.2 < ω < 14$ and $13 \times 10 - 3 N / m < σ < 59 \times 10 - 3 N / m$ .

Hibiki and Ishii³⁶ mechanistically modelled the active nucleation site density by accounting the size and cone angle distributions of cavities existing on the heated surface. The active nucleation site density was correlated a function of the critical cavity diameter, contact angle and wall superheat temperature. This correlation gave fairly good predictions over the flow conditions: $30 \circ < ϕ < 90 \circ$ , $1.0 \times 104 sites / m 2 < N a < 1.5 \times 1010 sites / m 2$ , $0.101 MPa < P < 19.8 MPa$ and $0 kg / m 2 s < G < 886 kg / m 2 s$ .

Basu et al.³⁷ conducted experiments of flow boiling at subcooled conditions using a nine-rod (zircalloy-4) bundle and a flat plate with a copper surface to correlate the partitioning of the heat flux at the heated surface. Using a CCD camera, individual activated nucleation sites were manually accounted. The active nucleation site density was correlated with a dependence on the contact angle and wall superheat temperature. The range of parameters covered for the flat plate test surface were: $2.5 W / cm 2 < Q w < 96 W / cm 2$ , $30 \circ < ϕ < 90 \circ$ , $6.6 K < Δ T sub, in < 52.5 K$ and $124 kg / m 2 s < G < 886 kg / m 2 s$ .

Bubble departure diameter

Most empirical correlations of bubbles departure diameter have been formulated as a function of the contact angle, Jacob number or thermo-hydraulic parameters in consideration of the particular boiling process. Table 2 depicts the various correlations for the bubble departure diameter.

Fritz³⁸ considered the static equilibrium between the adhesive and buoyancy forces acting on the bubble in order to predict the bubble departure diameter whereby the contact angle has been considered.

Zuber³⁹ assumed the existence of a thin superheated thermal layer near the surface for the occurrence of bubble growth. The correlation can be expressed as a function of the wall superheat temperature and wall heat flux.

Han and Griffith⁴⁰ derived a modified model based on existing bubble growth model. The bubble departure diameter could be obtained through the consideration of the bubble growth derivative set to zero, which has a similar form to the correlation of Fritz.³⁸

Cole and Shulman⁴¹ proposed that the bubble departure diameter being inversely proportional to the absolute pressure of the boiling system. Subsequently, Cole⁴² proposed another correlation as a function of the contact angle along with the consideration of the wall superheat temperature. Cole and Rohsenow⁴³ modified the correlation of Cole⁴² by replacing the wall superheat temperature with a critical temperature at saturation.

Kocamustafaogullari and Ishii³⁴ correlated the bubble departure diameter with their experimental water data around atmospheric pressure. This correlation represents a modified version for higher pressures based on the evaluation of the correlation of Fritz³⁸ at low pressures.

Kirichenko et al.⁴⁴ and Bald⁴⁵ formulated the correlation for bubble departure diameter as a function to the critical cavity diameter due to the edge-break mechanism since small contact angles ( $ϕ < 7.5 o$ ) have been observed in most metallic surfaces in cryogenic liquids.

Tolubinsky and Kostanchuk⁴⁶ employed the local bulk temperature for at subcooled conditions in the evaluation of the bubble departure diameter. Krepper et al.²⁹ employed this correlation for their modelling and defined the reference temperature ( $Δ T refD$ ) as 53 K and 45 K for atmospheric pressure and higher pressures, respectively. The reference diameter ( $d ref$ ) is obtained at atmospheric pressure. For experiments of Bartolomej and Chanturiya³² and DEBORA²⁹ $d ref$ ranged between 0.24 mm and 1.3 mm. Farajisarir⁴⁷ proposed a non-dimensional model for bubble departure diameter based on his examination of forced convective boiling flow at atmospheric pressure. The correlation is formulated as function of the Jakob number defined as $Ja = ρ l C p, l Δ T \sup / ρ g h fg$ and the local bulk temperature.

Ünal⁴⁸ investigated different sets of experiments to determine empirical correlations for bubble growth rate, maximum bubble diameter and maximum bubble growth time of subcooled boiling water. For calculating bubble departure diameter, the Rohsenow⁴⁹ superposition method was employed to determine the nucleate flow boiling heat transfer. At subcooled conditions. A similar method to Levenspiel⁵⁰ was then applied to calculate the heat transfer coefficient. The analysis showed that the proposed equation for bubble departure diameter has a good agreement with most of experiments within ±30% uncertainty. The correlation is valid for $0.1 MPa < P < 17.7 MPa$ , $0.47 W / cm 2 < Q w < 10.64 W / cm 2$ , $3 K < Δ T sub, in < 86 K$ , $0.08 m / s < u liquid < 9.15 m / s$ and $0.08 mm < D d < 1.24 mm$ . Morel et al.⁵¹ modified the correlation of Ünal⁴⁸ in order that it can be applied right up to the saturated condition.

Nucleation frequency

A number of investigations have been carried out to ascertain the nucleation frequency. Table 3 lists the various correlations for the nucleation frequency.

Cole⁵² investigated the boiling phenomena in the vicinity of the critical heat flux. On the basis of the assumption that when successive bubbles leave the surface as they touch and coalesce, the nucleation frequency multiply by the bubble departure diameter at lift off can be taken to be equal to the rate that bubbles leave the surface.

Zuber³⁹ assumed the quotient of bubble departure diameter divided by the growth time can be taken to equal to the bubble rise velocity in a gravitational field. In the development of the nucleation frequency, he further assumed that bubble waiting time is almost equal to bubble growth time. Kocamustafaogullari and Ishii³⁴ applied the same correlation for their investigative study in a boiling channel albeit with a higher nucleation frequency due to a larger constant being taken.

For cryogenic liquids, the correlations usually take the form similar to those of Cole⁵² or Zuber.³⁹ Correlation by McFadden and Grassman⁵³ has shown to agree well with the experimental data of liquid nitrogen. Kirichenko et al.⁴⁴ have developed a semi-empirical correlation, based on the experimental data, through correcting the contact angle of the detached bubble. In the developed pool nucleate boiling, Jin et al.⁵⁴ have nonetheless observed that the detachment frequency of coalesced bubbles remained almost unchanged with high wall heat fluxes and different heating surface materials in cryogenic liquids. For simplicity, the nucleation frequency can be regarded as a constant.

Hatton and Hall⁵⁵ correlated the nucleation frequency in terms of the square power of bubble departure diameter and critical cavity diameter with the experimental data from pool boiling of water at pressure of 0.0162, 0.0662 and 0.101 MPa. The correlation has been formulated based on the assumption that the waiting time was negligible in comparison to the growth time.

Ivey⁵⁶ proposed three correlations, relating the nucleation frequency for three different regions: (1) hydrodynamic region in which the buoyancy and drag forces predominate, correlated with water and methanol data; (2) transition region where the buoyancy, drag and surface tension forces are of the same order, correlated with water, methanol, isopropanol and carbon tetrachloride data; and (3) thermodynamic region where bubble growth dominates, correlated with water data.

Stephan⁵⁷ investigated the heat transfer coefficient in boiling systems and reported that nucleation frequency could be estimated from a correlation by Malenkov.⁵⁸

Mechanistic closure

It is recognised that not one single dedicated correlation is expected to be universally valid for all types of flow conditions to determine the three important parameters of active nucleation site density, bubble departure diameter and nucleation frequency in the wall heat flux partitioning model. In order to remove the application uncertainty of correlations empirically determined, a more mechanistic approach in determining the closure for the wall heat flux partitioning model may be required to improve the CFD predictions for cryogenic liquids.

Fractal model

The fractal distribution of sizes of active cavities on heated surfaces being considered is strongly dependent on not only the wall superheat as has been typified in most correlations for the active nucleation site density but also on other parameters including the liquid subcooled temperature, fractal dimension, minimum and maximum active cavity sizes, contact angle and physical properties of the adjacent fluid including the bulk velocity. From a mechanistic consideration, this model demonstrates enormous potential in appropriately determining the active nucleation site density distribution. Yeoh et al.^59,60 have recently applied the fractal model to predict the active nucleation site density on the heated surface for a specific vertical subcooled flow boiling configuration.

According to Xiao and Yu,⁶¹ the generated active cavities on heated surfaces for pool or flow nucleate boiling can be considered to be similar to the existence of pores in porous media. The cumulative number of active cavities with diameters equal to and greater than a particular active cavity diameter (D_c) can be described by

N a (D \geq D c) = (\frac{D c, max}{D c}) d f

(9)

The total number of active nucleation sites per unit area (sites/cm²) from the smallest to the largest active cavity diameter can thus be obtained as

N a, tot = (\frac{D c, max}{D c, min}) d f

(10)

where

D c, min \leq D c \leq D c, max

in which

D c, min

and

D c, max

are smallest and largest active cavity diameters, and

d f

is the fractal dimension.

The smallest and largest active cavity diameters for nucleation site distribution can be evaluated based on Hsu⁶² according to

D c, min = 2 \frac{δ}{C 1 (Δ T \sup + Δ T sub)} \times {Δ T \sup - \sqrt{Δ T_{\sup}^{2} - \frac{4 ζ C 3}{δ} (Δ T \sup + Δ T sub)}}

(11)

D c, max = 2 \frac{δ}{C 1 (Δ T \sup + Δ T sub)} \times {Δ T \sup + \sqrt{Δ T_{\sup}^{2} - \frac{4 ζ C 3}{δ} (Δ T \sup + Δ T sub)}}

(12)

where

ζ = 2 σ T sat / ρ g h fg

C 1 = (1 + \cos θ) / \sin θ

C 3 = 1 + \cos θ

and δ is the thickness of the thermal boundary layer which can be expressed as

δ = k l / h l

where h_l is the single-phase heat transfer coefficient for natural or force convection.

The fractal dimension $d f$ of nucleation sites for a boiling system can be derived from Yu and Cheng⁶³ as

d f = \frac{\ln [\frac{1}{2} (\frac{\bar{D} c, max}{D c, min}) 2]}{\ln (\frac{D c, max}{D c, min})}

(13)

where

\bar{D} c, max

is the averaged value over all the maximum active cavities, which is determined based on Xiao and Yu⁶¹ as

\bar{D} c, max = \frac{1}{Δ T \sup} \int_{T sat}^{T w} D c, \max (T w) d T w = \frac{1}{Δ T \sup} \sum_{i = 1}^{n} D c, \max (T w, i) δ T w = (1 / n) \sum_{i = 1}^{n} D c, \max (T w, i)

(14)

where

δ T w

is assumed as a constant (e.g.

δ T w = 0.1

);

T w, i = T sat + i (δ T w)

with i = 1, 2, … , n; and

n = T \sup / δ T w

Force balance model for bubble detachment from cavity site and bubble lift off

A force balance model developed in the context of CFD framework by Yeoh and Tu⁶⁴ has demonstrated that the bubble growth and its maximum detached bubble size contributed significantly towards the growth of void fraction and heat transfer within the bulk liquid. The development of the force balance model concentrates primarily on the various forces that influence the growth of a bubble in the directions parallel and normal to a heating surface. These forces are formulated according to the studies performed by Klausner et al.⁶⁵ and Zeng et al.^66,67 As illustrated in Figure 3, the forces acting on the bubble in the x- and y-directions for a vertical heating surface are

\sum F x = F sx + F dux + F L + F h + F cp

(15)

and

\sum F y = F sy + F duy + F qs + F b

(16)

while the forces acting on the bubble in the x- and y-directions for a horizontal heating surface are

\sum F x = F sx + F dux + F qs

(17)

and

\sum F y = F sy + F duy + F L + F h + F cp + F b

(18)

Figure 3.

Schematic views illustrating the forces acting on a growing vapour bubble along the vertical and horizontal heating surfaces.

In the above equations, F_s is the surface tension force,

F du

is the unsteady drag due to asymmetrical growth of the bubble and the dynamic effect of the unsteady liquid such as the history force and the added mass force, F_L is the lift force created by the wake of the preceding departed bubble in pool boiling or shear lift force in flow boiling, F_h is the force due to the hydrodynamic pressure in flow boiling which is absent in pool boiling,

F cp

is the contact pressure force accounting for the bubble being in contact with a solid rather than being surrounded by liquid,

F qs

is the quasi steady-drag in the flow direction and F_b is the buoyancy force. The various forces acting on the bubble during flow flowing along vertical and horizontal heating surfaces are illustrated in Table 4. It should be noted that g indicates the gravitational acceleration, α, β and θ_i are respectively the advancing, receding and inclination angles, d_w is the surface/bubble contact diameter, d is the bubble diameter at the wall and d_r is the curvature radius of the bubble at the reference point on the surface which is taken to be equivalent to five times of the bubble radius:

d r ~ 5 d

Table 4.

Forces acting on a growing bubble on the heated wall.

Bubble forces
$F sx = - d w σ \frac{π}{(α - β)} [\cos β - \cos α]$
$F sy = - d w σ \frac{π (α - β)}{π 2 - (α - β) 2} [\sin α + \sin β]$
$F l = \frac{1}{8} C L ρ l Δ U 2 π d 2$
$F h = \frac{9}{8} ρ l Δ U 2 \frac{π d_{w}^{2}}{4}$
$F cp = \frac{π d_{w}^{2}}{4} \frac{2 σ}{d r}$
$F qs = 3 C D μ l Δ U π d$
$F b = \frac{1}{6} π d 3 (ρ l - ρ g) g$
$F dux = - \frac{1}{4} ρ l π d 2 (\frac{3}{2} C s \overset{\cdot}{d} 2 + d \overset{··}{d}) \cos θ i, Verticalheatingsurface$ $F dux = - \frac{1}{4} ρ l π d 2 (\frac{3}{2} C s \overset{\cdot}{d} 2 + d \overset{··}{d}) \sin θ i, Horizontalheatingsurface$
$F duy = - \frac{1}{4} ρ l π d 2 (\frac{3}{2} C s \overset{\cdot}{d} 2 + d \overset{··}{d}) \sin θ i, Verticalheatingsurface$ $F duy = - \frac{1}{4} ρ l π d 2 (\frac{3}{2} C s \overset{\cdot}{d} 2 + d \overset{··}{d}) \cos θ i, Horizontalheatingsurface$

The drag and shear lift coefficients in Table 4 are determined according to the relationships proposed by Klausner et al.⁶⁵

C D = \frac{2}{3} + [(\frac{12}{Re}) 0.65 + 0.796 0.65] - \frac{1}{0.65}

(19)

C L = 3.877 G_{s}^{- 1 / 2} [\frac{1}{Re 2} + 0.014 G_{s}^{2}] \frac{1}{4}

(20)

where the dimensionless shear rate G_s is estimated as:

G s = (\frac{dU}{dx}) (\frac{r}{Δ U})

. For the force (

F du

) which is due to asymmetrical growth of the bubble, the constant C_s takes the value of

20 / 3

;

\overset{\cdot}{d} (t g)

and

\overset{··}{d} (t g)

represent the first- and second-order derivatives of bubble diameter with respect to bubble growth time (t_g).

While a bubble remains attached to the heated surface, the sum of the parallel and normal forces must satisfy the following conditions: $\sum F x = 0$ and $\sum F y = 0$ . In the flow boiling case, a bubble typically detaches from the nucleation site via sliding and before lifting off the heating surface. For a vertical heating surface, the bubble diameter departing from the cavity site (D_d) and the bubble diameter lifting off the heated surface (D_l) are evaluated by solving $\sum F y = 0$ and $\sum F x = 0$ respectively. For a horizontal heating surface, the bubble diameter departing from the cavity site (D_d) and the bubble diameter lifting off the heated surface (D_l) are nonetheless evaluated by solving $\sum F x = 0$ and $\sum F y = 0$ , respectively. In the pool boiling case, the departure and lift-off points coincide; this simply entails the evaluation of the bubble diameter at the departure point by solving $\sum F x = 0$ for the vertical heating surface and $\sum F y = 0$ for the horizontal heating surface.

Bubble growth time and bubble waiting time

Based on the description of an ebullition cycle during nucleate boiling, the nucleation frequency can be mechanistically determined from

f = \frac{1}{t g + t w}

(21)

where t_g is the growth time of the bubbles and t_w is the waiting time.

The bubble growth time can be determined through a diffusion-controlled bubble growth solution proposed by Zuber⁶⁸

d (t g) = \frac{4 b}{\sqrt{π}} Ja \sqrt{η l t g}

(22)

where

Ja

is the Jakob defined as

Ja = ρ l C p, l Δ T \sup / ρ g h fg

, η is the thermal diffusivity given by

η l = k l / ρ l C p, l

and b is an empirical constant that accounts for the asphericity of the bubble. Once the bubble diameter departing from the cavity site is determined from the force balance model, the bubble growth time can be immediately calculated from equation (22).

The bubble waiting time can be evaluated through the occurrence of transient conduction when a bubble slides or lifts off of which the boundary layer gets disrupted and cold liquid comes in contact with the heated wall. It is estimated as

t w = \frac{1}{π η} [\frac{(Δ T \sup + Δ T sub) C 1 r c}{Δ T \sup + 2 σ T sat / C 2 ρ g h fg r c}] 2

(23)

In equation (23), C₁ and C₂ are $C 1 = (1 + \cos θ) / \sin θ$ and $C 2 = 1 / \sin θ$ , and r_c is the cavity radius

r c = F [\frac{2 σ T sat k l}{ρ g h fg Q w}] \frac{1}{2}

(24)

where

F = [\frac{1}{C 1 C 2}] \frac{1}{2} [\frac{sin 2 θ}{1 + \cos θ}] \frac{1}{2}

(25)

According to Basu et al.,⁶⁹ the factor F indicates the degree of flooding of the available cavity size and the wettability of the surface. If the contact angle θ → 0, all the cavities will be flooded. Alternatively, as θ → 90°, F → 1, all the cavities will not be flooded (i.e. they contain traces of gas).

Alternatively, Podowski et al.⁷⁰ formulated mechanistic expressions for both the waiting time and the growth time. The bubble waiting time can be derived based on the balance of the transient heat transfer inside the wall and fluid near the heated surface

t w = [(- C 5 + \sqrt{C_{5}^{2} - 4 C 4 C 6}) / 2 C 5] 2

(26)

where

C 4 = 2 Q w / R π R = \frac{k w}{\sqrt{π η w}} + \frac{k l}{\sqrt{π η l}} C 5 = (\frac{k w T w}{\sqrt{η w}} + \frac{k l T l}{\sqrt{η l}}) / (\frac{k w}{\sqrt{η w}} + \frac{k l}{\sqrt{η l}}) - T sat - \frac{Q w r c}{R \sqrt{π η l}} - \frac{2 σ}{r c} \frac{T sat (1 / 1 ρ g ρ g - 1 / 1 ρ l ρ l)}{h fg} C 6 = - [(\frac{k w T w}{\sqrt{η w}} + \frac{k l T l}{\sqrt{η l}}) / (\frac{k w}{\sqrt{η w}} + \frac{k l}{\sqrt{η l}}) - T l] r c / \sqrt{π η l}

The bubble growth time is however obtained through the consideration of the critical bubble size and an energy balance during the bubble growth

t g = [(- A 2 + \sqrt{A_{2}^{2} - 4 A 1 A 3}) / 2 A 1] 2

(27)

where

A 1 = Q w / k w A 2 = 2 [(k w T w / \sqrt{η w} + k l T l / \sqrt{η l}) / (k w / \sqrt{η w} + k l / \sqrt{η l}) + 2 Q w \sqrt{t w} / R π - T sat] / \sqrt{π η w} A 3 = d (t g) ρ g h fg / 2 k w

Two-fluid model for nucleate boiling

Away from the heated surface, ensemble-averaged of mass, momentum and energy transport equations are considered for each phase in the Eulerian–Eulerian framework to resolve the two-phase boiling and bubble behaviours in the bulk liquid. The liquid is generally taken as to be the continuum phase while the vapour (bubbles) is treated as the disperse phase. These equations can be derived as

Mass conservation of liquid phase

\frac{\partial ρ l α l}{\partial t} + \nabla . (ρ l α l u l) = Γ \lg

(28)

Mass conservation of vapour phase

\frac{\partial ρ g α g}{\partial t} + \nabla . (ρ g α g u g) = - Γ \lg

(29)

Momentum conservation of liquid phase

\frac{\partial ρ l α l u l}{\partial t} + \nabla . (ρ l α l u l u l) = - α l \nabla P + α l ρ l g + \nabla [α l μ_{l}^{eff} (\nabla u l + (\nabla u l) T)] + (Γ \lg u g - Γ gl u l) + F \lg

(30)

Momentum conservation of vapour phase

\frac{\partial ρ g α g u g}{\partial t} + \nabla . (ρ g α g u g u g) = - α g \nabla P + α g ρ g g + \nabla [α g μ_{g}^{eff} (\nabla u g + (\nabla u g) T)] + (Γ gl u l - Γ \lg u g) + F gl

(31)

Energy conservation of liquid phase

\frac{\partial ρ l α l h l}{\partial t} + \nabla . (ρ l α l u l h l) = \nabla [α l λ_{l}^{eff} (\nabla T l)] + (Γ gl h l - Γ \lg h g)

(32)

Energy conservation of vapour phase

\frac{\partial ρ g α g h g}{\partial t} + \nabla . (ρ g α g u g h g) = \nabla [α g λ_{g}^{eff} (\nabla T g)] + (Γ \lg h l - Γ gl h g)

(33)

The interfacial force

F \lg

in equation (30) is usually formulated through appropriate consideration of different sub-forces affecting the interface between the liquid and vapour phases. For the liquid phase, the interfacial momentum exchange consists of the sum of drag, lift, virtual mass, wall lubrication and turbulent dispersion forces which are summarised in Table 5. It should be noted that

F gl = - F \lg

for the vapour phase. Detail descriptions of these forces can be found in Anglart and Nylund,⁷¹ Lahey and Drew⁷² and Yao and Morel.⁷³

Table 5.

Interfacial momentum forces.

Interfacial forces	Expressions
Drag	$F_{\lg}^{drag} = \frac{1}{8} C D a if ρ l \| u g - u l \| (u g - u l)$
Lift	$F_{\lg}^{lift} = C L α g ρ l (\nabla \times u l) \times (u g - u l)$
Virtual mass	$F_{\lg}^{virtual mass} = C VM α g ρ l (\frac{D u g}{Dt} - \frac{D u l}{Dt})$
Wall lubrication	$F_{\lg}^{wall lubrication} = \frac{α g ρ l [(u g - u l) - (u g - u l) Δ n w] 2}{D s} (C w 1 + C w 2 \frac{D s}{x w}) n w$
Turbulent dispersion	$F_{\lg}^{turbulent dispersion} = - C TD [\frac{1}{8} C D a if ρ l \| u g - u l \|] \frac{μ Tg}{ρ g Sc b} (\frac{\nabla α g}{α g} - \frac{\nabla α l}{α l}), Burnsetal . 82$ $F_{\lg}^{turbulent dispersion} = - C TD ρ l k l \nabla α l, Antaletal . 80$

In Table 5, the drag force is a result of the shear and form drag of the fluid flow. Coefficient C_D appearing in the drag force, normally known as the drag coefficient, can be determined based on the correlations proposed by Ishii and Zuber⁷⁴ for different flow regimes (bubbly, slug or churn-turbulent). As bubbles rise in a liquid, they are subjected to a lateral lift force which acts perpendicular to the direction of relative motion between the two phases. Lopez de Bertodano⁷⁵ and Takagi and Matsumoto⁷⁶ have proposed a value of $0.1$ for the lift coefficient C_L. Drew and Lahey⁷⁷ suggested a value of $0.5$ based on the consideration of an inviscid flow around a sphere. Wang et al.⁷⁸ employed a value of $0.01$ which was deemed to be appropriate for viscous flows. Nevertheless, Tomiyama⁷⁹ developed an Eotvos number dependent correlation that allows negative coefficients for bubble diameter larger than $5.5$ mm in an air–water system; a negative lateral lift force results in large bubbles to be emigrated towards the centre of the flow channel. In contrast to the lateral lift force, wall lubrication force represents another lateral force due to surface tension which is formed to prevent bubbles from attaching on the solid wall. On the basis of numerical experimentation, the wall lubrication constants are: $C w 1 = - 0.01$ and $C w 2 = 0.05$ according to Antal et al.,⁸⁰ which was subsequently revised by Krepper et al.⁸¹ to $C w 1 = - 0.0064$ and $C w 2 = 0.016$ . The virtual mass effect is significant when the gas phase density is much smaller than the liquid phase density. For an inviscid flow around an isolated sphere, the constant $C VM$ is set at a value of 0.5. Burns et al.⁸² derived a model for the turbulence dispersion force based on the consistency of Favre-averaging. The constant $C TD$ is normally set to a value of unity while the turbulent bubble Schmidt number $Sc b$ adopts a value of $0.9$ . Alternatively, values of constant $C TD$ for the model developed by Antal et al.⁸⁰ ranges from $0.1$ to $0.5$ .

On the basis of the eddy viscosity hypothesis, the effective viscosity in the momentum liquid phase transport equation can be taken as the sum of the molecular viscosity and turbulent viscosity of each phase. With the inclusion model for bubble-induced turbulence by Sato et al.,⁸³ the effective viscosity is thus considered as the total of the shear-induced turbulent viscosity and the bubble-induced turbulent viscosity. The effective viscosity of the liquid phase in equation (30) can thus be expressed as

μ_{l}^{eff} = μ l + μ Tl + μ Tb

(34)

In equation (34), the liquid turbulent viscosity is given by

μ Tl = ρ l C μ \frac{k_{l}^{2}}{ɛ l}

(35)

while the extra bubble-induced turbulent viscosity is determined in accordance with

μ Tb = ρ l C μ b α g D s | u g - u l |

(36)

For the vapour phase, the turbulent viscosity can be simply evaluated as

μ Tg = μ Tl \frac{ρ g}{ρ l}

(37)

In order to evaluate the liquid turbulent viscosity in equation (35), the two-phase two-equation k–ɛ turbulence model, which is mere extension of its single-phase counterpart, may be employed to resolve the turbulent flow associated with the continuous liquid phase.⁸⁴ The transport equations for the turbulent kinetic energy of the liquid phase (k_l) and dissipation of the turbulent kinetic energy of the liquid phase ( $ɛ l$ ) can be written as

\frac{\partial ρ l α l k l}{\partial t} + \nabla . (ρ l α l u l k l) = \nabla [α l \frac{μ Tl}{σ k} (\nabla k l)] + α l (P l + G l - ρ l ɛ l)

(38)

\frac{\partial ρ l α l ɛ l}{\partial t} + \nabla . (ρ l α l u l ɛ l) = \nabla [α l \frac{μ Tl}{σ ɛ} (\nabla ɛ l)] + α l \frac{ɛ l}{k l} (C ɛ 1 P l + C ɛ 3 G l - C ɛ 2 ρ l ɛ l)

(39)

The shear production of turbulence P_l is given by

P l = μ Tl \nabla u l . [(\nabla u l + (\nabla u l) T)] - \frac{2}{3} \nabla . u l + (μ Tl \nabla . u l - ρ l u l)

(40)

while the production due to the gravity G_l, which is only applicable for weakly compressible flows, can be expressed in the form

G l = - \frac{μ Tl}{ρ l σ ρ} g . \nabla ρ l

(41)

Appropriate values for $C μ$ , $C μ b$ , $C ɛ 1$ , $C ɛ 2$ , $C ɛ 3$ and $σ ρ$ can be found in Yeoh and Tu.⁸⁴ Nevertheless, it should be noted that Menter’s shear stress transport (SST)⁸⁵ could also be alternatively employed to resolve the two-phase turbulent boiling flow. In the recent works by Cheung et al.⁸⁶ and Frank et al.,⁸⁷ it has been observed that the ‘wall’ peaking behaviour of the void fraction distribution in an isothermal bubbly flow was better predicted through the SST turbulent model. This could be attributed to a more realistic prediction of the turbulent dissipation of the kinetic energy close to the wall through the two-equation k–ω formulation of the liquid phase. Further development should be focused on two-phase flow turbulence modeling in order to better understand or improve existing models.

In the case of subcooled boiling, the interfacial mass transfer rate due to condensation in the bulk liquid in equation (28) can be formulated as

Γ \lg = \frac{h if a if Δ T sub}{h fg}

(42)

where

h if

is the interfacial heat transfer coefficient and

a if

is the interfacial area concentration. The wall nucleation rate can be modelled in a mechanistic manner by considering the total mass of bubbles detaching from the heated surface in accordance with

\overset{\cdot}{m e} = \frac{Q e}{h fg + C p, l Δ T sub}

(43)

For nucleate boiling, the wall nucleation rate is generally treated as a specified boundary condition, which can be treated as the boundary mass flux imposed on the heating surface.

Population balance model

CFD coupled with an appropriate population balance model are increasingly being employed to predict the evolution of the particle (bubble) size distribution in the gas–liquid systems. This is particularly important especially in resolving the interfacial momentum transfer in the two-fluid model where the interfacial forces appearing in the momentum transport equations require the knowledge of the bubble size (Sauter mean diameter) or interfacial area concentration, which is crucial to the overall prediction of the local phase distribution. Some earlier studies have simplified the two-phase problem with a prescribed single bubble size (or so-called mono-dispersed bubble size).^88,89 Such assumption is obviously invalid beyond bubbly flow conditions such as can be found in practical gas–liquid systems. Mounting interests have spurred numerous studies in attempting to better synthesise the behaviour of bubbles and track their dynamical evolution within the two-phase flow domain.

In accordance with Cheung et al.⁹⁰ and Yeoh et al.,⁹¹ the population balance equation can be expressed in an integro-differential form as

\frac{\partial f (x, ξ, t)}{\partial t} + \nabla \cdot (U (x, ξ, t) f (x, ξ, t)) = S (x, ξ, t)

(44)

where

f (x, ξ, t)

is the particle (bubble) number density distribution per unit mixture and particle (bubble) volume,

U (x, ξ, t)

is velocity vector in external space dependent on the external variables x for a given time t and the internal space ξ whose components could be the characteristic dimensions such as volume, mass, etc. The term

S (x, ξ, t)

contains the particle (bubble) source/sink rates per unit mixture volume due to the particle (bubble) interactions such as coalescence, break-up as well as phase change.

In practice, the transport of particle (bubble) number density needs to be solved within the framework of CFD. Following the formulation proposed by Prince and Blanch,⁹² Luo⁹³ and Luo and Svendsen,⁹⁴ the continuous size range of particles represented by the particle (bubble) number density distribution can be discretised into a series number of discrete size classes. By taking the volume V_p as the independent internal coordinate, the population balance equation for the discrete particle (bubble) number density defined by

N i (x, t) = \int_{v i}^{v i + 1} f (x, ξ, t) d v

(45)

can be written in accordance with

\frac{\partial N i (x, t)}{\partial t} + \nabla \cdot (V (x, t) N i (x, t)) = S N i (x, t)

(46)

The population balance equation for the average particle (bubble) number density defined by

N (x, t) = \int_{0}^{v} f (x, ξ, t) d v

(47)

can nonetheless be written in a similar form as

\frac{\partial N (x, t)}{\partial t} + \nabla \cdot (V (x, t) N (x, t)) = S N (x, t)

(48)

Population changes in equations (46) and (48) are accommodated through mechanistic models being formulated for the source terms within these transport equations for particular gas–liquid systems.

Considerable attention has been focused on the description of the temporal and spatial evolution of the geometrical structure caused by the effects of coalescence and break-up through the interactions among bubbles as well as between bubbles and turbulent eddies in turbulent flows. For bubbly flow, appropriate mechanistic models for the source terms in equations (46) and (48) have been established: (i) Coalescence through random collision driven by turbulent eddies; (ii) Coalescence due to the acceleration of the following bubble in the wake of the preceding bubble; and (iii) Break-up due to the impact of turbulent eddies. Schematic illustrations of these mechanisms are illustrated in Figure 4. Table 6 tabulates three mechanistic models that have been considered for the consideration of the average particle (bubble) number density equation while Table 7 lists some of the more popular mechanistic models that have been applied to characterise the coalescence and break-up of bubbles for the discrete particle (bubble) number density equation. More details can be found in Yeoh et al.⁹¹ In Table 6, $Φ RC$ and $Φ TI$ are the coalescence and break-up mechanisms of bubbles due to random collision and turbulent impact for the average particle (bubble) number density equation. In Table 7, $a (M i, M j)$ and $β (f BV, 1)$ are also the coalescence and break-up mechanisms of bubbles due to random collision and turbulent impact but for the discrete particle (bubble) number density equation with $C 1'$ , $C 2'$ , $C 3'$ , $C 4'$ and $C 5'$ being the adjustable constants for the respective models.

Figure 4.

Bubbles undergoing coalescence and break-up in gas–liquid systems.

Table 6.

Coalescence and break-up mechanisms for averaged particle (bubble) number density.

Contributors	Mechanistic models
Wu et al.⁹⁵	$Φ RC = - C RC 1 \frac{α_{g}^{2} ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133} α_{g, \max}^{1 / 133} (α_{g, \max}^{1 / 133} - α_{g}^{1 / 133})} [1 - \exp (- C RC 2 \frac{α_{g, \max}^{1 / 133} α_{g}^{1 / 133}}{α_{g, \max}^{1 / 133} - α_{g}^{1 / 133}})]$ $Φ TI = - C TI \frac{α g ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133}} (1 - \frac{We cr}{We}) \exp (- \frac{We cr}{We})$ $C RC 1 = 0.021, C RC 2 = 3.0, C TI = 0.0945, We cr = 2.3, α g, \max = 0.8$
Hibiki and Ishii⁹⁶	$Φ RC = - C RC 1 \frac{α_{g}^{2} ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133} (α g, \max - α g)} \exp (- C RC 2 \frac{ρ_{l}^{1 / 122} ɛ_{l}^{1 / 133} D_{s}^{5 / 566}}{σ 1 / 122})$ $Φ TI = - C TI \frac{α g (1 - α g) ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133} (α g, \max - α g)} \exp (- C TI 2 \frac{σ}{ρ l ɛ_{l}^{2 / 233} D_{s}^{5 / 533}})$ $C RC 1 = 0.03, C RC 2 = 1.29, C TI 1 = 0.03, C TI 2 = 1.37, α g, \max = 0.8$
Yao and Morel⁷⁰	$Φ RC = - C RC 1 \frac{α_{g}^{2} ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133} (α g, \max - α g)} \frac{\exp (- C RC 2 (We / We cr) 1 / 122)}{α_{g, \max}^{1 / 133} (α_{g, \max}^{1 / 133} - α_{g}^{1 / 133}) + C RC 3 α g (We / We cr) 1 / 122}$ $Φ TI = - C TI \frac{α g (1 - α g) ɛ_{l}^{1 / 133}}{D_{s}^{11 / 1133}} \frac{\exp (- We cr / We)}{1 + C TI 2 (1 - α g) (We / We cr) 1 / 122}$ $C RC 1 = 2.86, C RC 2 = 1.017, C TI 1 = 1.922, C TI 2 = 1.6, We cr = 1.42, α g, \max = 0.52$

Table 7.

Coalescence and break-up mechanisms for discrete particle (bubble) number density.

Contributors	Mechanistic models
Coulaloglou and Tavlarides⁹⁷	$a (M i, M j) = C 1' (d i + d j) 2 (d_{i}^{2 / 233} + d_{j}^{2 / 233}) 1 / 122 ɛ_{l}^{1 / 133} \exp [- C 2' \times \frac{μ l ρ l ɛ l}{σ 2} (\frac{d i d j}{d i + d j}) 4]$
Prince and Blanch⁹²	$a (M i, M j) = C 3' (d i + d j) 2 (d_{i}^{2 / 233} + d_{j}^{2 / 233}) 1 / 122 ɛ_{l}^{1 / 133} \exp (\frac{t ij}{τ ij})$
Lehr et al.⁹⁸	$a (M i, M j) = C 4' (d i + d j) 2 (d_{i}^{2 / 233} + d_{j}^{2 / 233}) 1 / 122 ɛ_{l}^{1 / 133} \exp [- (\frac{α_{g, \max}^{1 / 133}}{α_{g}^{1 / 133}} - 1) 2] \min (\frac{u cr}{u'}, 1)$
Luo and Svendsen⁹⁴	$β (f BV, 1) = \frac{2 \int_{ξ \min}^{1} \frac{(1 + ξ) 2}{ξ 11 / 3} exp (- \frac{12 C f σ}{β ρ l ɛ_{l}^{2 / 233} d_{i}^{5 / 533} ξ 1 / 133}) d ξ}{\int_{0}^{1} \int_{ξ \min}^{1} \frac{(1 + ξ) 2}{ξ 11 / 3} exp (- \frac{12 C f σ}{β ρ l ɛ_{l}^{2 / 233} d_{i}^{5 / 533} ξ 1 / 133}) d ξ d f BV}$
Wang et al.⁹⁹	$β (f BV, 1) = \frac{\int_{λ \min}^{d i} \frac{(λ + d) 2}{λ 11 / 1133} \int_{0}^{\infty} \frac{1}{f BV, \max - f BV, \min} \frac{1}{\bar{e} (λ)} \exp (- \frac{e (λ)}{\bar{e} (λ)}) d e (λ) d λ}{\int_{0}^{1} \int_{λ \min}^{d i} \frac{(λ + d) 2}{λ 11 / 1133} \int_{0}^{\infty} \frac{1}{f BV, \max - f BV, \min} \frac{1}{\bar{e} (λ)} \exp (- \frac{e (λ)}{\bar{e} (λ)}) d e (λ) d λ d f BV}$
Martinez-Bazan et al.^100,101	$β (f BV, 1) = \frac{[\frac{1}{2} ρ l C 5' (ɛ l d j) 2 / 3 - 6 \frac{σ}{d i}] \times [\frac{1}{2} ρ l C 5' (ɛ l (d_{i}^{3} - d_{j}^{3}) 1 / 3) 2 / 3 - 6 \frac{σ}{d i}]}{\int_{0}^{1} [\frac{1}{2} ρ l C 5' (ɛ l d j) 2 / 3 - 6 \frac{σ}{d i}] \times [\frac{1}{2} ρ l C 5' (ɛ l (d_{i}^{3} - d_{j}^{3}) 1 / 3) 2 / 3 - 6 \frac{σ}{d i}] d f BV}$

Beyond bubbly flow, the requirement of accommodating other bubble sizes and shapes involves the intra-group and inter-group interactions for spherical and cap/slug bubbles. The probable mechanisms of these interactions that need to be considered are: (i) Coalescence due to random collisions driven by liquid turbulence; (ii) Coalescence due to wake entrainment; (iii) Break-up due to the impact of turbulent eddies; (iv) Shearing-off of small bubbles from cap/slug bubbles; and (v) Break-up of large cap bubbles due to flow instability on the bubble surface.

For intra-group mechanisms for spherical bubbles, the usual coalescence and break-up processes due to random collisions and turbulent impact can be adopted. The intra-group mechanisms for cap/slug bubbles results in a more complex consideration where the bubble coalescence is due to random collisions and wake entrainment while the bubble break-up is attributed by turbulent impact, shearing-off and surface instability. Referring to Figure 5, the wake entrainment for cap/slug bubbles predominantly governs the cap/slug bubble number, which significantly affects the flow structure and intensiveness of inter-group interactions. Here, when bubbles enter the critical wake region of a preceding bubble, they can accelerate and coalesce with this preceding bubble according to the mechanisms based upon the critical distance of a probability approach such as proposed by Sun et al.¹⁰² Nevertheless, the cap/slug bubble disintegration due to surface instability can be significantly enhanced by the high turbulent intensity and active eddy-bubble interaction in the wake region of the slug bubbles. This occurs when the volume of the resulting bubble from bubble coalescence exceeds the maximum bubble stable limit; it becomes unstable and disintegrates. The mechanism for shearing-off effect is reflected by the fact that the bubbles shearing-off from the cap/slug bubbles contribute as a primary source for the spherical bubbles and ultimately the interfacial area concentration. More details on the mechanisms for cap/slug bubbles can be found in Hibiki and Ishii,¹⁰³ Fu and Ishii,¹⁰⁴ Hibiki and Ishii¹⁰⁵ and Cheng et al.¹⁰⁶ In addition to the intra-group interactions, it has been demonstrated through experiments performed by Hibiki and Ishii¹⁰³ that the inter-group interactions also contribute significantly to the transport of the interfacial states. Similar mechanisms from the above considerations for intra-group interactions may be derived accordingly for the inter-group interactions.

Figure 5.

Schematic drawings of the mechanisms associated with the two-group bubble interactions.¹⁰⁵

Concluding remarks

The boiling behaviour of cryogenic liquids is considered to be different from the boiling behaviour of ordinary liquids. In comparison to an ordinary liquid such as water, cryogenic liquids can only remain in their liquid form within a much narrower atmospheric liquid range. A review on the use of empirical correlations and theoretical development of mechanistic models to handle nucleate boiling of cryogenic liquids is presented. In spite of the prospect of adopting a mechanistic approach in modelling the heat transfer mechanisms occurring during the complex boiling process, which can be readily employed for both the prediction of the wall heat flux as well as the partitioning of the wall heat flux between the liquid and gas phases, limitations still persist due to the inadequacy of empirical correlations in fully achieving the necessary closure. The proposal of the fractal and force balance models to mechanistically determine the active nucleation site density and bubble diameters lifting off from the heated wall can in some way elevate the deficiencies that are prevalent in the use of empirical correlations. In addition, the necessary closure for the inter-phase mass, momentum and energy exchanges that are required to be accounted in the two-fluid model requires constitutive relations which can be realised through the use of mechanistic population balance model in predicting the local changes of the Sauter mean diameter or interfacial area concentration. Nevertheless, most mechanistic models for bubble coalescence and break-up that have been developed thus far in the population balance model are established mainly for the bubbly flow regime. Extensive work is still required to identify proper bubble mechanisms and develop models for flow regimes that may be filled with complex shape bubbles such as slugs and columns as the wall heat flux increases to CHF.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Appendix

References

NASA.org. About the space launch system and multi-purpose crew vehicle, http://www.nasa.gov/exploration/new_space_enterprise/sls_mpcv/index.html.

Evans

. The Large Hadron Collider. Ann Rev Nucl Part Sci 2011; 61: 435–466.

ten Kate

HHJ

. The ATLAS superconducting magnet system at the Large Hadron Collider. Phys C: Supercond 2008; 468: 2137–2142.

Yeoh

. Modelling subcooled boiling flows, New York: Nova Science Publishers, 2009.

Warrier

Dhir

. Heat transfer and wall heat flux partitioning during subcooled flow nucleate boiling. ASME J Heat Transfer 2006; 128: 1243–1256.

Bergles

Rohsenow

. The determination of forced convection, surface boiling heat transfer. ASME J Heat Transfer 1964; 86: 365–372.

Chen

. A correlation for boiling heat transfer to saturated fluids in convective flows. Ind Eng Chem Process Des Dev 1966; 5: 322–329.

McAdams

Kennel

Minden

. Heat transfer at high rates to water with surface boiling. Ind Eng Chem 1949; 41: 1945–1953.

Kandlikar

. Heat transfer characteristics in partial boiling, fully developed boiling, and significant void flow regions of subcooled flow boiling. ASME J Heat Transfer 1998; 120: 395–401.

10.

Chatoorgoon

Dimmick

Carver

. Application of generation and condensation models to predict subcooled boiling void at low pressure. Nucl Technol 1992; 98: 366–378.

11.

Zeitoun O. Subcooled flow boiling and condensation. PhD Thesis, McMaster University, Canada, 1994.

12.

Kurul N and Podowski MZ. Multi-dimensional effects in forced convection sub-cooled boiling. In: Proceedings of the 9th heat transfer conference, Jerusalem, Israel, Hemisphere Publishing Corporation, 1990, pp. 21–26.

13.

Podowski

. Towards mechanistic modelling of boiling heat transfer. Nucl Eng Technol 2012; 44: 889–896.

14.

Yeoh

. On numerical modelling of low-pressure subcooled boiling flows. Int J Heat Mass Transfer 2002; 45: 1197–1209.

15.

Wang

Huang

. Numerical investigation of boiling flow of nitrogen in a vertical tube using the two-fluid model. Appl Therm Eng 2006; 26: 2425–2432.

16.

Wei

Huang

. Numerical and experimental investigation of heat transfer on heating surface during subcooled boiling flow of liquid nitrogen. Int J Heat Mass Transfer 2006; 52: 1510–1516.

17.

Yeoh

. Population balance modeling for bubbly flows with heat and mass transfer. Chem Eng Sci 2004; 59: 3125–3139.

18.

Yeoh

. Thermal hydraulic modeling of bubbly flows with heat and mass transfer. AIChE J 2005; 51: 8–27.

19.

Benjamin

Balakrishnan

. Nucleate pool boiling heat transfer of pure liquids at low to moderate heat fluxes. Int J Heat Mass Transfer 1996; 39: 2495–2504.

20.

Gaertner

Westwater

. Population of active sites in nucleate boiling heat transfer. Chem Eng Prog Stmp Ser 1960; 56: 39–48.

21.

Hsu

Y-Y

Graham

. Transport processes in boiling and two-phase systems, La Grange Park, IL: American Nuclear Society, 1986.

22.

Westwater

Kirby

. Bubble and vapor behavior on a heated horizontal plate during pool boiling near burnout. Chem Eng Prog Symp Ser 1965; 61: 238–248.

23.

Johov K. Nucleations number during steam production, aerodynamics and heat transfer in the working elements of the power facilities, Leningrad. In: Proceedings of CKTI, vol. 91, 1969, pp. 131–135.

24.

Mikic

Rohsenow

. A new correlation of pool boiling data including the effect of heating surface characteristics. ASME J Heat Transfer 1969; 91: 245–256.

25.

Bier K, Gorenflo D, Salem M, et al. Pool boiling heat transfer and size of active nucleation centers for horizontal plates with different roughness. In: Proceedings of 6th international heat transfer conference, Hemisphere, Washington, USA, 1978, pp. 151–156.

26.

Cornwell K and Brown E. Boiling surface topography. In: Proceedings of 6th international heat transfer conference, Hemisphere, Washington, USA, 1978, pp. 157–161.

27.

Lemmert

Chwala

. Influence of flow velocity on surface boiling heat transfer coefficient, New York and Washington: Academic Press and Hemisphere, 1977.

28.

Končar

Kljenak

Mavko

. Modeling of local two-phase flow parameters in upwards subcooled flow boiling at low pressure. Int J Heat Mass Transfer 2004; 47: 1499–1513.

29.

Krepper

Rzehak

. CFD for subcooled flow boiling: Simulation of DEBORA experiments. Nucl Eng Des 2011; 241: 3851–3866.

30.

Krepper

Rzehak

Lifante

. CFD for subcooled flow boiling: Coupling wall boiling and population balance models. Nucl Eng Des 2013; 255: 330–346.

31.

Kolev

. Uniqueness of the elementary physics driving heterogeneous nucleate boiling and flashing. Nucl Eng Technol 2006; 38: 33–42.

32.

Bartolomej

Chanturiya

. Experimental study of true void fraction when boiling subcooled water in vertical tubes. Therm Eng 1967; 14: 123–128.

33.

Kirichenko lu

Dolgoy

Levchenko

. Study of the boiling of cryogenic liquids. Heat Transfer Soviet Res 1976; 8: 63–72.

34.

Kocamustafaogullari

Ishii

. Interfacial area and nucleation site density in boiling systems. Int J Heat Mass Transfer 1983; 26: 1377–1387.

35.

Wang

Dhir

. Effect of surface wettability on active nucleation site density during pool boiling of water on a vertical surface. ASME J Heat Transfer 1993; 115: 659–669.

36.

Hibiki

Ishii

. Active nucleation site density in boiling systems. Int J Heat Mass Transfer 2003; 46: 2587–2601.

37.

Basu

Warrier

Dhir

. Wall heat flux partitioning during subcooled flow boiling: Part 1 – model development. ASME J Heat Transfer 2005; 127: 131–140.

38.

Fritz

. Maximum volume of vapor bubbles. Physik Zeitschr 1935; 36: 379–384.

39.

Zuber N. Hydrodynamic aspects of boiling heat transfer. PhD Thesis, Research Laboratory, Los Angeles and Ramo-Wooldridge Corporation, University of California, USA, 1959.

40.

Han

C-Y

Griffith

. The mechanism of heat transfer in nucleate pool boiling, Cambridge, MA: MIT Division of Sponsored Research, 1962.

41.

Cole

Shulman

. Bubble departure diameters at subatmospheric pressures. Chem Eng Prog Symp Ser 1967; 62: 6–16.

42.

Cole

. Bubble frequencies and departure volumes at subatmospheric pressures. AIChE J 1967; 13: 779–783.

43.

Cole

Rohsenow

. Correlation of bubble departure diameters for boiling of saturated liquids. Chem Eng Prog Symp Ser 1969; 65: 211–213.

44.

Kirichenko lu

Slobozhanin

Shcherbakova

. Analysis of quasi-static conditions of boiling onset and bubble departure. Cryogenics 1983; 23: 110–112.

45.

Bald

. Cryogenic heat transfer research at Oxford, Part 1 – nucleate pool boiling. Cryogenics 1973; 13: 457–469.

46.

Farajisarir D. Growth and collapse of vapour bubbles in convective subcooled boiling of water. Masters Thesis, University of British Columbia, Canada.

47.

Tolubinsky V and Kostanchuk D. Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling. In: Proceedings 4th international heat transfer conference, 5 Paris, France, 1970, Paper B-2.8.

48.

Ünal

. Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/m². Int J Heat Mass Transfer 1976; 19: 643–649.

49.

Rohsenow WM. Heat transfer with evaporation. Proceedings of heat transfer—A symposium, University of Michigan, USA, 1952.

50.

Levenspiel

. Collapse of steam bubbles in water. Ind Eng Chem 1959; 51: 787–790.

51.

Morel C, Mimouni S, Lavieville M, et al. R113 boiling bubbly flow in an annular geometry simulated with the NEPTUNE code. In: Proceedings of 11th international topical meeting nuclear reactor thermalhydraulics (NURETH-11), Avignon, France, 2005.

52.

Cole

. A Photographic study of pool boiling in the region of the critical heat flux. AIChE J 1960; 6: 533–538.

53.

McFadden

Grassmann

. The relation between bubble frequency and diameter during nucleate pool boiling. Int J Heat Mass Transfer 1962; 5: 169–173.

54.

Jin

Zhang

Tang

. Observation and analysis of the detachment frequency of coalesced bubbles in pool boiling liquid nitrogen. Cryogenics 2011; 51: 516–520.

55.

Hatton A and Hall I. Photographic study of boiling on prepared surfaces. In: Proceedings of 3rd international heat transfer conference, Chicago, USA, 1966, pp. 24–37.

56.

Ivey

. Relationships between bubble frequency, departure diameter and rise velocity in nucleate boiling. Int J Heat Mass Transfer 1967; 10: 1023–1040.

57.

Stephan

Green

. Heat transfer in condensation and boiling, Berlin: Springer, 1992.

58.

Malenkov

. The frequency of vapor-bubble separation as a function of bubble size. Fluid Mech Sov Res 1972; 1: 36–42.

59.

Yeoh

Cheung

SCP

. Modeling vertical subcooled boiling flows at low pressures. J ASTM Int 2011; 8: 1–20.

60.

Yeoh

Vahaji

Cheung

SCP

. Modeling subcooled flow boiling in vertical channels at low pressures – Part 2: Evaulation of mechanistic approach. Int J Heat Mass Transfer 2014; 75: 754–768.

61.

Xiao

. A fractal analysis of subcooled flow boiling heat transfer. Int J Multiphase Flow 2007; 33: 1126–1139.

62.

Hsu

. On the size range of active nucleation cavities on a heating surface. ASME J Heat Transfer 1962; 84: 207–213.

63.

Cheng

. Fractal models for the effective thermal conductivity of bidispersed porous media. J Thermophys Heat Transfer 2002; 16: 22–29.

64.

Yeoh

. A unified model considering force balances for departing vapour bubbles and population balance in subcooled boiling flow. Nucl Eng Des 2005; 235: 1251–1265.

65.

Klausner

Mei

Bernhard

. Vapor bubble departure in forced convection boiling. Int J Heat Mass Transfer 1993; 36: 651–662.

66.

Zeng

Klausner

Mei

. A unified model for the prediction of bubble detachment diameters in boiling systems – I. Pool boiling. Int J Heat Mass Transfer 1993; 36: 2261–2270.

67.

Zeng

Klausner

Mei

. A unified model for the prediction of bubble detachment diameters in boiling systems – II. Flow boiling. Int J Heat Mass Transfer 1993; 36: 2271–2279.

68.

Zuber

. The dynamics of vapor bubbles in nonuniform temperature fields. Int J Heat Mass Transfer 1961; 2: 83–98.

69.

Basu

Warrier

Dhir

. Onset of nucleate boiling and active nucleation site density during subcooled flow boiling. ASME J Heat Transfer 2002; 124: 717–728.

70.

Podowski RM, Drew DA, Lahey RT Jr, et al. A mechanistic model of the ebullition cycle in forced convection subcooled boiling. In: Proceedings of 8th international topical meeting nuclear reactor thermal-hydraulics (NURETH-8), 1997, pp. 1535–1542.

71.

Anglart

Nylund

. CFD application to prediction of void distribution in two-phase bubbly flows in rod bundles. Nucl Sci Eng 1996; 163: 81–98.

72.

Lahey

Jr Drew

. The analysis of two-phase flow and heat transfer using multidimensional, four field, two-fluid model. Nucl Eng Des 2001; 204: 29–44.

73.

Yao

Morel

. Volumetric interfacial area prediction in upward bubbly two-phase flow. Int J Heat Mass Transfer 2004; 47: 307–328.

74.

Ishii

Zuber

. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J 1979; 25: 843–855.

75.

Lopez de Bertodano M. Turbulent bubbly two-phase flow in a triangular duct. PhD Thesis, Rensselaer Polytechnic Institute, UK, 1992.

76.

Takagi S and Matsumoto Y. Numerical study on the forces acting on a bubble and particle. In: Proceedings of 3rd international conference multiphase flow, Lyon, France, 1998.

77.

Drew

Lahey

Jr . Application of general constitutive principles to the derivation of multi-dimensional two-phase flow equation. Int J Multiphase Flow 1979; 5: 243–264.

78.

Wang

Lee

Lahey

Jr . 3-D Turbulence structure and phase distribution measurements in bubbly two-phase flows. Int J Multiphase Flow 1987; 13: 327–343.

79.

Tomiyama A. Struggle with computational bubble dynamics. In: Proceedings of 3rd international conference on multiphase flow, Lyon, France, 1998.

80.

Antal

Lahey

Jr Flaherty

. Analysis of phase distribution and turbulence in dispersed particle/liquid flows. Chem Eng Commun 1991; 174: 85–113.

81.

Krepper

Lucas

Prasser

. On the modeling of bubbly flow in vertical pipes. Nucl Eng Des 2005; 235: 597–611.

82.

Burns AD, Frank T, Hamill I, et al. The Favre averaged drag model for turbulent dispersion in Eulerian multiphase flow. In: Proceedings of 5th international multiphase flow, Yokohama, Japan, 2004.

83.

Sato

Sadatomi

Sekoguchi

. Momentum and heat transfer in two-phase bubbly flow – I. Int J Multiphase Flow 1981; 7: 167–178.

84.

Yeoh

. Computational techniques for multiphase flows, Oxford, UK: Butterworth-Heinemann, Elsevier, 2005.

85.

Menter

. Two-equation eddy viscosity turbulence models for engineering applications. AIAA J 1994; 32: 1598–1605.

86.

Cheung

SCP

Yeoh

. On the modeling of population balance in isothermal vertical bubbly flows – Average bubble number density approach. Chem Eng Proc 2007; 46: 742–756.

87.

Frank T, Shi J and Burns AD. Validation of Eulerian multiphase flow models for nuclear safety application. In: Proceedings of 3rd symposium on two-phase modeling & experimentation, Pisa, Italy, 2004.

88.

Krishna

van Baten

Urseanu

. Three-phase Eulerian simulations of bubble column reactors operating in the churn-turbulent regime: A scale up strategy. Chem Eng Sci 2000; 55: 3275–3286.

89.

Pan Y and Dudukovic MP. CFD simulations of a bubble column-2D versus 3D. In: 6th world congress of chemical engineering, Melbourne, Australia, Nauka/Interperiodica, 2001.

90.

Cheung

SCP

Yeoh

. Review of population balance modelling for isothermal bubbly flows. J Comput Multiphase Flows 2009; 1: 161–199.

91.

Yeoh

Cheung

SCP

. Multiphase flow analysis using population balance modeling: bubbles, drops and particles, Oxford, UK: Butterworth-Heinemann, Elsevier, 2013.

92.

Prince

Blanch

. Bubble coalescence and break-up in air-sparged bubble column. AIChE J 1990; 36: 1485–1499.

93.

Luo H. Coalescence, break-up and liquid recirculation in bubble column reactors. PhD Thesis, Norwegian Institute of Technology, Norway, 1993.

94.

Luo

Svendsen

. Theoretical model for drop and bubble break-up in turbulent dispersions. AIChE J 1996; 42: 1225–1233.

95.

Kim

Ishii

. One-group interfacial area transport in vertical bubbly flow. Int J Heat Mass Transfer 1998; 41: 1103–1112.

96.

Hibiki

Ishii

. Development of one-group interfacial area transport equation in bubbly flow systems. Int J Heat Mass Transfer 2002; 45: 2351–2372.

97.

Coulaloglou

Tavlarides

. Description of interaction processes in agitated liquid-liquid dispersions. Chem Eng Sci 1977; 32: 1289–1297.

98.

Lehr

Millies

Mewes

. Bubble-size distributions and flow fields in bubble columns. AIChE J 2002; 48: 2426–2443.

99.

Wang

Jin

. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow. Chem Eng Sci 2003; 58: 4629–4637.

100.

Martinez-Bazan

Montanes

Lasheras

. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 1. Breakup frequency. J Fluid Mech 1999; 401: 157–182.

101.

Martinez-Bazan

Montanes

Lasheras

. On the breakup of an air bubble injected into a fully developed turbulent flow. Part 2. Size PDF of the resulting daughter bubbles. J Fluid Mech 1999; 401: 183–207.

102.

Sun

Kim

Ishii

. Model evaluation of two-group interfacial area transport equation for confined upward flow. Nucl Eng Des 2004; 230: 27–47.

103.

Hibiki

Ishii

. Two-group interfacial area transport equations at bubbly-to-slug flow transition. Nucl Eng Des 2000; 202: 39–76.

104.

Ishii

. Two-group interfacial area transport in vertical air-water flow. II. Model Evaluation. Nucl Eng Des 2003; 219: 169–190.

105.

Hibiki

Ishii

. Interfacial area transport equations for gas-liquid flow. J Comput Multiphase Flows 2009; 1: 1–22.

106.

Cheung

SCP

Yeoh

. Classification of bubbles in vertical gas-liquid flow: Part 2 – A model evaluation. Int J Multiphase Flow 2012; 39: 135–147.