Mechanistic Multidimensional Modeling of Forced Convection Boiling Heat Transfer

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𝜕(𝛼_𝑘𝜌_𝑘𝐮_𝑘)𝜕𝑡+∇⋅(𝛼_𝑘𝜌_𝑘𝐮_𝑘𝐮_𝑘)=−∇(𝛼_𝑘𝑝_𝑘)+∇⋅(𝛼_𝑘𝝉^tot_𝑘)+𝛼_𝑘𝜌_𝑘𝐠+𝐌^𝑖_𝑘, (2)

Energy𝜕(𝛼_𝑘𝜌_𝑘𝑒_𝑘)𝜕𝑡+∇⋅(𝛼_𝑘𝜌_𝑘𝐮_𝑘𝑒_𝑘)=−∇⋅(𝛼_𝑘𝐪^tot_𝑘)+∇⋅[𝛼_𝑘(−𝑝_𝑘𝐈+𝝉^tot_𝑘)⋅𝐮_𝑘]+𝛼_𝑘𝜌_𝑘𝐠⋅𝐮_𝑘+𝐸^𝑖_𝑘, (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-𝑘, 𝝉^tot_𝑘 is the total shear stress, 𝐌^𝑖_𝑘 is the total interfacial force, 𝐪^tot_𝑘 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.

(a) Two-field model

(b) Three-field model

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-𝑣₁), it is appropriate to use a separate field for the former, denoted here as field-𝑣₂.

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] 𝜇^𝑡_𝐵=𝐶_𝜇𝐵𝛼_𝑙𝜌_𝑙²_𝑘=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 Γ_𝑙=−Γ^evap+Γ^cond=−(Γ_𝑣1+Γ_𝑣2Γ),_𝑣1=Γ^evap_𝑣1−Γ^cond_𝑣1−Γ^clsc_𝑣1,Γ_𝑣2=Γ^evap_𝑣2−Γ^cond_𝑣2+Γ^clsc_𝑣1. (7)(8)(9)

In (7), Γ^evap=Γ^evap_𝑣1+Γ^evap_𝑣2 and Γ^cond=Γ^cond_𝑣1+Γ^cond_𝑣2 are the total evaporation and condensation rates per unit volume, respectively, and Γ^clsc_𝑣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: 𝐌^𝑖_𝑘=𝑝^𝑖_𝑘∇𝛼_𝑘−𝝉^𝑖_𝑘⋅∇𝛼_𝑘+Γ_𝑘𝐮^𝑖_𝑘+_𝑗𝐌^𝑖_𝑘,𝑗. (10)

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

Drag Force𝐌^𝐷_{𝑙,𝑣𝑘}=−𝐌^𝐷_{𝑣𝑘,𝑙}1=−8𝐶_𝐷,𝑘𝐴_𝑣^′′′_𝑘𝜌_𝑙(𝐮_𝑣𝑘−𝐮_𝑙)|𝐮_𝑣𝑘−𝐮_𝑙|, (11)

Virtual Mass Force𝐌^VM_{𝑙,𝑣𝑘}=−𝐌^VM_{𝑣𝑘,𝑙}=−𝐶_VM,𝑘𝛼_𝑣𝑘𝜌_𝑙(𝐷_𝑣𝑘𝐮_𝑣𝑘−𝐷𝐷𝑡^𝑙𝐮_𝑙𝐷𝑡), (12)

Lift Force𝐌^𝐿_{𝑙,𝑣𝑘}=−𝐌^𝐿_{𝑣𝑘,𝑙}=−𝐶_𝐿,𝑘𝛼_𝑣𝑘𝜌_𝑙(𝐮_𝑣𝑘−𝐮_𝑙)×∇×𝐮_𝑙, (13)

Turbulent Dispersion Force𝐌^TD_{𝑙,𝑣𝑘}=−𝐌^TD_{𝑣𝑘,𝑙}=𝐶_TD,𝑘𝛼_𝑣𝑘𝜌_𝑙∇(𝛼_𝑙𝜅_𝑙) (14)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 Γ^cond_𝑘=𝑞_𝑘^′′𝑖ℎ_𝑓𝑔(𝑘=𝑣₁,𝑣₂), (15) where the interfacial energy transport from vapor field-k to the continuous liquid field is given by 𝑞_𝑘^′′𝑖=𝐻^𝑖_𝑘−𝑙(𝑇_sat−𝑇_𝑙)𝐴_𝑘^′′′(𝑘=𝑣₁,𝑣₂). (16)

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.33],(𝑘=𝑣₁,𝑣₂). (17)

Equations (16) and (17) imply 𝑞_𝑙^′′𝑖=𝑞_𝑣^′′𝑖₁+𝑞_𝑣^′′𝑖₂=[𝐻^𝑖_𝑣1−𝑙𝐴_𝑣^′′′₁+𝐻^𝑖_𝑣2−𝑙𝐴_𝑣^′′′₂](𝑇_sat−𝑇_𝑙). (18)

The rate of bubble coalescence can be obtained from Γ^clsc_𝑣1=𝑏_{𝑣1−𝑣1}𝛼²_𝑣1+𝑏_{𝑣1−𝑣2}𝛼_𝑣1𝛼_𝑣2, (19)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]: 𝑞^′′_NB=𝑞^′′_1𝜙+𝑞_𝑒^′′+𝑞_𝑄^′′, (20) 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𝜙St_𝑃𝜌_𝑙𝑐_𝑝𝑙𝑢_𝑙,𝑃(𝑇_𝑤−𝑇_𝑙,𝑃), (21) where 𝐴^′′_1𝜙 is the fraction of the wall unaffected by the nucleation sites, St_𝑃 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𝑑³_det𝜌_𝑣𝑓_det𝑁^′′ℎ_𝑓𝑔, (22) where 𝑑_det is the critical bubble diameter at detachment, 𝑓_det 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 𝑞_𝑄^′′=Θ_𝐷𝑓_det𝐴^′′_2𝜙2𝑘_𝑙(𝑇_𝑤−𝑇_𝑙)𝜋Θ_𝐷𝑘_𝑙/(𝜌_𝑙𝑐_𝑝𝑙), (23) 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] 𝑑_det=0.0014⋅Exp(−Δ𝑇_sub45), (24) where Δ𝑇_sub is in ^∘C, and 𝑑_det 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] 𝑁^′′=185(𝑇_𝑤−𝑇_𝑙)^1.805. (25)

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 4 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 2. The cycle consists of two major stages: the quenching period and the bubble growth period.

Figure 2 

Surface temperature oscillations in the wall during the ebullition cycle.

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: 𝜕𝑇𝜕𝑡=𝑎∇²𝑇, (26) 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√𝜋^𝑡₀𝑦[𝑎_𝑤(𝑡−𝑡^′)]^1.5𝑦Exp[−²4𝑎_𝑤(𝑡−𝑡^′)]×[𝑇_𝑖(𝑡^′)−𝑇_𝑜]𝑑𝑡^′, (27) 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⁺)−𝑇_𝑜√𝑡]}, (28) 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, 𝑇_bl=constant, whereas the liquid temperature in contact with the heated wall is 𝑇_𝑖(𝑡), yields the following expression: 𝑡_𝑙(𝑦,𝑡)=𝑡_bl+𝑎_𝑙2√𝜋^𝑡₀𝑦[𝑎_𝑙(𝑡−𝑡^′)]^1.5𝑦exp[−²4𝑎_𝑙(𝑡−𝑡^′)][𝑡_𝑖(𝑡^′)−𝑡_bl]𝑑𝑡^′ (29)

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⁺)−𝑇_bl√𝑡}. (30)

Combining (28) and (30) and taking the limit as 𝑡→0⁺ yields the surface temperature at the beginning of the quenching (or dwell) period: 𝑇_𝑖(0⁺𝑘)=(^𝑤𝑇_𝑜√𝑎_𝑤+𝑘_𝑙𝑇_bl√𝑎_𝑙𝑘)(^𝑤√𝑎_𝑤+𝑘_𝑙√𝑎_𝑙)⁻¹. (31)

Substituting(31) back into the combined (28) and (30), one obtains ^𝑡₀𝑑𝑇_𝑖(𝑡^′)𝑑𝑡^′𝑑𝑡^′𝑡−𝑡^′=𝑞_𝐻^′′𝑘_𝑤/√𝜋𝑎_𝑤+𝑘_𝑙/√𝜋𝑎_𝑙. (32)

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

Substituting (33) into (29) and rearranging yield the following expression for the instantaneous surface heat flux during the dwell period: 𝑞_𝑤^′′(𝑡)=−𝑘_𝑙𝜕𝑇_𝑙|𝜕𝑦_𝑦=0=−𝑘_𝑙𝑇_𝑖(0⁺)−𝑇_bl𝜋𝑎_𝑙𝑡−𝑞_𝐻^′′1+(𝑘_𝑤/𝑘_𝑙)𝑎_𝑙/𝑎_𝑤. (34)

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⁺)−𝑇_bl]𝑟_𝑐𝜋𝑎_𝑙Θ_𝐷+𝑇_𝑖(0⁺)−𝑇_sat−𝑞_𝐻^′′𝑟_𝑐/𝑘_𝑙1+(𝑘_𝑤/𝑘_𝑙)𝑎_𝑙/𝑎_𝑤−2𝜎𝑟_𝑐𝑇_sat𝜐_𝑓𝑔ℎ_𝑓𝑔=0, (35) 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⁺)−𝑇_bl]𝜋𝑎_𝑙Θ_𝐷+𝑞_𝐻^′′1+(𝑘_𝑤/𝑘_𝑙)𝑎_𝑙/𝑎_𝑤. (36)

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𝑑(𝑑_𝐵)𝑑𝑡=|𝑞_𝑤^′′(𝑡)|. (37)

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: 𝑞_𝑤^′′(𝑡)=−𝑞_𝐻^′′−𝑘_𝑤𝑇_𝑖(Θ_𝐷)−𝑇_sat𝜋𝑎_𝑤(𝑡−Θ_𝐷). (38)

Substituting (38) into (37) and integrating between Θ_𝐷 and Θ_𝐷+Θ_𝐺 yield the following expression for the maximum bubble diameter at detachment: 𝑑_det=2𝜌_𝑣ℎ_𝑓𝑔^{Θ𝐷+Θ𝐺}_Θ𝐷|𝑞_𝑤^′′=2(𝑡)|𝑑𝑡𝜌_𝑣ℎ_𝑓𝑔[𝑞_𝐻^′′Θ_𝐺+2𝑘_𝑤𝑇_𝑖(Θ_𝐷)−𝑇_sat√𝜋𝑎_𝑤Θ_𝐺], (39) where (see (33)) 𝑇_𝑖(Θ_𝐷)=𝑇_𝑖(0⁺)+2𝑞_𝐻^′′/√𝜋𝑘_𝑤/√𝑎_𝑤+𝑘_𝑙/√𝑎_𝑙Θ_𝐷. (40)

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: 𝑑_det=0.5{[3𝜏_𝑤2(𝜌_𝑙−𝜌_𝑣])𝑔²+3𝜎𝐹(𝛽_𝑐)(𝜌_𝑙−𝜌_𝑣−)𝑔3𝜏_𝑤2(𝜌_𝑙−𝜌_𝑣)𝑔}, (41) 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𝑘_𝑤𝑇_𝑖(Θ_𝐷)−𝑇_sat𝜋𝑎_𝑤Θ_𝐺. (42)

Finally, the bubble frequency of detachment and the total detachment time can be obtained from 𝑓_det=1Θ_det=1Θ_𝐷+Θ_𝐺. (43)

Naturally, summing up (34) and (42) yields |𝑞_𝐷^′′|+|𝑞_𝐺^′′|=𝑞_𝐻^′′, (44) 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 𝛼≤𝛼_max=𝜋3√2≈0.74. (45)

(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𝑑³_det𝜌_𝑣𝑓_det𝑁^′′≤𝜌_𝑣𝑢_det,min𝛼, (46) where 𝑢_det,min=𝑚𝑑_det𝑓_det. (47)

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𝑑²_det𝑁^′′≤𝑚𝛼. (48)

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).

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: 𝑞_𝐻^′′=𝑞_𝑣^′′₁+𝑞_𝑣^′′₂. (49)

In (49), 𝑞_𝑣^′′₁ is the nucleate boiling component corresponding to small dispersed bubbles: 𝑞_𝑣^′′₁=𝑞^′′_NB𝐴_𝑣^′′₁, (50) where 𝑞^′′_NB is given by (20), and 𝐴_𝑣^′′₁ is the area density (wall area fraction) for the small bubble region.

The heat flux component corresponding to elongated bubbles, 𝑞_𝑣^′′₂, is given by[9, 11] 𝑞_𝑣^′′₂=𝑞^′′_ls𝐴^′′_LB𝐿_ls𝐿_LB+𝑞^′′_dry𝐴^′′_LB𝐿(1−^ls𝐿_LB), (51) where 𝐴^′′_LB is the wall area density for elongated bubbles, 𝐿_ls is the length of the liquid sublayer underneath an elongated bubble, 𝑞^′′_ls is the corresponding average heat flux across the liquid layer, and 𝑞^′′_dry 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 𝐿_LB=2𝜋𝜎𝜌_𝑙(𝑣_𝑣−𝑣_𝑙)²𝜌(1+^𝑙𝜌_𝑣). (52)

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 𝐿_LB=10𝑣_𝑙𝜌_𝑙𝜏_𝑤, (53) 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, 𝑞_𝑣^′′₂, will go down and that in the dispersed bubble region, 𝑞_𝑣^′′₁, 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𝑑²_det𝑁^′′≤𝑚𝛼_NB𝛼=𝑚^𝑣1𝐴_𝑣^′′₁, (54) where 𝛼_NB 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. (55)

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.

Table 1 

Predicted bubble diameter at detachment.

Table 2 

Predicted nucleation frequency.

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.73⋅10³kW/m², a subcooling of 10^∘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⋅10³kW/m². 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 15^∘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 241kW/m², 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.75, 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 238kW/m². Since the uniform heat flux of the heater used in the calculations (241kW/m²) 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 ±15%.

Figure 5 

Color contours of the calculated phasic velocities, temperature, and void fraction in subcooled boiling: massflux=512kg/m2s, heatflux=241kW/m2, inlet subcooling = 30 K.

(a)

(b)

Figure 6 

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).

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].

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.