Mechanistic Prediction of Nucleate Boiling Heat Transfer–Achievable or a Hopeless Task?

Numerous studies of different modes of boiling have been reported in the literature during the last half of the twentieth century. Of the three modes of boiling—film, nucleate, and transition boiling—film boiling is most amenable to analysis. Nucleate and transition boiling are most complex, and invariably empirical correlations or correlations having a mechanistic basis have been proposed. The correlations serve a useful purpose in their application to engineered systems. However, because of their limited range of applicability, they can rarely be applied with confidence to new situations. In addition, many of these correlations for the global process, are rarely consistent at the subprocess level. Mechanistic models based on basic principles can alleviate this problem. For this purpose, effort has continued to be made by researchers to develop mechanistic models. In this work, we first review the past work on nucleate boiling heat transfer and thereafter provide results from a new approach that is based on direct numerical simulation of the process.

Past studies of nucleate boiling have resulted in either empirical correlations for the dependence of wall heat flux on wall superheat and other fluid and solid surface properties or correlations that are based on modeling of the underlying subphenomena.

It should be noted that all of the correlations suggest that wall heat flux varies as $ΔTw3$⁠. Rohsenow’s correlation suggests heat flux to vary as the square root of gravity, whereas the other two correlations indicate heat flux to be virtually independent of gravity. Equation 5 accounts for surface roughness, but does not account for the variation in the degree of surface wettability. Rohsenow’s correlation implicitly accounts for surface wettability through the proportionality constant $Cs$⁠, whereas the generalized correlation of Stephan and Abdelsalam bounds the data for a wide range of contact angles. Similarly, while Eq. 5 accounts for heater geometry, the other two correlations are independent of heater geometry.

It should be noted that none of the above correlations explicitly account for surface wettability. Wang and Dhir (11) systematically studied the effect of surface wettability by using the same surface-liquid (polished copper surface and water) combination, while controlling the degree of oxidation of the surface. They found the number density of active sites for a given cavity diameter to decrease by a factor of about 25 when the contact angle was reduced from $90to18deg$⁠. Subsequently, Basu et al. (12) have correlated nucleation site density data for a variety of liquids and contact angles, for both pool and flow boiling, as

It is generally believed that heterogeneous nucleation on a surface occurs on defects or cavities that trap gas/vapor. Bankoff (13) was the first to propose that an advancing liquid front will trap gas in a wedge-shaped cavity provided the advancing contact angle was greater than the wedge angle. In their study, Wang and Dhir (14) found that any conical or wedge-shaped cavities on the surface were too shallow to trap any gas. They noted that cavities that trapped gas were spherical in nature. By minimizing the Helmholtz free energy for a liquid/gas/solid system, they developed a criterion for gas entrapment according to which a cavity will trap gas only if the contact angle was greater than the cavity mouth angle. A key weakness in their approach is that it excluded effect of gas diffusion into liquid. Recently, Pesse et al. (15) experimentally studied filling of rectangular microchannels closed at one end. They found that the rate of filling of the channels was dependent on surface forces and the diffusion of gas into liquid. For low contact angles, multiple interfaces were formed in the channel. It was found that given sufficient time, the microchannels were filled with liquid for all contact angles. Thus, time for which surface remains in contact with liquid is another variable.

Despite numerous efforts, we do not yet have a model/correlation that can be used to predict a priori the number of cavities that will become active at a particular superheat for a given liquid-solid combination. There are several impediments in our ability to have a comprehensive predictive model for density of active nucleation sites. These include the degree to which the surface must be characterized with respect to the number of cavities that are present on the surface, as well as their shape and size. This can be an extremely tedious task. The degree of surface wettability greatly influences how rapidly the cavities are filled with liquid. As such, it is important to know the duration for which the surface is exposed to liquid prior to an experiment. A number of investigators have noted that thermal response of the substrate affects nucleation of cavities at neighboring sites. Judd (16) has found that for active cavities spaced between $1∕2$ and 1 bubble departure diameter, the formation of a bubble at the initiating site promotes the formation of bubbles at the adjacent sites (site seeding). However, for separation distance between 1 and 3 bubble departure diameters, the formation of a bubble at the initiating site inhibits the formation of bubbles at an adjacent site (deactivation of sites).

Recently, Dinh et al. (17) have claimed that the presence of cavities on the surface is not required for heterogeneous nucleation and that it is the solid surface energy that determines the propensity of nucleation.

Despite substantial efforts over a period spanning almost $70years$⁠, we do not yet have a generalized correlation or comprehensive model for bubble departure diameter. The key impediments to these efforts have been the lack of knowledge of temperature and velocity fields that vary both temporally and spatially. The temperature and velocity fields, in turn, affect the growth rate and forces that act on the vapor bubble, respectively. The bubble shape changes continuously, hence, the magnitude of forces acting on it. Surface wettability, contribution of microlayer, and lateral and vertical merger of vapor bubbles influence the bubble departure diameter.

Predictions from these correlations are valid only for the limited range over which supporting data have been obtained. Difficulties arise in developing a model for bubble release frequency because it is tied to bubble diameter at departure and growth rate, which, in turn, are dependent on mechanisms responsible for heat transfer into and out of the bubble. Waiting and growth periods are influenced by the temperature field in the solid and/or liquid and, as such, require knowledge of the prevailing phenomena. Cavity-cavity interaction, microlayer evaporation contribution, and bubble merger also influence bubble release frequency.

At this juncture, it appears that mechanistic prediction of nucleate boiling is a hopeless task. However, before accepting this conclusion, it is imperative that we attempt an approach that is vastly different from what we have employed in the past. This approach and the results obtained from this approach are discussed next.

In order to have a credible predictive model of nucleate boiling, one must address four subprocesses as indicated in Fig. 1 and their interactions, which tend to be nonlinear. These subprocesses are density of active nucleation sites, bubble dynamics (which includes bubble growth, merger, and departure), and several mechanisms of heat transfer, such as transient conduction into liquid replacing the space originally occupied by a departing bubble, evaporation at bubble base, and bubble boundary. Convection resulting from surface-tension gradient along the interface and that induced by density difference must be included. The convective motion can be altered by agitation created by vapor bubbles. In most experiments, power to the test surface is controlled. Because of spatial and temporal variation in heat transfer on the fluid side, wall temperature will oscillate. To determine the thermal response of the solid, one must solve a conjugate problem for heat conduction in the substrate. At present, we assume the surface temperature to remain constant. As such, the solid is thermally decoupled. In the same vein, a microfabricated surface is employed on which there is control on the cavities that become active. Thus, our focus will be mainly on bubble dynamics and associated heat transfer mechanisms. For this purpose, we rely on complete numerical simulation of the process—a tool that has been developed only recently.

In simulating the process, the region of interest is divided into micro- and macro regions (Fig. 2) as initially proposed by Son et al. (24). The microregion is the ultrathin liquid film that forms between the solid surface and the evolving liquid-vapor interface. On the inner edge, the microlayer has a thickness of the order of a nanometer, i.e., few molecules of liquid are adsorbed on the surface and do not evaporate. As such, the region further radially inward is considered to be dry. On the outer edge, the microlayer has a thickness of the order of several microns. Heat is conducted across the film and is utilized for evaporation. Lubrication theory is used to solve for the radial variation of microlayer thickness. The macroregion is the region occupied by vapor and liquid, except the microlayer. For the macroregion, complete conservation equations of mass, momentum, and energy are solved for both phases. Interface position is captured through a level-set function.

Governing equations for micro- and macroregions are given in Son et al. (24) and are not repeated here. In analyzing the microregion, continuum assumption is considered to hold good until the film is a few molecules thick. Long-range forces are evaluated through the Hamaker constant, which is also related to the static contact angle. Ramanujapu and Dhir (25) have shown from experiments that in pool boiling, advancing and receding contact angles are within only $±5deg$ of the static contact angle. As such, the use of the static contact angle throughout the bubble evolution process is justified. Capillary pressure gradient is related to change in the curvature of the interface. Recoil pressure resulting from the momentum of vapor leaving the interface being higher than the liquid approaching the interface is included. Inertia terms are neglected in the momentum equation, and convection terms are ignored in the energy equation.

Quasi-static analysis is carried out, and a two-dimensional model for the microlayer is used even in three-dimensional (3D) situations under the assumption that no crossflow occurs in the azimuthal direction.

For the macroregion, it is assumed that fluids are incompressible, flows are laminar, and properties are evaluated at the mean temperature. Vapor is assumed to remain at saturation temperature corresponding to pressure in the bubble. As such, the energy equation is not solved in the vapor space. Finite difference scheme is used in which diffusion terms are solved implicitly, whereas explicit scheme is used for convection terms. Projection method is used to solve for pressure that conserves mass. To prevent numerical oscillations, second-order ENO scheme is adopted for advection terms when solving for the level-set function. In order to increase the rate of convergence, the multigrid method is used. The computational domain is assumed to have a width that is twice the characteristic length, $lo=[σ∕g(ρl−ρv)]1∕2$ and symmetry conditions are imposed on the bounding walls.

Numerical simulation results have been validated with experiments conducted on polished silicon wafers. On this wafer, cylindrical cavities were strategically microfabricated. The wafer is heated on the backside with several strain-gage-type heaters. By controlling power to these heaters, any of the cavities could be nucleated. Details of the experiments are given in Qiu et al. (26).

Son et al. (24) showed a good agreement between predictions from complete numerical simulations and data from experiments with respect to time-dependent vapor bubble shape, growth rate, bubble departure diameter, and growth period. The effect of increased wall superheat, consistent with the data, was to increase the growth rate and reduce the growth period, but increase the bubble diameter at departure. A higher growth rate is due to increased heat transfer at the vapor-liquid interface, whereas the diameter at departure increases because of the increased contribution of liquid inertia. With the increase in liquid subcooling (27), the bubble growth rate decreases and growth period increases, while the bubble diameter at departure decreases.

Fluid properties as well as surface wettability or contact angle are important variables that influence bubble dynamics. Bubble growth histories for water and PF-5060 with widely differing liquid properties as predicted from numerical simulations (28) are plotted in Fig. 3. In the experiments, contact angle for water on silicon was $35deg$⁠, whereas that for PF-5060 was $∼10deg$⁠. Bubble diameter at departure and the growth period for PF-5060 are much smaller than those for water. Numerical predictions do quite well in predicting the data for the two fluids. In order to explicitly demonstrate the effect of contact angle, results of numerical simulations carried out by parametrically varying the contact angle for the two fluids are shown in Fig. 4. It was found that results for the two fluids nearly overlap when bubble diameter at departure and growth period are normalized with their corresponding values for a contact angle of $90deg$⁠. Bubble diameter at departure is seen to decrease linearly as the contact angle is reduced up to about $20deg$⁠. This behavior is similar to that predicted from Fritz’s (10) correlation. However, for contact angles less than $20deg$ when the surface becomes well wetting, the dependence of bubble departure diameter on contact angle is nonlinear. The growth period shows nonlinear behavior for all contact angles. Good agreement of predictions with data obtained using water while varying systematically the surface wettability through oxidation of silicon and with PF-5060 data is seen.

During the growth and departure cycle of a bubble, the wall heat transfer rate as well as contributions of various mechanisms vary both spatially and temporally. Wall heat flux and contributions of individual mechanisms integrated over the heater area supporting the computational domain are shown in Fig. 5. During the waiting period, which was taken to be $16ms$⁠, heat transfer by transient conduction dominates. Its contribution decreases during the early growth period, but increases again as the bubble base shrinks during the detachment phase of the bubble. Natural convection contribution is highest during the early growth period of the bubble when liquid is pushed out radically. Microlayer contributes only during the period the bubble is attached to the heater surface. Heat transfer rate from the surface peaks just before the vapor bubble lifts off from the heater. Time-integrated values suggest that about 50% of the energy from the wall is transferred by transient conduction, 35% by natural convection, and 15% by microlayer evaporation.

Partitioning of wall energy into vapor and liquid, as determined from numerical simulations, is shown in Fig. 6 for one growth and departure cycle of a bubble. The highest rate at which energy is utilized for vapor production occurs when the vapor bubble base diameter is nearly maximum. Conversely, this corresponds to the lowest energy transfer rate into superheating of the liquid. Time integrated values suggest that about 30% of energy from the wall is utilized in vapor production, whereas 70% goes into superheating of liquid.

Vapor bubble merger in the vertical direction at a single nucleation site occurs when the growth rate of a bubble formed at the nucleation site exceeds the rate at which the lower interface of the preceding bubble moves away from the heater surface (29). After merger, the combined vapor mass may detach from the heater surface before the process repeats itself. Such a merger is referred to here as a two-bubble merger process. If, on the other hand, after merger the vapor mass merges with a second succeeding bubble before departure, we call it a three-bubble merger process. The shift from a two- to a three-bubble merger, and so on, depends on the wall superheat and on when vapor leaves the heater in almost a continuous column. Figure 7 shows the results of visual observations and those from numerical simulations for one cycle of the merger of three consecutive bubbles in the vertical direction. The upper set of frames are from visual observations, whereas the lower set of frames are results from numerical simulations.

The individual frames in each figure are from left to right and from top to bottom. After the merger of the departed bubble with the succeeding bubble, the larger vapor mass causes the vapor bubble at the nucleation site to prematurely depart. Thereafter, the second succeeding bubble merges with the vapor mass hovering over the surface. The combined vapor mass goes through several shape changes and departs as a cylindrical bubble. The departing bubble creates a wall jet that impinges on the lower interface of the bubble and forms a dimple. Thereafter, the vapor mass tries to acquire a spherical shape as it moves away from the wall. The rapid acceleration of the vapor mass breaks down the merger process before the cycle repeats itself. The bubble shapes as well as the merger behavior predicted from the numerical simulations is in startling agreement with the visual observations.

As the wall superheat is increased, the number density of the active nucleation sites also increases. With increased nucleation site density, vapor bubbles start to merge laterally as well. Mukherjee and Dhir (30) have carried out numerical simulations and experiments for merger of two bubbles nucleating at adjacent sites. In Fig. 8, the upper set of frames show the results of visual observations, whereas the results of numerical simulations are depicted in the lower set of frames. After the merger of two neighboring bubbles, a mushroom-type bubble with two stems attached to the solid surface forms. A liquid bridge exists between the two stems. As the merged vapor mass tries to acquire a spherical shape as a result of surface tension, vapor tails are formed. The vapor bubble oscillates in size in the plane of the photographs and normal to it before detaching. Numerical simulations capture the essential physics of the process—formation of vapor tails and oscillations in the size of the bubble prior to departure. However, the extent of the trapped liquid in numerical simulations is less than that in experiments.

A quantitative comparison of the predictions from numerical simulations of the growth history of the two bubble merger case with data from experiments is made in Fig. 9. Bubble equivalent diameter in Fig. 9 is the diameter of a perfect sphere having the same volume as the volume of two single bubbles or the merged vapor mass. Numerical simulation results, shown by the solid line, appear to describe well not only the growth history, but also the bubble departure diameter and the growth period.

Abarajith et al. (31) have systematically investigated through numerical simulations, the merger of three inline bubbles and that of three and five bubbles in a plane under low- and microgravity conditions. They found that generally, the predicted vapor bubble departure diameter after merger is smaller than that of a single bubble. In Fig. 10, the numerically predicted growth histories of merger of three inline bubbles and that of a single bubble are compared for earth normal gravity. Equivalent bubble departure diameter and the corresponding growth period for the three-bubble merger case is smaller than that for the single bubble. In order to determine the cause of the premature departure of vapor mass after merger of three bubbles, net force acting on the vapor mass was calculated. This force is taken to be positive when acting upwards and negative when acting downward. Figure 11 shows the variation of the net force on vapor mass during growth and merger of three bubbles and during the growth of a single bubble. For the three-bubble merger case, force changes sign from negative to positive at about $7ms$⁠. This is about the time when the bubble base (Fig. 10) starts to shrink after reaching its maximum value and the vapor bubble enters the detachment phase. However, the single bubble continues to experience negative force and steadily grows for a substantial period beyond when the merged vapor mass starts to detach. The difference between the force acting on the merged vapor mass and the single bubble when the vapor mass, after merger, starts to detach is termed as “lift force.” This force is responsible for the premature departure of the vapor mass after merger. The importance of this force increases as the level of gravity decreases.

Numerical simulations have been carried out (32,33) to predict nucleate boiling heat flux as a function of wall superheat on a surface simulating a commercial surface. On this polished silicon wafer surface, $4cm×4cm$ in area, cylindrical cavities of 10, 7, and $4μm$ dia were fabricated as shown in Fig. 15. Different size cavities were chosen, so that smaller cavities will become active with an increase in wall superheat in a manner similar to that for a commercial surface. In order to accelerate the computations on a large domain, subdomains, as depicted in Fig. 15, were defined around clusters of cavities. While carrying out the computations, no interactions were allowed between neighboring domains. Heat flux in the regions falling between the domains was obtained by interpolating across the values that exist at the boundaries of the domains. The top portion of Fig. 16 shows the visual observation of the boiling phenomenon on the above-described surface when 6–7 cavities were active. The lower figure shows the results of computations at one instant of time. Observed and predicted vapor removal configurations are in a qualitative agreement. Figure 17 shows a comparison of wall heat flux as a function of wall superheat predicted from numerical simulations with data from experiments (26). The solid line is the prediction from numerical simulations when the number of cavities that are active, and their location is given as an input to the simulations. The reported results are post several bubble growth and departure cycles. Predictions are seen to compare reasonably well with the data. The lower curve in Fig. 17 represents the wall energy that goes into production of vapor. Although the fraction of wall energy utilized in production of vapor varies with wall superheat, almost 40% of the energy is utilized for vapor production at a wall superheat of $12°C$⁠. The remaining energy goes into superheating of liquid.

A comparison of the observed and numerically cumputed bubble release pattern when the level of gravity is reduced by almost a factor of 100 is shown in Fig. 18. The same surface was used in the low-gravity experiments as in the experiments conducted at earth normal gravity. The photograph shown in the top portion of Fig. 18 was obtained from the experiments performed in the KC-135 aircraft. Vapor removal configurations predicted from numerical simulations appear to be in general agreement with those found in the experiments. In comparison to earth normal gravity, the departing vapor bubbles are now much larger in size. Predictions and low-gravity heat flux data are compared in Fig. 19. The predictions generally compare well with data, but large scatter in data is seen. Large scatter in data is due to the uncertainty in calculation of heat loss in the KC-135 environment (26). It is found that in comparison to earth normal gravity, the highest heat flux at the highest wall superheat investigated is almost 25% lower. The lower curve represents partitioning of wall energy into vapor. At $12°C$ wall superheat, $∼80%$ of energy is predicted to be consumed in vapor production, whereas the remaining 20% goes into superheating of liquid. It should be stressed that, thus far, the numerically analyzed heat fluxes are low. Further efforts in extending these results to high heat fluxes are needed.

The approach developed here can be extended to commercial surfaces provided the location and number density of active nucleation sites is known. The determination of cavities that become active depends on the shape and size of the defects or cavities present on the surface, which, in turn, requires a detailed knowledge of the topography of the surface. The latter information combined with a criterion such as that developed by Wang and Dhir (14) for entrapment of gas can be used as a first-order approximation to determine the number of active cavities at a particular wall superheat. Alternatively, we could use the correlation of Basu et al. (12) to determine the number density of active sites as a function of wall superheat and contact angle, and randomly distribute these cavities before carrying out numerical simulations of the process.

• 1

  Various correlations and mechanistic models for nucleation site density, bubble departure diameter, and frequency have not always included the effect of all of the variables that influence these subprocesses.

• 2

  Our ability to predict nucleate boiling heat fluxes mechanistically and theoretically is negatively impacted by the fact that in the past the interacting subprocesses have been treated as disjoint processes.

• 3

  Complete numerical simulation of bubble dynamics and associated heat transfer processes appear to represent a potent approach.

• 4

  As further advances are made in computations of all of the physical processes associated with rapidly evolving vapor-liquid interfaces, it should be possible in the near future to predict mechanistically nucleate boiling even at high heat fluxes on microfabricated surfaces simulating a real surface.

• 5

  In the not-too-distant future, there should be an opportunity to carryout a standard problem in which prediction of dependence of nucleate boiling heat flux on wide range of wall superheats is made prior to experiments.

Efforts of all of my graduate students and postdoctoral fellows are gratefully acknowledged. The financial support was provided by NASA under the Microgravity Fluid Physics Program.

$C1$

empirical constant

$Cs$

empirical constant

$cp$

specific heat

$D$

diameter

$Ds$

maximum diameter of largest cavity

$f$

bubble release frequency

$g$

gravitational acceleration

$ge$

earth normal gravity

$h¯$

average heat transfer coefficient

$hfg$

latent heat of vaporization

Ja

Jacob number, $ρlcplΔT∕ρvhfg$

$Ja*$

modified Jacob number, $ρlcplTsat∕ρvhfg$

$k$

thermal conductivity

$K$

ratio of area of influence of bubble to cross-sectional area of bubble

$lo$

characteristic length scale

$M$

molecular weight

$m$

empirical constant

$n$

empirical constant

$Na$

active nucleation site density

$Nmax$

maximum value of $Na$

Nu

Nusselt number

$P$

pressure

$Pc$

critical pressure

Pr

Prandtl number

$q$

heat flux

Re

Reynolds number

$Rp$

surface roughness

$T$

temperature

$t$

time

$td$

growth period

$to$

characteristic time scale

$tw$

waiting period

$ΔT$

temperature difference

$V$

volume

Greek

$α$

thermal diffusivity

$ϕ$

contact angle

$μ$

dynamic viscosity

$ρ$

density

$σ$

surface tension

Subscript

$b$

bubble

$c$

cavity

$d$

departure

ev

evaporation

$l$

liquid

nc

natural convection

ONB

onset of nucleate boiling

sat

saturation

sub

subcooling

$v$

vapor

$w$

wall