Sage Journals: Discover world-class research

Abstract

In order to accurately predict the thermal hydraulic of two-phase gas–liquid flows with heat and mass transfer, special numerical considerations are required to capture the underlying physics: characteristics of the heat transfer and bubble dynamics taking place near the heated wall and the evolution of the bubble size distribution caused by the coalescence, break-up, and condensation processes in the bulk subcooled liquid. The evolution of the bubble size distribution is largely driven by the bubble coalescence and break-up mechanisms. In this paper, a numerical assessment on the performance of six different bubble coalescence and break-up kernels is carried out to investigate the bubble size distribution and its impact on local hydrodynamics. The resultant bubble size distributions are compared to achieve a better insight of the prediction mechanisms. Also, the void fraction, bubble Sauter mean diameter, and interfacial area concentration profiles are compared against the experimental data to ensure the validity of the models applied.

Keywords

Population balance coalescence break-up multiphase heat and mass transfer subcooled boiling flow wall heat partitioning

Introduction

Two-phase gas–liquid flows with heat and mass transfer, such as subcooled boiling flows in heated channels, are prevalent in various industrial applications. Subcooled boiling flow can be characterized by the thermodynamic non-equilibrium that persists between the liquid and vapor phases inside the two-phase flow. A high-temperature two-phase region exists near the heated wall, while a low-temperature single-phase liquid generally occurs away from the heated surface. Heterogeneous bubble nucleation ensues within the active nucleation sites on the heated surface when the surface temperature exceeds the saturated liquid temperature at local pressure. At the onset of nucleate boiling (ONB), boiling occurs and bubbles remain attached to the heater surface. As the bulk temperature liquid temperature increases, the bubbles grow larger and detach from the heater surface. The void fraction increases sharply designated at the point of net vapor generation (NVG) – a low void fraction region upstream is followed by another region of which the void fraction increases significantly downstream.

In the past decades, modeling subcooled boiling flow based on the computational fluid dynamic techniques have been largely focused on capturing the fundamental considerations of: (i) heat and mass transfer in terms of heat flux partitioning during subcooled boiling flow at the heated wall and (ii) and the two-phase vapor bubble behaviors and its size evolution due to bubble interaction in the bulk subcooled flow away from the heated wall.

For the modeling of heat transfer at the heated wall, empirical correlations for determining the bubble departure diameter, bubble frequency, and active nucleation site density have been heavily utilized to achieve the necessary closure for the wall heat flux partitioning model.^1–6A detailed assessment of the extent of applicability of these correlations has been performed in Cheung et al.⁷ of which their limitations to certain flow conditions were demonstrated. Nevertheless, a further assessment in Vahaji et al.⁸ and Yeoh et al.⁹ on the mechanistically developed wall heat flux partitioning model to circumvent the use of empirical correlations has shown the applicability of such a model to a wider range of heating and flow conditions.

For the modeling of two-phase and bubble behaviors in the bulk flow, majority of studies have thus far concentrated on the study of coalescence and break-up kernels in two-phase isothermal flows.^10,11 Recently, comparative analysis of different coalescence and break-up kernels has been performed by Deju et al.¹² for two-phase isothermal flows in a large bubble column. Nevertheless, such investigations on subcooled boiling flows remain elusive.

The major objectives of the paper are thus twofold: (i) assessing the applicability of existing bubble coalescence and break-up kernels in the framework of population balance modeling for subcooled boiling flows and (ii) utilize the mechanistic wall heat flux partitioning model with specific emphasis at elevated pressures for subcooled boiling flows. A numerical assessment on the performance of six different bubble coalescence and break-up kernels has been performed to investigate the bubble size distribution and its impact on local hydrodynamics based on different heating and flow conditions. For the mechanistic break-up kernels, two widely adopted models with different predictions for daughter size distribution (DSD) proposed by Luo and Svendsen¹³ and Wang et al.¹⁴ are assessed. These break-up kernels are then coupled with three different mechanistic coalescence kernels by Coulaloglou and Tavlarides,¹⁵ Prince and Blanch,¹⁶ and Lehr et al.¹⁷ to form the six different combinations as presented in Table 1. Resulting bubble size distributions in the form of bubble Sauter mean diameter are compared to provide insights on the predictive capability of the bubble coalescence and break-up mechanisms. In addition, void fraction profiles are compared against the experimental data of Ozar et al.¹⁸ to ensure the validity of the model predictions of subcooled boiling flows at elevated pressures.

Mathematical formulation

Two-fluid model

For two-phase subcooled boiling flows, the ensemble-averaged mass, momentum, and energy transport equations for continuous and dispersed phases are modeled using the Eulerian modeling framework. Considering the liquid ( $α_{l}$ ) as continuous phase and the bubbles ( $α_{g}$ ) as disperse phase, numerical simulations are carried out via the two-fluid model Eulerian–Eulerian approach.¹⁹

Mass transport of continuous phase

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

(1)

Mass transport of disperse phase

\partial (ρ_{g} α_{g}) \partial t + \nabla . (ρ_{g} α_{g} u_{g}) = Γ_{gl}

(2)

Momentum transport of continuous phase

\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}

(3)

Momentum transport of disperse phase

\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}

(4)

Energy transport of continuous phase

\partial (ρ_{l} α_{l} H_{l}) \partial t + \nabla . (ρ_{l} α_{l} u_{l} H_{l}) = \nabla . [α_{l} λ_{l} \nabla T_{l} + μ_{Tl} \Pr_{Tl} \nabla H_{l}] + (Γ_{\lg} H_{g} - Γ_{gl} H_{l})

(5)

Energy transport of disperse phase

\partial (ρ_{g} α_{g} H_{g}) \partial t + \nabla . (ρ_{g} α_{g} u_{g} H_{g}) = \nabla . [α_{g} λ_{g} \nabla T_{g} + μ_{Tg} \Pr_{Tg} \nabla H_{g}] + (Γ_{gl} H_{l} - Γ_{\lg} H_{g})

(6)

From the above transport equations, ρ is the density, $u$ is the velocity vector, P is the pressure, $g$ is the gravitational vector, μ is the viscosity, H is the enthalpy, T is the temperature, and λ is the thermal conductivity. The source term $Γ_{\lg}$ represents the mass transfer rate due to condensation in the bulk subcooled liquid, which can be expressed by

Γ_{\lg} = {ha}_{if} (T_{sat} - T_{l}) h_{\lg}

(7)

where h is the inter-phase heat transfer coefficient determined from the Nusselt number correlation by Ranz and Marshall,²⁰

a_{if}

is the interfacial area between phases per unit volume,

T_{sat}

is the saturation temperature and

h_{\lg}

is the latent heat. Note that

Γ_{gl} = - Γ_{\lg}

. The wall vapor generation rate, which is accounted as a boundary condition for the mass transport of the disperse phase, can be evaluated in accordance with

Γ_{wg} = Q_{e} h_{\lg}

(8)

where Q_e is the heat transfer due to evaporation, C_p is the specific heat of constant pressure and

T_{sub}

is the local subcooled temperature.

Interfacial transfer terms in the momentum and energy transport equations ( $Γ_{\lg}$ and $F_{\lg}$ ) represent the transfer terms from the continuous phase to the disperse phase. The total interfacial force $F_{\lg}$ considered in the present study includes the drag, lift, wall lubrication, and turbulent dispersion forces:

F_{\lg} = F_{\lg}^{drag} + F_{\lg}^{lift} + F_{\lg}^{wall lubrication} + F_{\lg}^{turbulent dispersion}

(9)

Note that $F_{gl} = - F_{\lg}$ . More detail descriptions of these forces can be found in Anglart and Nylund.²¹

By the inclusion of Sato et al.²² model for bubble-induced turbulence, the effective viscosity is considered as the total of the shear-induced turbulent viscosity and the bubble-induced turbulent viscosity. The viscosity of the liquid phase can be expressed as

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

(10)

where

μ_{l}

is the dynamic viscosity of the continuous phase,

μ_{Tl}

is the turbulent viscosity of the continuous phase which is given by

μ_{Tl} = ρ_{l} C_{μ} k_{l}^{2} ɛ_{l}

(11)

and the extra bubble-induced turbulent viscosity

μ_{Tb}

is evaluated according to

μ_{Tb} = ρ_{l} C_{μ b} α_{g} D_{s} | u_{g} - u_{l} |

(12)

where D_s is the local bubble Sauter mean diameter. The constants

C_{μ}

and

C_{μ b}

have values of 0.09 and 0.6, respectively.

The effective viscosity of the disperse phase can be obtained as

μ_{g}^{eff} = (μ_{Tl} + μ_{Tb} σ_{g} ρ_{g} ρ_{l}) + μ_{g}^{lam}

(13)

where

σ_{g}

, holding the value of unity, is the turbulent Prandtl number of the disperse phase.

The turbulence variables k and ɛ are determined via the extended version of the two-equation single-phase standard $k - ɛ$ model.²³

Population balance model

Particle (bubble) size distribution can be determined in accordance with the population balance equation expressed in integro-differential form via

\partial f (x, ξ, t) \partial t + \nabla . (V (x, ξ, t) f (x, ξ, t)) = S (x, ξ, t)

(14)

where

f (x, ξ, t)

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

V (x, ξ, t)

is the velocity vector in external space dependent on the external variables x for a given time t and the internal space ξ whose components could be characteristic dimensions such as volume, mass etc. On the right hand side, 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, and phase change.

By taking mass M as the independent coordinate, the discrete particle number density can be defined by

N_{i} (x, t) = \int_{v_{i}}^{v_{i + 1}} f (x, M, t) d v i = 0, 1, 2, \dots, N

(15)

From above, the particle number density equation can be alternatively expressed in terms of size fraction f_i of N bubble size groups as

\partial (ρ_{g} α_{g} f_{i}) \partial t + \nabla . (ρ_{g} α_{g} u_{g} f_{i}) = S_{i}

(16)

In the above equation, S_i represents the net change in the number density distribution due to coalescence and break-up processes. This entails the use of a fixed non-uniform volume distribution along a grid that allows a range of bubble sizes to be covered with a small number of bins and offers good resolution. Such discretization of the population balance equation has been found to allow accurate determination of the desired characteristics of the number density distribution. The interaction term $S_{i} = (B_{C} + B_{B} + D_{C} + D_{B})$ contains the source rates of birth rates due to coalescence (B_C) and break-up (B_B) and death rates due to coalescence (D_C) and break-up (D_B) of bubbles, respectively.

Coalescence kernels

For coalescence between fluid particles, the coalescence rate $a (M_{i}, M_{j})$ could be calculated as a product of collision frequency $h (M_{i}, M_{j})$ and coalescence efficiency $λ (M_{i}, M_{j})$ :

a (M_{i}, M_{j}) = h (M_{i}, M_{j}) λ (M_{i}, M_{j})

(17)

The coalescence kernels that are adopted in this paper are described in the following.

Coulaloglou and Tavlarides ¹⁵ : This model has been developed based on the consideration of turbulent random motion-induced collisions as primary source of bubble coalescence. The collision frequency was defined as the effective swept volume rate of sizes of d_i and d_j

h (M_{i}, M_{j}) = π 4 (d_{i} + d_{j})^{2} (u_{ti}^{2} + u_{tj}^{2})^{12}

(18)

The turbulent velocity u_t in the inertial sub range of isotropic turbulence given by $u_{t} = C_{1} (ɛ_{l} d)^{13}$ where C₁ is a constant. As only a fraction of collisions will lead to coalescence, the coalescence efficiency based on the film drainage model for deformable particle with immobile surfaces can be written as

λ (M_{i}, M_{j}) = exp [- C_{2} \times μ_{l} ρ_{l} ɛ_{l} σ^{2} (d_{i} d_{j} d_{i} + d_{j})^{4}]

(19)

where C₂ is a constant.

Prince & Blanch ¹⁶ : The coalescence process can been described in three steps. Firstly, the bubbles trap small amount of liquid between them. Secondly, the liquid drains out until the liquid film thickness reaches a critical thickness. Thirdly, the film ruptures and bubbles coalesce together. Analogously to Coulaloglou and Tavlarides,¹⁵ the collision frequency is calculated as

h (M_{i}, M_{j}) = C_{3} (d_{i} + d_{j})^{2} (d_{i}^{23} + d_{j}^{23})^{12} ɛ_{l}^{13}

(20)

where C₃ is a constant. The coalescence efficiency for deformable particle with mobile surfaces can be written in accordance with

λ (M_{i}, M_{j}) = \exp (- t_{ij} τ_{ij})

(21)

where

t_{ij}

is the time due to intervening film reaching a critical thickness before rupturing and

τ_{ij}

is the time due to contact of bubble–bubble collision.

Lehr et al. ¹⁷ : In contrast to the above, the coalescence frequency has been considered based on the critical approach velocity model. A critical velocity is defined as the maximum velocity of bubbles resulting in coalescence which has no dependency on the size of the bubbles. Collisions will result in coalescence only when the relative approach velocity of bubbles perpendicular to the surface of contact is lower than the critical approach velocity. The collision frequency function is given by

h (M_{i}, M_{j}) = π 4 (d_{i} + d_{j})^{2} min (u', u_{critical}) \times \exp [- (α_{max}^{13} α_{g}^{13} - 1)^{2}] α_{max} = 0.6

(22)

The characteristic velocity $u'$ is equivalent to the turbulent eddy velocity with the similar length scale of the bubbles. Smaller eddies would not have sufficient energy to have significant impact on bubbles to collide. On the other hand, larger eddies would end up transporting the bubbles. For the larger eddies, characteristic velocity has been defined as the difference in rise velocities of the bubbles, viz.

u' = max (\sqrt{2} ɛ^{13} \sqrt{d_{i}^{23} + d_{j}^{23}}, | u_{i} - u_{j} |)

(23)

where

u_{i}

and

u_{j}

are the velocity vectors of the discrete bubbles travelling in the liquid flow. The coalescence efficiency is given in accordance with

λ (M_{i}, M_{j}) = min (u_{critical} u', 1)

(24)

The critical velocity $u_{critical}$ can be determined according to

u_{critical} = \sqrt{{We}_{critical} σ ρ_{l} d_{eq}}

(25)

where

{We}_{critical}

is the critical Weber number and

d_{eq}

is the equivalent diameter given by

d_{eq} = 2 d_{i} d_{j} / (d_{i} + d_{j})

Break-up kernels

For the break-up of fluid particles, the break-up rate $r (M_{i}, M_{j})$ is a function of total break-up frequency $r (M_{i})$ and DSD $β (M_{i}, M_{j})$

β (M_{i}, M_{j}) = r (M_{i}, M_{j}) r (M_{i})

(26)

Luo and Svendsen ¹³ : The assumption of bubble binary break-up is invoked under isotropic turbulence situation. Break-up event is determined by the energy level of arriving eddy with smaller or equal length scale compared to the bubble diameter to induce the oscillation. The DSD is accounted using a stochastic break-up volume fraction

f_{BV}

. The break-up rate in terms of mass can be obtained as

r (M_{i}, M_{j}) = 0.923 (1 - α_{g}) n (ɛ_{l} d_{j})^{13} \times \int_{ξ_{min}}^{1} (1 + ξ)^{2} ξ^{113} P_{b} (f_{BV} | d_{j}, λ) d ξ

(27)

The breakage probability, $P_{b} (f_{BV} | d_{j}, λ)$ is calculated by using the energy distribution of turbulent eddies with size λ. Hence, the breakage rate becomes

r (M_{i}, M_{j}) = 0.923 (1 - α_{g}) n (ɛ_{l} d_{j})^{13} \times \int_{ξ_{min}}^{1} (1 + ξ)^{2} ξ^{113} exp (- 12 c_{f} σ β ρ_{l} ɛ_{l}^{23} d_{j}^{53} ξ^{113}) d ξ

(28)

The break-up frequency represents the break-up bubble with mass of M_i into fraction of $f_{BV}$ and $f_{BV} + {df}_{BV}$ for a continuous $f_{BV}$ function. The total break-up frequency can thus be obtained by integrating over the whole interval of 0 to 1 as

r (M_{i}) = 12 \int_{0}^{1} r (M_{i}, M_{j}) d f_{BV}

(29)

with a daughter bubble size distribution expressed as

β (f_{BV}, 1) = \int_{ξ_{min}}^{1} (1 + ξ)^{2} ξ^{113} \exp (- 12 c_{f} σ β ρ_{l} ɛ_{l}^{23} d_{j}^{53} ξ^{113}) d ξ 12 \int_{0}^{1} \int_{ξ_{min}}^{1} (1 + ξ)^{2} ξ^{113} \exp (- 12 c_{f} σ β ρ_{l} ɛ_{l}^{23} d_{j}^{53} ξ^{113}) d ξ d f_{BV}

(30)

Wang et al. ¹⁴ : Considering the energy constraint of the break-up kernel developed by Luo and Svendsen,¹³ the model is further extended by adding the capillary constraint to calculate the break-up of bubbles. In this model, the dynamic pressure of the turbulent eddy is taken to be larger than the capillary pressure resulting in minimum break-up fraction. On the other hand, the eddy kinetic energy also needs to be larger than the increase of the surface energy resulting in maximum break-up. The advantage of this model is to have no adjustable parameter and provide the DSD directly by normalizing the partial break-up frequency by the total frequency:

r (M_{i}, M_{j}) = 0.923 (1 - α_{g}) n ɛ_{l}^{13} \times \int_{λ_{min}}^{d_{i}} P_{b} (f_{BV} | d_{j}, λ) (λ + d)^{2} λ^{113} d λ

(31)

The total break-up frequency is also calculated according to equation (29) but with a different daughter bubble size distribution expressed as where $e (λ)$ is the kinetic energy in the inertial sub range and $\bar{e} (λ)$ is the mean kinetic energy.

The DSD profiles for Luo and Svendsen¹³ and Wang et al.¹⁴ for various bubble sizes are illustrated in Figure 1.

Figure 1.

Daughter size distribution profiles: (a) Luo and Svendsen¹³ and (b) Wang et al.¹⁴

Fractal wall heat flux partitioning model

A fractal model for the wall heat flux partitioning model is derived based on the concept that the nucleation site size distribution follows the fractal power law. The number of active cavities for the sizes of cavities between D_c and D_c+dD_c can be obtained as $- {dN}_{a} = d_{f} (D_{c, max}^{d_{f}} / D_{c}^{d_{f} + 1}) {dD}_{c}$ where ${dD}_{c} > 0$ and $- {dN}_{a} > 0$ . The different heat flux components for the modified heat flux partitioning model are derived in the following.

Surface quenching

The process of surface quenching or transient conduction occurs in regions that are swept by sliding bubbles, Q_tcsl, or in regions at the point of inception, Q_tc. A fractal model for the heat flux from the smallest site D_c_,min to the largest site D_c_,max is given by

no Q_{tc} = \int_{D_{c, min}}^{D_{c, max}} 2 \sqrt{k_{l} ρ_{l} C_{pl} π (t_{w} + t_{sl})} (T_{\sup} + T_{sub}) \times R_{f} (K π D_{sl}^{2} 4) t_{w} f (- {dN}_{a}) + \int_{D_{c, min}}^{D_{c, max}} 2 \sqrt{k_{l} ρ_{l} C_{pl} π (t_{w} + t_{sl})} (T_{\sup} + T_{sub}) \times R_{f} (π D_{sl}^{2} 4) (- t_{w} f) (- {dN}_{a})

(33)

where t_w is the waiting time,

t_{sl}

is the sliding time of departed bubbles,

T_{\sup}

is the superheat temperature, R_f is the ratio of the actual number of bubbles lifting off per unit area of the heater surface to the number of active nucleation sites per unit area, K is the ratio of the area of influence (the area from where the liquid is drawn in when the bubble leaves the heater surface) to the projected area of the bubble,

D_{sl}

is the size of the sliding bubbles, f is the bubble frequency, and N_a is the active nucleation site density.

The transient conduction that takes place during the sliding phase and the area occupied by the sliding bubble at any instant of time from the smallest site D_c_,min to the largest site D_c_,max is nonetheless given by the fractal model for the heat flux as

Q_{tcsl} = \int_{D_{c, min}}^{D_{c, max}} 2 \sqrt{k_{l} ρ_{l} C_{pl} π (t_{w} + t_{sl})} (T_{\sup} + T_{sub}) \times R_{f} l_{s} {KDt}_{w} f (- {dN}_{a}) + \int_{D_{c, min}}^{D_{c, max}} 2 \sqrt{k_{l} ρ_{l} C_{pl} π (t_{w} + t_{sl})} (T_{\sup} + T_{sub}) \times R_{f} {ft}_{sl} (π D^{2} 4) (1 - t_{w} f) (- {dN}_{a})

(34)

where D is the average bubble diameter given by

D = (D_{sl} + D_{l}) / 2

Evaporation

The heat flux due to vapour generation occurs at the nucleate boiling region which is calculated by the energy carried away by the bubbles lifting off from the heated surface. A fractal model for this heat flux from the smallest site D_c_,min to the largest site D_c_,max is given by

Q_{e} = \int_{D_{c, min}}^{D_{c, max}} R_{f} f (π D_{l}^{3} 6) ρ_{g} h_{fg} (- {dN}_{a})

(35)

where D_l is the size of the bubbles leaving the heated wall.

Turbulent convection

Based on fractal characteristics, there is an expression that could relate the pore volume fraction to fractal dimension, minimum, and maximum pore size in porous media. This expression is given by $ψ = (D_{c, min} / D_{c, max})^{d - d_{f}}$ where d = 2 in a two-dimensional space. The forced convection will always prevail at all times in areas of the heater surface that are not influenced by the stationary and sliding bubbles. This heat flux can be obtained as

Q_{c} = St ρ_{l} C_{pl} u_{l} (1 - K ψ) (T_{\sup} + T_{sub})

(36)

where

St

is the Stanton number and u_l is the adjacent liquid velocity.

The total wall heat flux Q_w is thus obtained as sum of the aforementioned heat flux components: $Q_{w} = Q_{tc} + Q_{tcsl} + Q_{e} + Q_{c}$ . More details of the model can be referred in Yeoh et al.⁹

Experimental details

In order to assess the vapor distribution in the radial direction for low and medium pressures, three experiments are investigated. Experimental conditions for low pressure (Cases P143) and elevated pressure (Cases P218, P497 and P949) data are presented in Table 2. These cases cover a range of different flow conditions including pressure, inlet liquid velocity, wall heat flux, and inlet subcooling temperature that play important roles on vapor phase distribution and wall heat flux partitioning. The authors tried to illustrate the underlying physics through the results obtained by simulations. For each case, simulation results are validated against available data of these experiments. To help the readers understand the experimental conditions investigated in this paper, the details of experiments are given as follows. For more details refer to the references cited below.

Table 1.

Consideration of different kernel combinations.

No.	Coalescence kernel	Break-up kernel
1	Prince and Blanch¹⁶	Luo and Svendsen¹³
2	Prince and Blanch¹⁶	Wang et al.¹⁴
3	Coulaloglou and Tavlarides¹⁵	Luo and Svendsen¹³
4	Coulaloglou and Tavlarides¹⁵	Wang et al.¹⁴
5	Lehr et al.¹⁷	Luo and Svendsen¹³
6	Lehr et al.¹⁷	Wang et al.¹⁴

Table 2.

Experimental conditions for specific cases.

Case	P_inlet (kPa)	T_inlet (℃)	T_sub@inlet (℃)	Q_w (kW/m²)	G (kg/m²s)
P143	143	92.1	17.9	251.5	1059.2
P218	218	110.3	12.7	237.9	1843.8
P497	497	136.7	14.8	190.9	942.3
P949	949	167.6	10.0	208.5	964.4

Low pressure experiment performed by Yun et al.²⁴ and Lee et al.²⁵ consisted of a vertical concentric annulus with an inner diameter of 37.5 mm for the outer wall, and outer diameter of 19 mm for the inner heating rod as the test section; the working fluid was demineralized water. The heated section was 1.67 m long and entire rod was heated by a 54 kW DC power supply. Radial measurements of phasic parameters were done at 1.61 m downstream of the start of the heated section. A two-conductivity probe method was used to measure local gas phase parameters such as local void fraction, bubble frequency, and bubble velocity. The bubble Sauter mean diameters (assuming spherical bubbles) were determined through the interfacial area concentration (IAC), calculated using the measured bubble velocity spectrum and bubble frequency. The uncertainties in the measurement of local void fraction, velocity, volumetric flow rate, temperature, heat flux, and pressure are estimated to be within ±3.0%, ±3.3%, ±1.9%, ±0.2℃, ±1.7%, and ±0.0005 MPa, respectively.

Ozar et al.¹⁸ performed medium pressure experiments where a vertical concentric annulus was employed. The outer wall’s inner diameter was 38.1 mm, and the inner heating rod had 19.1 mm outer diameter. The annulus was designed between the pipes and the cartridge heater. The heated section was 2.845 m long which was followed by a 1.632 m long unheated section. The heater could produce a maximum heat flux of 260 kW/m². The measurements presented in this paper, were performed at 1.06, 2.05, and 2.83 m downstream of the start of the heated section. The uncertainties in the measurement of local void fraction (done through a 4-sensor conductivity probe), gas velocity, flow rate, temperature, and pressure are estimated to be less than 10%, less than 10%, within ±0.75%, ±2.2℃, and less than ±0.2%, respectively.

Results and discussion

In order to discretize the conservation equations of mass, momentum and energy, the finite volume method is employed. Mentioned equations for each phase along with 15 extra set of transport equations for capturing coalescence, break-up and condensation of the bubbles for the MUSIG boiling model are solved. Since a uniform wall heat flux is applied, only a 60° section of the annulus is modeled as the computational domain for all the cases. Grid independence is inspected for 45, 90, 180, 240, and 300 cells along the vertical direction, and 5, 10, 20, and 30 cell in the radial direction; the mean velocity profiles of liquid and gas, and the volume fraction distribution did not change significantly by further grid refinement of 180 cells in the vertical direction and 10 cells in the radial direction. The proposed mechanistic approach along with some of the existing empirical correlations is compared against experimental data of Yun et al.²⁴ and Lee et al.²⁵ for Case P143 and Ozar et al.¹⁸ for Cases P218–P949. The proposed mechanistic model consists of fractal wall heat flux partitioning model. For the break-up kernels, two widely adopted models with different predictions for DSD proposed by Luo and Svendsen¹³ and Wang et al.¹⁴ are selected. These break-up kernels are then coupled with three different coalescence kernels by Coulaloglou and Tavlarides,¹⁵ Prince and Blanch,¹⁶ and a more recent one by Lehr et al.¹⁷ to form six different combinations of kernels.

Bubble sauter mean diameter profiles

In Figures 2 to 4, the predicted bubble Sauter mean diameter profiles in the radial direction for six aforementioned kernels are presented against the experimental data of Yun et al.²⁴ and Lee et al.²⁵ for Case P143 and experiments of Ozar et al.¹⁸ for Cases P218–P949.

Figure 2.

Bubble Sauter mean diameter profiles for Case P143 at Z/Dh = 84.7 and Cases P218–P949 at Z/Dh = 108.

Figure 3.

Bubble Sauter mean diameter profiles for Cases P218–P949 at Z/Dh = 56.

Figure 4.

Bubble Sauter mean diameter profiles for Cases P218–P949 at Z/Dh = 149.

In general, one could notice that the coalescence kernels pose an insignificant contribution in the prediction of the bubble size. Among the three coalescence kernels, Coulaloglou and Tavlarides¹⁵ model tends to predict a higher rate of bubbles coalescence Lehr et al.¹⁷ model predicts a lower rate.

All the kernels predicting the bubble size closely agreed with the experimental measurement at the near heated wall region; however, moving away from the heated wall into the bulk liquid region, the kernels 2, 4, 6 with the same break-up kernel of Wang et al.¹⁴ produce considerably different predictions to the kernels 1, 3, 5 with break-up kernel of Luo and Svendsen.¹³

For the lower pressure cases (Cases P143–P497), the break-up kernel of Wang et al.¹⁴ tends to over-predict the bubble size in the subcooled region. In other words, the rate of bubble break-up predicted by the model is lower than that of Luo and Svendsen.¹³ Nevertheless, for the Case P949 where two group of bubbles (i.e. spherical and cap bubbles) are present, the Wang et al.¹⁴ kernel predicts a considerably better agreement with the experimental data. On the other note, one should also notice that the break-up kernels is not the only influential parameter dictating the bubble size evolution. The condensation in the subcooled region as well as the influence of different bubble shapes (rather than spherical) could also be significant but subject to further investigations.

Void fraction profiles

Figures 5 to 7 present the comparison of the predicted void fraction profiles in the radial direction for the six aforementioned kernels against the experimental data of Yun et al.²⁴ and Lee et al.²⁵ (i.e. Case P143) and the experimental data of Ozar et al.¹⁸ (i.e. Cases P218–P949).

Figure 5.

Void fraction profiles for Case P143 at Z/Dh = 84.7 and Cases P218–P949 at Z/Dh = 108.

For all cases, the trend of void fraction distribution is captured and well agreed with the experimental data. A higher void fraction near the heated wall is due to the vapor generation at the surface of the heated wall. In the bulk liquid region, where bubbles are exposed to the subcooled liquid, condensation of bubble becomes dominant, vanishing the local void fraction. However, an over-prediction of void fraction near the heated wall is observed. All six kernels predict closely with the measurements for the lower pressure cases (Cases P143–P497); yet, the kernels 2, 4, and 6 pose better predictions for the elevated pressure case (Case P949). In this Case, two groups of bubbles are present which leads to higher void fractions compared to other Cases. The lower break-up rate that is predicted by Wang et al.¹⁴ poses a better agreement with the experimental results in such case.

IAC profiles

IAC profiles in the radial direction for six kernels are depicted and compared against the experimental data of Yun et al.²⁴ and Lee et al.²⁵ (i.e. Case P143) and the experiment data of Ozar et al.¹⁸ (i.e. Cases P218–P949) in Figures 8 to 10. As depicted, the influence of different coalescence kernels is insignificant in the prediction of IAC profile for all cases.

Figure 6.

Void faction profiles for Cases P218–P949 at Z/Dh = 56.

Figure 7.

Void faction profiles for Cases P218–P949 at Z/Dh = 149.

Figure 8.

Interfacial area concentration profiles for Case P143 at Z/Dh = 84.7 and Cases P218–P949 at Z/Dh = 108.

Figure 9.

Void faction profiles for Cases P218–P949 at Z/Dh = 56.

Figure 10.

Void faction profiles for Cases P218–P949 at Z/Dh = 56.

On the other hand, the Kernels 1, 3, 5 with Luo and Svendsen’s¹³ break-up model tend to over-predict the IAC at the near heated wall region. Meanwhile, the Kernels 2, 4, 6 with Wang et al.’s¹⁴ break-up model pose a considerably better predictions of IAC at the vicinity of the heated wall. The over-prediction of IAC in Luo and Svendsen’s¹³ model in conjunction with the over-prediction of void fraction (as was observed in Figure 2, especially for the Case P497), leads to a better prediction of the bubble size (as was observed in Figure 1) compared to the Kernels with Wang et al.’s¹⁴ break-up model.

Similar to other radial profiles, the Wang et al.’s¹⁴ model poses better prediction for the IAC profile at the elevated pressure case (Case P949). This could be attributed to the formulation of the Wang et al.’s model where bubble break-up only occurs if the dynamic pressure of the approaching turbulent eddy is higher than the capillary pressure of bubbles. Therefore, the influence of pressure could be better captured in this model which leads to better prediction of all radial profiles of bubble Sauter mean diameter, void fraction, and IAC for the Case P949.

Conclusion

The performance of different coalescence and break-up kernels is investigated through numerical simulations of subcooled boiling flows at elevated pressures. Numerical predictions are validated against the experimental data of experiments of Yun et al.²⁴ and Lee et al.²⁵ for Case P143 and experiments of Ozar et al.¹⁸ for Cases P218, P497, and P949. Overall, the bubble size and void fraction profiles’ trends are reasonably captured through these kernels. The influence of different coalescence kernels investigated in this study has been found to be insignificant in comparison to the break-up kernels. The model by Luo and Svendsen¹³ seems to predict a higher rate of break-up, resulting in a better prediction of bubble size and void fraction for the lower pressure cases. Nonetheless, the consideration of capillary pressure in the break-up model by Wang et al.¹⁴ gave better predictions for the elevated pressure case.

Footnotes

Acknowledgements

This paper is an extension of a conference paper that was presented at the Eleventh International Conference on CFD in the Minerals and Process Industries (CFD2015) and was nominated for invitation into the CFD2015 Special Issue of Journal of Computational Multiphase Flows based on its designation as having very good scope of extension into a high-quality paper of relevance to multiphase flows.

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The financial support provided by the Australian Research Council, Australia (ARC project ID DP130100819) is gratefully acknowledged.

Appendix

References

Basu

Warrier

Dhir

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

Basu

Warrier

Dhir

. Wall heat flux partitioning during subcooled flow boiling: part II – model validation. ASME J Heat Transf 2005b; 127: 141–148.

Cole

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

Lemmert

Chawla

. Influence of flow velocity on surface boiling heat transfer coefficient. Heat Transf Boil 1977; 237: 247.

Tolubinsky V and Kostanchuk D. Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling. In: Proceedings of the 4th international heat transfer conference, vol. 5, Paper no. B2.8, Paris-Versailles, France, 1970.

Ü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/m2. Int J Heat Mass Trans 1976; 19: 643–649.

Cheung

SCP

Vahaji

Yeoh

. Modeling subcooled flow boiling in vertical channels at low pressures – part 1: assessment of empirical correlations. Int J Heat Mass Transf 2014; 75: 736–753.

Vahaji S, Cheung SCP, Tu JY, et al. Modeling subcooled boiling flow in low and medium pressure – evaluation of mechanistic approach. OECD/NEA & IAEA Workshop, CFD4NRS-5, 2014.

Yeoh

Vahaji

Chueng

SCP

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

10.

Liao

Lucas

. A literature review of theoretical models for drop and bubble breakup in turbulent dispersions. Chem Eng Sci 2009; 64: 3389–3406.

11.

Liao

Lucas

. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem Eng Sci 2010; 65: 2851–2864.

12.

Deju

Cheung

SCP

Yeoh

. Comparative analysis of coalescence and breakage kernels in vertical gas-liquid flow. Canad J Chem Eng 2015; 93: 1295–1310.

13.

Luo

Svendsen

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

14.

Wang

Jin

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

15.

Coulaloglou

Tavlarides

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

16.

Prince

Blanch

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

17.

Lehr

Millies

Mewes

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

18.

Ozar

Brooks

Hibiki

. Interfacial area transport of vertical upward steam-water two-phase flow in an annular channel at elevated pressures. Int J HeatMass Trans 2013; 57: 504–518.

19.

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.

20.

Ranz

Marshall

. Evaporation from droplets: Parts I and II. Chem Eng Prog 1952; 48: 141–148.

21.

Anglart

Nylund

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

22.

Sato

Sadatomi

Sekoguchi

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

23.

Yeoh

Chueng

SCP

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

24.

Yun BJ, Park G, G.-C., C. H. Song and M. K. Chung, Measurements of local two-phase flow parameters in a boiling flow channel. In: Proceedings of the OECD/CSNI specialist meeting on advanced instrumentation and measurement techniques, Santa Barbara, CA, USA. 17–20 March 1997.

25.

Lee

Park

Lee

. Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus. Int J Multiphase Flow 2002; 28: 1351–1368.

Numerical investigation on the performance of coalescence and break-up kernels in subcooled boiling flows in vertical channels

Abstract

Keywords

Introduction

Mathematical formulation

Two-fluid model

Population balance model

Coalescence kernels

Break-up kernels

Fractal wall heat flux partitioning model

Surface quenching

Evaporation

Turbulent convection

Experimental details

Results and discussion

Bubble sauter mean diameter profiles

Void fraction profiles

IAC profiles

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

Appendix

References