Science and Technology of Nuclear Installations
Volume 2009 (2009), Article ID 387020, 10 pages
doi:10.1155/2009/387020
Research Article

Mechanistic Multidimensional Modeling of Forced Convection Boiling Heat Transfer

Michael Z. Podowski1 and Raf M. Podowski2

1Center for Multiphase Research, Rensselaer Polytechnic Institute, Troy, NY 12180, USA
2SeeqPod, Inc., Emeryville, CA 94608, USA

Received 1 June 2008; Accepted 22 August 2008

Academic Editor: Iztok Tiselj

Copyright © 2009 Michael Z. Podowski and Raf M. Podowski. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Due to the importance of boiling heat transfer in general, and boiling crisis in particular, for the analysis of operation and safety of both nuclear reactors and conventional thermal power systems, extensive efforts have been made in the past to develop a variety of methods and tools to evaluate the boiling heat transfer coefficient and to assess the onset of temperature excursion and critical heat flux (CHF) at various operating conditions of boiling channels. The objective of this paper is to present mathematical modeling concepts behind the development of mechanistic multidimensional models of low-quality forced convection boiling, including the mechanisms leading to temperature excursion and the onset of CHF.

1. Introduction

Because of the complexity of phenomena governing boiling heat transfer in general, and subcooled boiling in particular, the predictions of CHF have traditionally been based on correlating data obtained from numerous experimental measurements. At the same time, our understanding of the local physical mechanisms governing the near-wall phase change and transport has been steadily improving. Given the progress made in the computational fluid dynamics methods, the possibility of using complete multidimensional models to predict boiling heat transfer before, and up to, the onset of CHF becomes an attractive option complementing the traditional phenomenological approach used in the past.

The objective of this paper is to discuss various physical and mathematical modeling concepts for local heat transfer phenomena in boiling systems, and to show that the proposed approach can be combined with mechanistic multidimensional models of two-phase flow and used to predict various parameters characterizing low-quality forced-convection boiling, including the mechanisms leading to temperature excursion and the onset of CHF.

2. Multidimensional Multifield Model of Forced-Convection Boiling

The multifield modeling concept is based on coupling a complete mechanistic multidimensional model of two-phase flow with the models governing local heat transfer and phase change phenomena in heated channels. The major components of the overall model are the following: (a)conservation equations for the individual fields,(b)the model of turbulence,(c)closure laws for the interfacial mass transfer,(d)closure laws for the interfacial momentum transfer,(e)closure laws for the interfacial energy transfer,(f)kinematic boundary conditions for phasic velocities and turbulence,(g)thermal boundary conditions for the near-wall heat transfer.

2.1. Multifield Conservation Equations

The multifield model of two- and multiphase flows assumes that each phase can be represented by either one or several fields. In the latter case, the individual fields represent topologically different flow structures within a given phase (such as liquid film and liquid droplets or groups of gas bubbles of different sizes). Typical ensemble-averaged conservation equations for a multifield model can be written as follows:

Mass 𝜕 𝛼 𝑘 𝜌 𝑘 𝛼 𝜕 𝑡 + ∇ ⋅ 𝑘 𝜌 𝑘 𝐮 𝑘 = Γ 𝑘 , ( 1 )

Momentum 𝜕 𝛼 𝑘 𝜌 𝑘 𝐮 𝑘 𝛼 𝜕 𝑡 + ∇ ⋅ 𝑘 𝜌 𝑘 𝐮 𝑘 𝐮 𝑘 𝛼 = − ∇ 𝑘 𝑝 𝑘 𝛼 + ∇ ⋅ 𝑘 𝝉 t o t 𝑘 + 𝛼 𝑘 𝜌 𝑘 𝐌 𝐠 + 𝑖 𝑘 , ( 2 )

Energy 𝜕 𝛼 𝑘 𝜌 𝑘 𝑒 𝑘 𝛼 𝜕 𝑡 + ∇ ⋅ 𝑘 𝜌 𝑘 𝐮 𝑘 𝑒 𝑘 𝛼 = − ∇ ⋅ 𝑘 𝐪 t o t 𝑘 𝛼 + ∇ ⋅ 𝑘 − 𝑝 𝑘 𝐈 + 𝝉 t o t 𝑘 ⋅ 𝐮 𝑘 + 𝛼 𝑘 𝜌 𝑘 𝐠 ⋅ 𝐮 𝑘 + 𝐸 𝑖 𝑘 , ( 3 ) where 𝛼 𝑘 , 𝜌 𝑘 , 𝐮 𝑘 , and 𝑒 𝑘 are the volume fraction, density, velocity, and specific internal energy, respectively, of field- 𝑘 , Γ 𝑘 is the net mass transfer rate (this and all other terms are per unit volume) of field- 𝑘 , 𝝉 t o t 𝑘 is the total shear stress, 𝐌 𝑖 𝑘 is the total interfacial force, 𝐪 t o t 𝑘 is the local heat flux, and 𝐸 𝑖 𝑘 is the interfacial heat transfer rate, all for field- 𝑘 .Details concerning the various terms in (1)–(3) can be found in [1].

In the modeling of low-quality nucleate boiling flows, a two-field model can normally be used, which accounts for the continuous liquid ( 𝑙 ) and dispersed vapor ( 𝑣 ) fields. This is shown in Figure 1(a), where upon departing from the heated wall, bubbles move into the subcooled liquid region and condense.

fig1
Figure 1: An illustration of the multifield model of two-phase flow in a boiling channel.

For high wall heat fluxes, vapor concentration near the heated wall increases, so that the individual departing bubbles may coalesce and form elongated bubbles, as shown in Figure 1(b). Since the mechanisms governing the motion and heat transfer of elongated bubbles may be quite different from those for small spherical bubbles (field- 𝑣 1 ), it is appropriate to use a separate field for the former, denoted here as field- 𝑣 2 .

The multidimensional multifield conservation equations, (1)–(3), must be complemented by appropriate boundary conditions and closure laws. The boundary conditions of particular interest to this work are associated with the effect of heated channel wall on local flow and phase change in the near-wall region. The necessary closure laws and models are those for turbulence, interfacial forces, interfacial mass transfer, and heat transfer. A brief description of the major closure laws and boundary conditions pertinent to low-quality two-phase flows in boiling channels is given below.

2.2. Multifield Model of Turbulence Model

Two-phase flow turbulence is normally modeled using the 𝜅 - 𝜀 model [2] for the continuous liquid field ( 𝑘 = 𝑙 ), modified to include the effect of bubble-induced turbulence. The turbulent kinetic energy, 𝜅 𝑙 , and the energy dissipation, 𝜀 𝑙 , of the liquid field are, respectively, given by 𝜕 𝛼 𝑙 𝜅 𝑙 𝛼 𝜕 𝑡 + ∇ ⋅ 𝑙 𝐮 𝑙 𝜅 𝑙 𝛼 = ∇ ⋅ 𝑙 𝜇 𝑡 𝑙 𝜌 𝑙 𝜎 𝜅 ∇ 𝜅 𝑙 + 𝑆 𝜅 𝑙 , 𝜕 𝛼 𝑙 𝜀 𝑙 𝛼 𝜕 𝑡 + ∇ ⋅ 𝑙 𝐮 𝑙 𝜀 𝑙 𝛼 = ∇ ⋅ 𝑙 𝜇 𝑡 𝑙 𝜌 𝑙 𝜎 𝜀 ∇ 𝜀 𝑙 + 𝑆 𝜀 𝑙 . ( 4 )

For given 𝜅 𝑙 𝜀 𝑙 , the turbulent viscosity, 𝜇 𝑡 𝑙 , can be expressed as 𝜇 𝑡 𝑙 = 𝐶 𝜇 𝜌 𝑙 𝜅 𝑒 𝑙 𝜀 𝑙 + 𝜇 𝑡 𝐵 , ( 5 ) where the bubble-induced turbulence can be written as [3] 𝜇 𝑡 𝐵 = 𝐶 𝜇 𝐵 𝛼 𝑙 𝜌 𝑙 2 𝑘 = 1 𝐮 𝑣 𝑘 − 𝐮 𝑙 . ( 6 )

2.3. Interfacial Mass Transfer

Since in nucleate boiling, evaporation occurs only at (or near) the heated wall; the interfacial mass transfer between the liquid and vapor fields is practically only due to vapor condensation in contact with subcooled liquid. In addition, after elongated bubbles start forming, another possible interfield mass exchange mechanism is due to small bubble coalescence. Hence, the volumetric mass transfer terms can be expressed as Γ 𝑙 = − Γ e v a p + Γ c o n d Γ = − 𝑣 1 + Γ 𝑣 2 , Γ ( 7 ) 𝑣 1 = Γ e v a p 𝑣 1 − Γ c o n d 𝑣 1 − Γ c l s c 𝑣 1 , Γ ( 8 ) 𝑣 2 = Γ e v a p 𝑣 2 − Γ c o n d 𝑣 2 + Γ c l s c 𝑣 1 . ( 9 )

In (7), Γ e v a p = Γ e v a p 𝑣 1 + Γ e v a p 𝑣 2 and Γ c o n d = Γ c o n d 𝑣 1 + Γ c o n d 𝑣 2 are the total evaporation and condensation rates per unit volume, respectively, and Γ c l s c 𝑣 1 is the coalescence rate of small bubbles per unit volume (which becomes a source term for the elongated bubble field).

2.4. Interfacial Momentum Transfer

In general, the total interfacial force on phase- 𝑘 can be expressed as a superposition of several component forces: 𝐌 𝑖 𝑘 = 𝑝 𝑖 𝑘 ∇ 𝛼 𝑘 − 𝝉 𝑖 𝑘 ⋅ ∇ 𝛼 𝑘 + Γ 𝑘 𝐮 𝑖 𝑘 + 𝑗 𝐌 𝑖 𝑘 , 𝑗 . ( 1 0 )

In dispersed bubbly flows, the major interfacial forces are the following (for details and source references, see [1]):

Drag Force 𝐌 𝐷 𝑙 , 𝑣 𝑘 = − 𝐌 𝐷 𝑣 𝑘 , 𝑙 1 = − 8 𝐶 𝐷 , 𝑘 𝐴 𝑣 𝑘 𝜌 𝑙 𝐮 𝑣 𝑘 − 𝐮 𝑙 | | 𝐮 𝑣 𝑘 − 𝐮 𝑙 | | , ( 1 1 )

Virtual Mass Force 𝐌 V M 𝑙 , 𝑣 𝑘 = − 𝐌 V M 𝑣 𝑘 , 𝑙 = − 𝐶 V M , 𝑘 𝛼 𝑣 𝑘 𝜌 𝑙 𝐷 𝑣 𝑘 𝐮 𝑣 𝑘 − 𝐷 𝐷 𝑡 𝑙 𝐮 𝑙 𝐷 𝑡 , ( 1 2 )

Lift Force 𝐌 𝐿 𝑙 , 𝑣 𝑘 = − 𝐌 𝐿 𝑣 𝑘 , 𝑙 = − 𝐶 𝐿 , 𝑘 𝛼 𝑣 𝑘 𝜌 𝑙 𝐮 𝑣 𝑘 − 𝐮 𝑙 × ∇ × 𝐮 𝑙 , ( 1 3 )

Turbulent Dispersion Force 𝐌 T D 𝑙 , 𝑣 𝑘 = − 𝐌 T D 𝑣 𝑘 , 𝑙 = 𝐶 T D , 𝑘 𝛼 𝑣 𝑘 𝜌 𝑙 ∇ 𝛼 𝑙 𝜅 𝑙 ( 1 4 ) for 𝑘 = 1 , 2 , where 𝛼 𝑙 = 1 − ( 𝛼 𝑣 1 + 𝛼 𝑣 2 ) .

As it can be seen from (11)–(14), direct momentum exchange occurs between the continuous liquid and each vapor field, whereas bubble-bubble interactions are due to the coalescence of small bubbles to form elongated bubbles (see Section 2.5 below) and they mainly involve mass and energy transfers.

2.5. Interfacial Energy Transfer

Given a very low-vapor superheat, the interfacial volumetric condensation rates can be expressed as Γ c o n d 𝑘 = 𝑞 𝑘 𝑖 ℎ 𝑓 𝑔 𝑘 = 𝑣 1 , 𝑣 2 , ( 1 5 ) where the interfacial energy transport from vapor field-k to the continuous liquid field is given by 𝑞 𝑘 𝑖 = 𝐻 𝑖 𝑘 − 𝑙 𝑇 s a t − 𝑇 𝑙 𝐴 𝑘 𝑘 = 𝑣 1 , 𝑣 2 . ( 1 6 )

In (16), 𝐴 𝑘 is the interfacial area density for vapor field- 𝑘 , and 𝐻 𝑖 𝑘 − 𝑙 is the local interfacial heat transfer coefficient that can be obtained from 𝐻 𝑖 𝑘 − 𝑙 = 𝑘 𝑙 𝑑 𝑘 𝜌 2 + 0 . 6 𝑙 𝑑 𝑘 𝑣 𝑟 , 𝑘 𝜇 𝑙 0 . 5 𝑐 𝑝 𝑙 𝜇 𝑙 𝑘 𝑙 0 . 3 3 , 𝑘 = 𝑣 1 , 𝑣 2 . ( 1 7 )

Equations (16) and (17) imply 𝑞 𝑙 𝑖 = 𝑞 𝑣 𝑖 1 + 𝑞 𝑣 𝑖 2 = 𝐻 𝑖 𝑣 1 − 𝑙 𝐴 𝑣 1 + 𝐻 𝑖 𝑣 2 − 𝑙 𝐴 𝑣 2 𝑇 s a t − 𝑇 𝑙 . ( 1 8 )

The rate of bubble coalescence can be obtained from Γ c l s c 𝑣 1 = 𝑏 𝑣 1 − 𝑣 1 𝛼 2 𝑣 1 + 𝑏 𝑣 1 − 𝑣 2 𝛼 𝑣 1 𝛼 𝑣 2 , ( 1 9 ) where 𝑏 𝑣 1 − 𝑣 1 and 𝑏 𝑣 1 − 𝑣 2 are the corresponding rate coefficients.

3. Near-Wall Heat Transfer in Nucleate Boiling

In the nucleate subcooled boiling in a heated channel, the wall heat is partially used to form bubbles and the remaining portion is transferred to the liquid. The heat transfer from the wall in the vicinity of a nucleation site occurs during two distinct periods: the bubble growth time and the waiting time. The total convective heat flux from the wall is the sum of three models [4]: 𝑞 N B = 𝑞 1 𝜙 + 𝑞 𝑒 + 𝑞 𝑄 , ( 2 0 ) where 𝑞 1 𝜙 is the single-phase convective heat flux, 𝑞 𝑒 is the heat flux associated with phase change (evaporation), and 𝑞 𝑄 is the so-called quenching heat flux, which is transferred to the liquid phase during the waiting time.

Outside of the influence area of the bubbles, the heat transfer from wall to the liquid can be determined from 𝑞 1 𝜙 = 𝐴 1 𝜙 S t 𝑃 𝜌 𝑙 𝑐 𝑝 𝑙 𝑢 𝑙 , 𝑃 𝑇 𝑤 − 𝑇 𝑙 , 𝑃 , ( 2 1 ) where 𝐴 1 𝜙 is the fraction of the wall unaffected by the nucleation sites, S t 𝑃 is the Stanton number calculated from a heat transfer correlation in terms of the local liquid velocity and Prandtl number, 𝑇 𝑤 , is the wall temperature, and 𝑇 𝑙 , 𝑃 is the local liquid temperature near the heated wall. In numerical calculations, the local fluid properties are normally determined at the center of the near-wall computational cell. As long as a consistent turbulence model is used (see Section 2.2), the calculated single-phase heat flux component is independent of the distance between point- 𝑃 and the wall, or, in other words, is grid-independent.

The evaporation heat flux is given by 𝑞 𝑒 = 𝜋 6 𝑑 3 d e t 𝜌 𝑣 𝑓 d e t 𝑁 ℎ 𝑓 𝑔 , ( 2 2 ) where 𝑑 d e t is the critical bubble diameter at detachment, 𝑓 d e t is the frequency of nucleation, and 𝑁 is the number of nucleation sites per unit area (nucleation site density).

The quenching heat flux has been analytically calculated by Del Valle and Kenning [5] as 𝑞 𝑄 = Θ 𝐷 𝑓 d e t 𝐴 2 𝜙 2 𝑘 𝑙 𝑇 𝑤 − 𝑇 𝑙 𝜋 Θ 𝐷 𝑘 𝑙 / 𝜌 𝑙 𝑐 𝑝 𝑙 , ( 2 3 ) where Θ 𝐷 is the waiting (or dwell) time elapsed between the detachment of a bubble and the nucleation of a subsequent one. The term 𝐴 2 𝜙 is the fraction of the wall area participating in the quenching heat flux.

Additional relationships needed to close the model include the following [6].

(i) Bubble diameter at departure [7] 𝑑 d e t − = 0 . 0 0 1 4 ⋅ E x p Δ 𝑇 s u b 4 5 , ( 2 4 ) where Δ 𝑇 s u b is in ∘ C , and 𝑑 d e t is in m. It is known that the liquid velocity has a significant effect on the bubble diameter at detachment.

(ii) Nucleation site density [8] 𝑁 𝑇 = 1 8 5 𝑤 − 𝑇 𝑙 1 . 8 0 5 . ( 2 5 )

Two very important terms in the model described above are the bubble waiting time and the nucleation frequency. A mechanistic approach to determine both parameters is described in Section based on the analysis of the bubble ebullition cycle.

4. Bubble Ebullition Cycle

The mechanisms of nucleation at a cavity on a heated wall and the subsequent bubble growth have been investigated very extensively before. However, most existing results are normally given in the form of relationships between selected parameters, rather than as complete analytical models. It turns out that using a rigorous analytical approach to combine transient heat transfer solutions for the heated wall and for the liquid filling the space vacated by departing bubbles, a consistent model can be derived for bubble ebullition in forced convection-subcooled boiling. A schematic of the ebullition cycle is shown in Figure . The cycle consists of two major stages: the quenching period and the bubble growth period.

The quenching (or dwell) period is initiated when the subcooled liquid fills the space near the heated wall vacated by the departing bubble.  During that time, the space inside a given cavity is gradually filled with vapor which eventually forms a hemispherical bubble on the top of the cavity. The quenching period is followed by a bubble growth period during which the bubble quickly expands forming a thin liquid sublayer separating the bubble from the heated surface. Since the liquid sublayer is very thin, the temperature drop across this layer is very small, so that the local wall temperature quickly drops to the saturated vapor temperature which itself only exceeds slightly the saturation temperature corresponding the system pressure.

Accounting to the periodic nature of the process, a closed-form solution can be obtained for the parameters characterizing the timing of bubble ebullition. The model is based on using coupled solutions to the transient heat conduction equations for the heated wall and the fluid laminar sublayer near the wall: 𝜕 𝑇 𝜕 𝑡 = 𝑎 ∇ 2 𝑇 , ( 2 5 ) where 𝑎 = 𝑘 / 𝜌 𝑐 𝑝 is the thermal diffusivity of the respective material.

Since the characteristic time of surface temperature fluctuations during nucleation is very short and, thus, the distance across the wall affected by a change in the surface temperature is small compared to the size of the surface area exposed to quenching by cold water, the wall heat conduction can be approximated by a one-dimensional model. Using a steady-state as a reference, the time-dependent temperature distribution across the heated wall during the quenching period (i.e., between 0 and Θ 𝐷 in Figure 2) becomes [9] 𝑇 𝑤 ( 𝑦 , 𝑡 ) = 𝑇 𝑜 + 𝑞 𝐻 𝑘 𝑤 𝑦 + 𝑎 𝑤 2 √ 𝜋 𝑡 0 𝑦 𝑎 𝑤 𝑡 − 𝑡 1 . 5 − 𝑦 E x p 2 4 𝑎 𝑤 𝑡 − 𝑡 × 𝑇 𝑖 𝑡 − 𝑇 𝑜 𝑑 𝑡 , ( 2 7 ) where 𝑞 𝐻 > 0 is the constant heating rate per unit surface area of the heated wall, 𝑇 𝑜 is the average (constant) temperature of the heated wall surface in contact with the fluid, and 𝑦 is the distance across the heater away from the wetted surface.

Using (27), the following expression can be derived for the time-dependent surface heat flux during the quenching (dwell) period 𝑞 𝑤 ( 𝑡 ) = − 𝑞 𝐻 + 𝑘 𝑤 √ 𝜋 𝑎 𝑤 𝑡 0 𝑑 𝑇 𝑖 𝑑 𝑡 𝑑 𝑡 √ 𝑡 − 𝑡 + 𝑇 𝑖 0 + − 𝑇 𝑜 √ 𝑡 , ( 2 8 ) where 𝑇 𝑖 ( 0 + ) is the wall temperature at the beginning of the dwell period (the “+” superscript is used to indicate that a temperature discontinuity occurs at 𝑡 = 0 , as shown in Figure 2).

Since the near-wall space vacated by the departing bubble is immediately filled by subcooled liquid, the liquid temperature starts quickly increasing in contact with the heated wall. The time- and position-dependent liquid temperatures across the laminar boundary layer (where the effect of flow-driven convection is negligible) can be determined by solving the transient heat conduction equation, (26). Assuming that the liquid temperature outside the boundary layer is constant, 𝑇 b l = c o n s t a n t , whereas the liquid temperature in contact with the heated wall is 𝑇 𝑖 ( 𝑡 ) , yields the following expression: 𝑡 𝑙 ( 𝑦 , 𝑡 ) = 𝑡 b l + 𝑎 𝑙 2 √ 𝜋 𝑡 0 𝑦 𝑎 𝑙 𝑡 − 𝑡 1 . 5 − 𝑦 e x p 2 4 𝑎 𝑙 𝑡 − 𝑡 𝑡 𝑖 𝑡 − 𝑡 b l 𝑑 𝑡 ( 2 9 )

In a manner similar to that used for the wall, the following expression is obtained from (29) for the instantaneous time-dependent surface heat flux into the fluid: 𝑞 𝑤 ( 𝑡 ) = − 𝑘 𝑤 𝜕 𝑇 𝑤 | | | 𝜕 𝑦 𝑦 = 0 𝑘 = − 𝑙 √ 𝜋 𝑎 𝑙 𝑡 0 𝑑 𝑇 𝑖 𝑑 𝑡 𝑑 𝑡 √ 𝑡 − 𝑡 + 𝑇 𝑖 0 + − 𝑇 b l √ 𝑡 . ( 3 0 )

Combining (28) and (30) and taking the limit as 𝑡 → 0 + yields the surface temperature at the beginning of the quenching (or dwell) period: 𝑇 𝑖 0 + = 𝑘 𝑤 𝑇 𝑜 √ 𝑎 𝑤 + 𝑘 𝑙 𝑇 b l √ 𝑎 𝑙 𝑘 𝑤 √ 𝑎 𝑤 + 𝑘 𝑙 √ 𝑎 𝑙 − 1 . ( 3 1 )

Substituting(31) back into the combined (28) and (30), one obtains 𝑡 0 𝑑 𝑇 𝑖 𝑡 𝑑 𝑡 𝑑 𝑡 √ 𝑡 − 𝑡 = 𝑞 𝐻 𝑘 𝑤 / √ 𝜋 𝑎 𝑤 + 𝑘 𝑙 / √ 𝜋 𝑎 𝑙 . ( 3 2 )

Solving (32) for 𝑇 𝑖 ( 𝑡 ) with the initial condition 𝑇 𝑖 ( 𝑡 ) = 𝑇 𝑖 ( 0 + ) yields 𝑇 𝑖 ( 𝑡 ) = 𝑇 𝑖 0 + + 2 𝑞 𝐻 √ 𝜋 𝑘 𝑤 / √ 𝑎 𝑤 + 𝑘 𝑙 / √ 𝑎 𝑙 √ 𝑡 . ( 3 3 )

Substituting (33) into (29) and rearranging yield the following expression for the instantaneous surface heat flux during the dwell period: 𝑞 𝑤 ( 𝑡 ) = − 𝑘 𝑙 𝜕 𝑇 𝑙 | | | 𝜕 𝑦 𝑦 = 0 = − 𝑘 𝑙 𝑇 𝑖 0 + − 𝑇 b l √ 𝜋 𝑎 𝑙 𝑡 − 𝑞 𝐻 𝑘 1 + 𝑤 / 𝑘 𝑙 √ 𝑎 𝑙 / 𝑎 𝑤 . ( 3 4 )

4.1. Dwell Time Θ 𝐷

The time of transition from the dwell period to the growth period is reached when the bubble radius becomes equal to the wall cavity radius. At this time instant, the steam temperature inside the hemispherical bubble outside the cavity must be equal to the temperature of the surrounding liquid at a distance from the wall equal to the cavity radius. Using the Clausius-Clapeyron equation to determine the steam superheat, and combining the resultant expression with (34), one arrives at the following relationship: 2 𝑞 𝐻 / √ 𝜋 𝑘 𝑤 / √ 𝑎 𝑤 + 𝑘 𝑙 / √ 𝑎 𝑙 √ Θ 𝐷 − 𝑇 𝑖 0 + − 𝑇 b l 𝑟 𝑐 √ 𝜋 𝑎 𝑙 Θ 𝐷 + 𝑇 𝑖 0 + − 𝑇 s a t − 𝑞 𝐻 𝑟 𝑐 / 𝑘 𝑙 𝑘 1 + 𝑤 / 𝑘 𝑙 √ 𝑎 𝑙 / 𝑎 𝑤 − 2 𝜎 𝑟 𝑐 𝑇 s a t 𝜐 𝑓 𝑔 ℎ 𝑓 𝑔 = 0 , ( 3 5 ) where 𝑟 𝑐 is the cavity radius and σis the surface tension.

As can be seen, for given fluid and wall properties and thermal conditions, (35) becomes an algebraic (quadratic) equation for the bubble dwell time, Θ 𝐷 .

Interestingly, if the dwell time is known, (35) can be used to determine the average surface heat flux during the quenching period. Specifically, integrating (34) from 0 to Θ 𝐷 yields | | 𝑞 𝐷 | | = 2 𝑘 𝑙 𝑇 𝑖 0 + − 𝑇 b l √ 𝜋 𝑎 𝑙 Θ 𝐷 + 𝑞 𝐻 𝑘 1 + 𝑤 / 𝑘 𝑙 √ 𝑎 𝑙 / 𝑎 𝑤 . ( 3 6 )

4.2. Growth Time Θ 𝐺

As soon as the bubble starts growing outside the cavity, the vapor is produced mainly from the evaporating liquid sublayer between the wall and the bubble. Since the temperature drop across the thin sublayer is very small, the local surface temperature during this period remains close to the saturation temperature. The energy balance for the growing bubble can be written as 𝜌 𝑔 ℎ 𝑓 𝑔 2 𝑑 𝑑 𝐵 = | | 𝑞 𝑑 𝑡 𝑤 | | ( 𝑡 ) . ( 3 7 )

The time-dependent wall heat flux during the bubble growth time can be obtained from (28) by noticing that during this period, that is, for Θ 𝐷 < 𝑡 < Θ 𝐷 + Θ 𝐺 , the wall temperature remains approximately constant and equal to the saturation temperature: 𝑞 𝑤 ( 𝑡 ) = − 𝑞 𝐻 − 𝑘 𝑤 𝑇 𝑖 Θ 𝐷 − 𝑇 s a t 𝜋 𝑎 𝑤 𝑡 − Θ 𝐷 . ( 3 8 )

Substituting (38) into (37) and integrating between Θ 𝐷 and Θ 𝐷 + Θ 𝐺 yield the following expression for the maximum bubble diameter at detachment: 𝑑 d e t = 2 𝜌 𝑣 ℎ 𝑓 𝑔 Θ 𝐷 + Θ 𝐺 Θ 𝐷 | | 𝑞 𝑤 | | = 2 ( 𝑡 ) 𝑑 𝑡 𝜌 𝑣 ℎ 𝑓 𝑔 𝑞 𝐻 Θ 𝐺 + 2 𝑘 𝑤 𝑇 𝑖 Θ 𝐷 − 𝑇 s a t √ 𝜋 𝑎 𝑤 √ Θ 𝐺 , ( 3 9 ) where (see (33)) 𝑇 𝑖 Θ 𝐷 = 𝑇 𝑖 0 + + 2 𝑞 𝐻 / √ 𝜋 𝑘 𝑤 / √ 𝑎 𝑤 + 𝑘 𝑙 / √ 𝑎 𝑙 √ Θ 𝐷 . ( 4 0 )

The bubble diameter at detachment can be evaluated using one of several different models that have been developed to date. In particular, using the force balance for a single bubble, Staub [10] developed the following expression: 𝑑 d e t ⎧ ⎪ ⎨ ⎪ ⎩ ⎷ = 0 . 5 3 𝜏 𝑤 2 𝜌 𝑙 − 𝜌 𝑣 𝑔 2 + 𝛽 3 𝜎 𝐹 𝑐 𝜌 𝑙 − 𝜌 𝑣 𝑔 − 3 𝜏 𝑤 2 𝜌 𝑙 − 𝜌 𝑣 𝑔 ⎫ ⎪ ⎬ ⎪ ⎭ , ( 4 1 ) where 𝜏 𝑤 is the wall shear stress and 𝐹 ( 𝛽 𝑐 ) is an experimentally determined function of the contact angle.

For given bubble diameter at detachment and bubble dwell time, (39) can be readily solved for the bubble growth time Θ 𝐺 .

Interestingly, if the bubble growth time is known, (39) can be used again, this time to determine the average wall heat flux during the bubble growth period: | | 𝑞 𝐺 | | = 1 Θ 𝐺 Θ 𝐷 + Θ 𝐺 Θ 𝐷 | | 𝑞 𝑤 | | ( 𝑡 ) 𝑑 𝑡 = 𝑞 𝐻 + 2 𝑘 𝑤 𝑇 𝑖 Θ 𝐷 − 𝑇 s a t √ 𝜋 𝑎 𝑤 Θ 𝐺 . ( 4 2 )

Finally, the bubble frequency of detachment and the total detachment time can be obtained from 𝑓 d e t = 1 Θ d e t = 1 Θ 𝐷 + Θ 𝐺 . ( 4 3 )

Naturally, summing up (34) and (42) yields | | 𝑞 𝐷 | | + | | 𝑞 𝐺 | | = 𝑞 𝐻 , ( 4 4 ) where 𝑞 𝐻 is the surface heat flux of the heater (see (27)).

5. Wall Temperature Excursion (CHF)

5.1. Physical Concept

The wall temperature excursion (or CHF) in low-quality boiling is associated with the ability of the bubbles formed at the nucleation sites to depart from the wall, so that the vacated space can be filled with fresh liquid, so that the quenching and bubble growth processes can continue [9, 11]. Thus, two conditions must be satisfied simultaneously to avoid wall temperature excursion, one concerned with dispersed bubble concentration in the near-wall region, the other with the velocity at which the bubble departs from the wall to make room for the next generation bubble formed at the same nucleation site. Since the size of bubbles departing from the nucleation sites is normally very small (of the order of 1 mm or less), their shape is nearly spherical. Thus, the two conditions mentioned above can be formulated as follows.

(1) The bubble maximum void fraction of dispersed bubbles in the near-wall region cannot exceed the maximum packing factor for spherical particles, the theoretical value of which (obtained from geometrical considerations) is 𝛼 ≤ 𝛼 m a x = 𝜋 3 √ 2 ≈ 0 . 7 4 . ( 4 5 )

(2) Taking into account that to form a new undeformed bubble at a given cavity at the time of detachment, the distance between the previously formed bubble at the same cavity and the heated wall at the same time instant must be at least equal to the bubble diameter, the maximum evaporation heat flux to avoid temperature excursion must satisfy the following condition: 𝑞 𝑒 ℎ 𝑓 𝑔 = 𝜋 6 𝑑 3 d e t 𝜌 𝑣 𝑓 d e t 𝑁 ≤ 𝜌 𝑣 𝑢 d e t , m i n 𝛼 , ( 4 6 ) where 𝑢 d e t , m i n = 𝑚 𝑑 d e t 𝑓 d e t . ( 4 7 )

The coefficient, 𝑚 ≥ 1 , reflects the fact that in reality the axial motion of bubbles may require a distance from the wall larger than one bubble diameter to avoid the coalescence with other bubbles formed in the adjacent upstream cavities.

It is interesting to notice that (46) can be rewritten as 𝜋 6 𝑑 2 d e t 𝑁 ≤ 𝑚 𝛼 . ( 4 8 )

5.2. Near-Wall Heat Transfer in the Presence of Elongated Bubbles

A schematic illustrating the near-wall conditions in the presence of elongated bubbles is shown in Figure 3. It follows from the previous discussion that at any axial location along the channel the time-dependent wall temperature fluctuates periodically, increasing in the poor heat transfer (elongated bubble) region and decreasing, where heat transfer is more efficient (in the nucleate boiling region).

387020.fig.003
Figure 3: Phase distribution in the near-wall region.

According to Figure 3, averaging the local wall heat flux over both regions always yields the given constant heating rate of the heater per unit surface area. In particular, the overall average heat flux can be partitioned into the following terms: 𝑞 𝐻 = 𝑞 𝑣 1 + 𝑞 𝑣 2 . ( 4 9 )

In (49), 𝑞 𝑣 1 is the nucleate boiling component corresponding to small dispersed bubbles: 𝑞 𝑣 1 = 𝑞 N B 𝐴 𝑣 1 , ( 5 0 ) where 𝑞 N B is given by (20), and 𝐴 𝑣 1 is the area density (wall area fraction) for the small bubble region.

The heat flux component corresponding to elongated bubbles, 𝑞 𝑣 2 , is given by[9, 11] 𝑞 𝑣 2 = 𝑞 l s 𝐴 L B 𝐿 l s 𝐿 L B + 𝑞 d r y 𝐴 L B 𝐿 1 − l s 𝐿 L B , ( 5 1 ) where 𝐴 L B is the wall area density for elongated bubbles, 𝐿 l s is the length of the liquid sublayer underneath an elongated bubble, 𝑞 l s is the corresponding average heat flux across the liquid layer, and 𝑞 d r y is the heat transfer rate per unit wall in the dry region.

The length of large bubbles has been measured by Gersey and Mudawar [12] and found to agree very well with the critical wavelength of the Helmholtz instability vapor/liquid interface. This length can be calculated from 𝐿 L B = 2 𝜋 𝜎 𝜌 𝑙 𝑣 𝑣 − 𝑣 𝑙 2 𝜌 1 + 𝑙 𝜌 𝑣 . ( 5 2 )

Another parameter characterizing elongated bubbles is their distance from the wall. Assuming that the initial distance (at the tip of the bubbles) corresponds to the viscous sublayer thickness [13], we obtain 𝐿 L B = 1 0 𝑣 𝑙 𝜌 𝑙 𝜏 𝑤 , ( 5 3 ) where 𝜏 𝑤 is the wall shear stress.

5.3. Wall Temperature Excursion Criterion

It can be readily noticed that as the total wall heat flux increases, the elongated bubble area density will also increase and so will the dry region under the elongated bubbles. Consequently, the local heat flux in the elongated bubble region, 𝑞 𝑣 2 , will go down and that in the dispersed bubble region, 𝑞 𝑣 1 , will go up. With the corresponding area density decreasing, nucleate boiling in the region between elongated bubbles will intensify, eventually leading to a situation when replenishment of the near-wall region with liquid will be no longer possible. As a result, a sudden wall temperature excursion will occur. The critical condition beyond which small bubbles cannot be removed away from the wall and replaced by fresh liquid can be written as 𝜋 6 𝑑 2 d e t 𝑁 ≤ 𝑚 𝛼 N B 𝛼 = 𝑚 𝑣 1 𝐴 𝑣 1 , ( 5 4 ) where 𝛼 N B is the volumetric fraction of dispersed bubbles in the nucleate boiling region and 𝛼 𝑣 1 is the local near-wall volume fraction (per unit volume of the channel) of the dispersed bubbles.

Naturally, the local near-wall void fraction is the sum of the partial vapor volume fractions of the dispersed bubbles and the elongated bubbles: 𝛼 = 𝛼 𝑣 1 + 𝛼 𝑣 2 . ( 5 5 )

6. Model Testing and Validation

6.1. Nucleate Boiling

The theoretical model of bubble ebullition discussed in Sections 3 and 4 was used to evaluate various parameters characterizing nucleate boiling, and compare the result of predictions against the results of measurements. Typical results are shown in Tables 1 and 2 and in Figure 4.

tab1
Table 1: Predicted bubble diameter at detachment.
tab2
Table 2: Predicted nucleation frequency.
387020.fig.004
Figure 4: A comparison between the predicted and experimental evaporation heat flux [14].

The predictions for the bubble diameter at departure and the nucleation frequency are shown in Tables 1 and 2. First, the predicted bubble diameter has been compared against the expression given by (24), obtained for a liquid velocity of 0.2 m/s. The calculations using the present models have been performed for a wall heat flux of 𝑞 𝑤 = 4 . 7 3 ⋅ 1 0 3 k W / m 2 , a subcooling of 1 0 ∘ C , and three different pressures. The results are shown in Table 1 [9].

As can be seen, whereas the predicted bubble diameter at detachment compares well against the simple correlation of Tolubinsky and Kostanchuk [7], the current model also allows for quantifying the effect of system pressure (and other parameters, such as velocity) which is ignored in the expression given by (24).

Next, the calculated bubble departure frequency was compared against the correlation of Ceumern-Lindenstjerna [15]. The calculations were performed for 𝑞 𝑤 = 1 . 4 ⋅ 1 0 3 k W / m 2 . The results for three different pressures are shown in Table 2.

The predicted evaporation heat flux in a boiling channel used by Bartolomei and Chanturia [13] is shown in Figure 4. In this comparison, the nucleation sites density was obtained from the expression proposed by Lemmert and Chawla [8] (see (25)). The measured wall superheat was about 9 . 5 ∘ C but changed slightly along the channel. In order to account for the uncertainties in the wall temperature measurements, the results for two fixed values are shown in Figure 4. As can be seen, a very good agreement has been obtained, especially for liquid subcooling between 5 ∘ C and 1 5 ∘ C .

It can be noticed that the current model is capable of quantifying the effect of various physical parameters which are normally not accounted for when a phenomenological approach, based solely on correlating experimental data, is used.

6.2. Temperature Excursion

To perform predictions of the temperature excursion, the three-dimensional three-field model of two-phase flow discussed in Section 2 has been combined with both models discussed in Sections 3–5. The calculations have been performed for the experimental conditions of Hino and Ueda [16], in which R113 was used at a pressure of 147 kPa. The heated test section was 0.357 m long and 0.018 m ID tube, with a centrally located heated rod, 0.008 m diameter. The outer tube wall was insulated, and there was an unheated section installed upstream of the annulus-shaped heated section, allowing the flow to reach fully developed conditions at the entrance to the heater.

Typical radial and axial distributions of various local flow parameters are shown in Figure 5. The axial distributions of the near-wall concentrations of both small bubbles and elongated bubbles are shown in Figure 6. As it can be seen, the rate of bubble concentration increase in the dispersed-bubble region accelerates as the elongated bubbles start being formed. This is due to a dramatic reduction in the heat transfer removal rate in the elongated bubble region, which in turn increases the local heat flux in the nucleate-boiling region. As it can be seen in Figure 6, whereas the average wall heat flux is fixed at a level of 2 4 1 k W / m 2 , the nucleate boiling heat flux experiences a dramatic growth as the concentration of elongated bubbles stars increasing. Indeed, due to the combined effects of increasing evaporation rate and increasing fraction of the wall area occupied by elongated bubbles, the local void fraction approached the value, 𝛼 = 0 . 7 5 , which already exceeded the limit given by (44). This, in turn, stopped the wall replenishment by liquid phase and caused a sudden temperature excursion. Thus, one concludes that the assumed wall heat flux was just beyond the onset of CHF. Converting the difference between the actual volume fraction value of 0.75 and the critical value of 0.74 into the corresponding power level difference yields the critical heat flux (CHF) of about 2 3 8 k W / m 2 . Since the uniform heat flux of the heater used in the calculations ( 2 4 1 k W / m 2 ) corresponded to the experimental onset of temperature excursion, the estimated prediction error in the present case was less than 2%. Similar calculations performed for other conditions produced errors of the order of ± 1 5 %.

387020.fig.005
Figure 5: Color contours of the calculated phasic velocities, temperature, and void fraction in subcooled boiling: m a s s fl u x = 5 1 2 k g / m 2 s , h e a t fl u x = 2 4 1 k W / m 2 , inlet subcooling = 30 K.
fig6
Figure 5: Axial distributions of near-wall channel parameters for the conditions shown in Figure 5: (a) the concentrations of small bubbles and elongated bubbles, (b) a comparison between the heat transfer removal rate in the elongated bubble region, average wall heat flux, and critical heat flux.

Figures 7 and 8 show the calculated CHF for various flow conditions, in comparison with typical data based on experimental correlations. The predicted critical heat flux is within the range of the measured values for heated channels operating at subcooled boiling or low-quality boiling conditions. What is particularly important is that the predicted trends in the critical heat flux agree well with the existing experimental evidence. Specifically, 𝑞 𝑐 gradually decreases with increasing velocity (and, thus, flow rate) and decreases with increasing vapor concentration (flow quality).

387020.fig.007
Figure 7: Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].
387020.fig.008
Figure 8: Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].

7. Conclusions

Several aspects of mechanistic multidimensional modeling and computer simulations of two-phase flows and boiling heat transfer have been discussed. The specific models included the mechanisms of local-subcooled boiling heat transfer in forced convection flows, a mechanistic approach to bubble ebullition cycle, and criteria for temperature excursion (CHF) in low-quality flows.

The results of model testing indicate that the proposed physical mechanisms are consistent with the existing experimental evidence regarding the phase and temperature distributions, and wall temperature excursion. The current approach is particularly suitable for implementation in general CFD computational models using a multifield concept of two-phase flow.

References

  1. M. Z. Podowski, “On the mechanistic modeling of interfacial phenomena in gas/liquid two-phase flows,” in Proceedings of the 6th International Conference on Multiphase Flow (ICMF '07), Leipzig, Germany, July 2007.
  2. B. E. Launder and D. B. Spalding, “The numerical computation of turbulent flows,” Computer Methods in Applied Mechanics and Engineering, vol. 3, no. 2, pp. 269–289, 1974. View at Publisher · View at Google Scholar
  3. Y. Sato, M. Sadatomi, and K. Sekoguchi, “Momentum and heat transfer in two-phase bubble flow—I. Theory,” International Journal of Multiphase Flow, vol. 7, no. 2, pp. 167–177, 1981. View at Publisher · View at Google Scholar
  4. N. Kurul and M. Z. Podowski, “On the modeling of multidimensional effects in boiling channels,” in Proceedings of the 27th National Heat Transfer Conference, Minneapolis, Minn, USA, July 1991.
  5. V. H. Del Valle and D. B. R. Kenning, “Subcooled flow boiling at high heat flux,” International Journal of Heat and Mass Transfer, vol. 28, no. 10, pp. 1907–1920, 1985. View at Publisher · View at Google Scholar
  6. A. Alajbegovic, N. Kurul, M. Z. Podowski, D. Drew, and R. T. Lahey, Jr., “A new mechanistic model of critical heat flux in forced-convection subcooled boiling,” in Proceedings of the 8th International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-8), pp. 973–980, Kyoto, Japan, September-October 1997.
  7. V. I. Tolubinsky and D. M. Kostanchuk, “Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling,” in Proceedings of the 4th International Heat Transfer Conference, vol. 5, Paris, France, August 1970, paper no. B-2.8.
  8. M. Lemmert and L. M. Chawla, “Influence of flow velocity on surface boiling heat transfer coefficient,” in Heat Transfer in Boiling, E. Hahne and U. Grigull, Eds., Academic Press and Hemisphere, New York, NY, USA, 1977.
  9. R. M. Podowski, D. A. Drew, R. T. Lahey, Jr., and M. Z. Podowski, “A mechanistic model of the ebullition cycle in forced-convection subcooled boiling,” in Proceedings of the 8th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-8), vol. 3, pp. 1530–1537, Kyoto, Japan, September-October 1997.
  10. F. W. Staub, “The void fraction in subcooled boiling: prediction of the initial point of net vapor generation,” Journal of Heat Transfer, vol. 90, pp. 151–157, 1968.
  11. M. Z. Podowski and S. P. Antal, “CFD predictions of temperature excursion (CHF) in low-quality boiling in heated channels,” in Proceedings of the 12th International Heat Transfer Conference, Grenoble, France, August 2002.
  12. C. O. Gersey and I. Mudawar, “Effects of heater length and orientation on the trigger mechanism for near-saturated flow boiling critical heat flux—I. Photographic study and statistical characterization of the near-wall interfacial features,” International Journal of Heat and Mass Transfer, vol. 38, no. 4, pp. 629–641, 1995. View at Publisher · View at Google Scholar
  13. E. Moursali, J. L. Marié, and J. Bataille, “An upward turbulent bubbly boundary layer along a vertical flat plate,” International Journal of Multiphase Flow, vol. 21, no. 1, pp. 107–117, 1995. View at Publisher · View at Google Scholar
  14. G. G. Bartolomei and V. M. Chanturia, “Experimental study of true void fraction when boiling subcooled water in vertical tubes,” Thermal Engineering, vol. 14, pp. 123–128, 1967.
  15. W.-C. B. Ceumern-Lindenstjerna, “Bubble departure diameter and release frequencies during nucleate pool boiling of water and aqueous NaCl solution,” in Heat Transfer in Boiling, E. Hahne and U. Grigull, Eds., Academic Press and Hemisphere, New York, NY, USA, 1977.
  16. R. Hino and T. Ueda, “Studies on heat transfer and flow characteristics in subcooled flow boiling—part 2: flow characteristics,” International Journal of Multiphase Flow, vol. 11, no. 3, pp. 283–297, 1985. View at Publisher · View at Google Scholar
  17. E. Janssen and S. Levy, “Burnout limit curves for boiling water reactors,” General Electric, San Jose, Calif, USA, 1962.