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

[+] Author and Article Information
Vijay K. Dhir

Henry Samueli School of Engineering and Applied Science, University of California, Los Angeles, 420 Westwood Plaza, Los Angeles, CA 90095-1597vdhir@seas.ucla.edu

J. Heat Transfer 128(1), 1-12 (Oct 10, 2005) (12 pages) doi:10.1115/1.2136366 History: Received August 31, 2005; Revised October 10, 2005

Over the last half of the twentieth century, a number of purely empirical and mechanism-based correlations have been developed for pool nucleate boiling. Empirical correlations differ from each other substantially with respect to the functional dependence of heat flux on fluid and surface properties, including gravity. The mechanism-based correlations require knowledge of the number density of active sites, bubble diameter at departure, and bubble-release frequency. However, because of the complex nature of the subprocesses involved, it has not been possible to develop comprehensive models or correlations for these parameters. This, in turn, has led to the pessimistic view that mechanistic prediction of nucleate boiling is a hopeless task. However, there is an alternative to the past approaches—complete numerical simulation of the boiling process. Value of this approach for bubble dynamics and associated heat transfer is shown through excellent agreement of predictions with data obtained on microfabricated surfaces on which active nucleation sites can be controlled.

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.

Empirical Correlations

The earliest correlation for pool nucleate boiling is that due to Rohsenow (1). In developing the correlation, Rohsenow reduced the phase-change problem to a single-phase forced convection problem. He proposed that heat from the wall was removed by liquid motion induced by vapor bubbles rising in the liquid pool. Thus, he related the Nusselt number to Reynolds and Pandtl numbers, with the characteristic length being the diameter of the departing bubble. As such, he correlated the data asDisplay Formula

orDisplay Formula
The exponents m and n were found empirically. It was also found that a value of 13 for m would fit the data well. The exponent n was given a value of 1 for water and 1.7 for all other liquids. The proportionality constant Cs depended on heater material and fluid combination, but no attempt was made to relate the value of Cs to surface conditions. Liaw and Dhir (2) have shown that Cs varies with contact angle and its value increases as the contact angle is decreased. Equation 2 can be rewritten asDisplay Formula

Almost three decades later, Stephan and Abdelsalam (3) developed a comprehensive correlation for saturated pool nucleate boiling of different liquids. In developing the correlation, they divided the data into four liquid groups: water, hydrocarbons, cryogenic liquids, and refrigerants. The dimensionless heat flux was related to several dimensionless parameters that depended on fluid and solid properties. The important fluid property groups were identified through regression analysis, and the values of exponents of the property group were obtained by matching predictions with data. In developing the correlation, data from different heater geometrics (such as flat plates, horizontal cylinders, vertical cylinders, etc.) were used. In addition, surface roughness was assumed to be 1μm for all heaters. They also developed a generalized correlation applicable to all liquids,Display Formula

where Dd is the bubble diameter at departure.

Cooper (4) has proposed a much simpler correlation for saturated pool nucleate boiling. His correlation employs reduced pressure, molecular weight, and surface roughness as the correlating parameters. The correlation for a flat plate isDisplay Formula

In Eq. 5, Rp is the surface roughness measured in microns. Cooper suggested that for application of the correlation to the horizontal cylinders, the lead constant on the right-hand side should be increased to 95. Also, q is given in watts per square meters.

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.

Mechanism-Based Correlations

In partial nucleate boiling or in the isolated bubble regime, transient conduction into liquid adjacent to the heater surface is one of the important mechanisms for heat transfer from an upward-facing horizontal surface. After bubble inception, the superheated liquid layer is pushed outward and mixes with the bulk liquid. During the growth phase, the bubble acts like a pump in removing the superheated liquid from the heater surface and replacing it with cold liquid upon departure. This mechanism was originally proposed by Forster and Greif (5). Combining the contribution of transient conduction on and around nucleation sites, microlayer evaporation underneath the bubbles and natural convection on inactive areas of the heater, an expression for partial nucleate boiling heat flux is obtained asDisplay Formula

Only the first two terms in Eq. 6 were included in the original model proposed by Mikic and Rohsenow (6). The addition of the last term on the right-hand side of Eq. 6 was suggested by Judd and Hwang (7). This term accounts for evaporative heat transfer from the microlayer underneath the bubble. For Eq. 6 to serve as a predictive tool, three key parameters, aside from constant K, accounting for the ratio of the area of influence of a bubble to the cross-sectional area of the bubble that must be known, are: number density of active sites Na, bubble diameter at departure Dd, and bubble release frequency f. During the last four decades, a number of attempts have been made to develop correlations or models for these parameters, but with limited success.

Number Density of Nucleation Sites

Mikic and Rohsenow (6) have proposed that on commercial surfaces, the cumulative number of active sites per unit heater area can be assumed to vary in partial nucleate boiling as,Display Formula

where Ds is the diameter of the largest cavities present on the surface and m1 is an empirical constant. The nucleating cavity diameter Dc is related to wall superheat asDisplay Formula
Bier et al. (8), on the other hand, have deduced an expression for active site density from their heat transfer data asDisplay Formula
where Nmax is the maximum value of Na, which occurs at Dc=0. The value of m2 was found to depend on the manner in which the surface was prepared. Kocamustafaogullari and Ishii (9) have correlated the cumulative nucleation site density reported by various investigators for water boiling on a variety of surfaces at pressures varying from 1to198bars asDisplay Formula
whereDisplay Formula
andDisplay Formula
In Eq. 11, the bubble diameter at departure Dd is obtained by using Fritz’s correlation (10), and at low to moderate pressures, Dc is given by Eq. 8.

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

According to their correlation, the density of active nucleation sites depends only on wall superheat and contact angle.

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

Bubble Diameter at Departure

A large number of correlations and mechanistic models have been proposed in the literature for diameter to which a bubble grows before departure. Fritz (10) correlated the bubble departure diameter Dd by balancing buoyancy, which acts to lift the bubble away from the surface, with surface tension force, which tends to hold the bubble to the wall. The resulting correlation is given asDisplay Formula

where ϕ is the contact angle measured in degrees. Significant deviations with respect to prediction from the above equation have been reported in the literature, especially at high pressures. Several other expressions, either obtained analytically or by considering various forces acting on a bubble, have been reported for bubble departure diameter. These expressions, however, are not always consistent with one another with respect to dependence of bubble departure diameter on independent variables. Some investigators report an increase in bubble diameter at departure with wall superheat, whereas others find the bubble diameter at departure to be insensitive to or decrease with wall superheat. In addition, contrary to the commonly held view, the work of Cooper et al. (18) suggests that surface tension may assist bubble departure. Cole and Rohsenow (19) correlated bubble diameter at departure at low pressures asDisplay Formula
andDisplay Formula
where Ja*=ρlcplTsatρvhfg.

Kocamustafaogullari (20) has developed a correlation for bubble departure diameter that includes high-pressure dataDisplay Formula

Gorenflo et al. (21) have proposed an expression for bubble departure diameter at high heat fluxes asDisplay Formula
Different values of constant C1 were used for different liquids. It should be noted that correlations of Fritz (10), Cole and Rohsenow (19), and Kocamustafaogullari (9) suggest that bubble diameter at departure varies at g12, where as according to Eq. 18, it varies as g13. In addition, the correlation of Gorenflo et al. (21) suggests thermal diffusivity of liquid as well as Jakob number to be important parameters. According to the correlation of Kocamustafaogullari (9), the liquid to vapor density ratio is another important variable.

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.

Bubble Release Frequency

Bubble release frequency is the inverse of the sum of waiting and growth periods. Both the waiting time tw and growth period td are difficult to determine quantitatively. Waiting time depends on the temperature field in the solid, as well as in the liquid in the vicinity of the nucleation site. Similarly, growth period depends on the rate of evaporation at the bubble base and around the bubble and the bubble diameter at departure. Often data for bubble release frequency are correlated under the assumption of twtd or vice versa. Correlations have also been developed for vapor superficial velocity or the product of bubble diameter at departure and frequency. Zuber’s (22) expression for vapor superficial velocity is often used. According to it,Display Formula

Malenkov (23) has proposed thatDisplay Formula
whereDisplay Formula

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

Single Bubble Dynamics

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.

Bubble Merger and Formation of Vapor Columns

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.

Bubble Merger and Formation of Mushroom Type Bubbles

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

The computed shapes of the vapor mass during merger of three bubbles placed at the corners of an equilateral triangle are shown in Fig. 1. Departure diameter and the vapor bubble after merger and corresponding growth period depend on the spacing and orientation of nucleating cavities. Figure 1 shows predictions from numerical simulations of bubble diameter at departure as a function of spacing. Both the bubble diameter at departure and the spacing have been normalized with the characteristic length defined asDisplay Formula

Bubble departure diameter decreases as the spacing between cavities increases until the spacing is equal to Dds4, where Dds is the departure diameter of a single bubble. For a contact angle of 54deg, with water as the test liquid at one atmosphere pressure, Dds is given asDisplay Formula
Thereafter, the departure diameter increases until the spacing is large enough so that bubbles do not merge. The limiting values of departure diameter being equal to single bubble diameter, when spacing between cavities is zero, and is equal to 33Dds, when the spacing is 1.4lo.

The corresponding variation of the growth period is shown in Fig. 1. The growth period is nondimensionalized with the characteristic time defined asDisplay Formula

Growth period is found to show a behavior similar to the bubble diameter at departure. It is noted that depending on the magnitude of the spacing, bubble departure diameter can vary by a factor of two, whereas the growth period by a factor of three.

Nucleate Boiling Heat Flux

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. 1. 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. 1, 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. 1 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 1 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. 1 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. 1. 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. 1 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. 1. 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.

  • 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.
  • 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.
  • Complete numerical simulation of bubble dynamics and associated heat transfer processes appear to represent a potent approach.
  • 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.
  • 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.


empirical constant


empirical constant


specific heat




maximum diameter of largest cavity


bubble release frequency


gravitational acceleration


earth normal gravity


average heat transfer coefficient


latent heat of vaporization


Jacob number, ρlcplΔTρvhfg


modified Jacob number, ρlcplTsatρvhfg


thermal conductivity


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


characteristic length scale


molecular weight


empirical constant


empirical constant


active nucleation site density


maximum value of Na


Nusselt number




critical pressure


Prandtl number


heat flux


Reynolds number


surface roughness






growth period


characteristic time scale


waiting period


temperature difference




thermal diffusivity


contact angle


dynamic viscosity




surface tension












natural convection


onset of nucleate boiling









Copyright © 2006 by American Society of Mechanical Engineers
View article in PDF format.



Grahic Jump Location
Figure 1

Mechanistic prediction of nucleate boiling heat transfer

Grahic Jump Location
Figure 2

Macro and micro regions of the mathematical model for numerical simulation

Grahic Jump Location
Figure 3

Comparison of bubble departure diameter and bubble growth time for water and PF5060

Grahic Jump Location
Figure 4

Variation of the normalized departure diameter and departure time with contact angle (Fluids: water and PF5060, ΔTw=8°C, ΔTsub=0°C, g=1.0ge).

Grahic Jump Location
Figure 5

Partitioning of wall heat transfer into natural convection, transient conduction, and microlayer heat transfer rates

Grahic Jump Location
Figure 6

Variation of the heat transfer rates from wall, into liquid, and into vapor with time

Grahic Jump Location
Figure 7

Comparison of numerical and experimental bubble shapes during vertical merger (29) (Fluid: water, ΔTw=10°C, ΔTsub=0.0°C, g=1.0ge, ϕ=38deg).

Grahic Jump Location
Figure 8

Comparison of numerically predicted bubble growth with experimental data of Mukherjee and Dhir (30) (Fluid: water, ΔTw=5°C, ΔTsub=0°C, g=1.0ge, spacing=1.5mm)

Grahic Jump Location
Figure 9

Comparison of numerically predicted bubble growth with experimental data of Mukherjee and Dhir (30)

Grahic Jump Location
Figure 10

The variation of equivalent bubble diameter with time for single and three inline bubbles (Fluid: water, ΔTw=10°C, ΔTsub=0°C, g=1.0ge, ϕ=54deg, spacing=1.25mm)

Grahic Jump Location
Figure 11

The variation of force acting in the vertical direction with time for single and three inline bubbles (Fluid: saturated water, ΔTw=10°C, g=1.0ge, ϕ=54deg, spacing=1.25mm)

Grahic Jump Location
Figure 12

Growth, merger, and departure of three bubbles located at the corners of an equilateral triangle (Fluid: saturated water, spacing=1.25mm, ΔTw=10°C, g=1.0ge, ϕ=54deg)

Grahic Jump Location
Figure 13

The variation of nondimensional bubble departure diameter with cavity spacing for bubbles placed on the corners of equilateral triangle (Fluid: saturated water, ΔTw=10°C, g=1.0ge)

Grahic Jump Location
Figure 14

The variation of nondimensional time period of growth with cavity spacing for bubbles placed on the corners of equilateral triangle (Fluid: saturated water, ΔTw=10°C, g=1.0ge)

Grahic Jump Location
Figure 15

Simulated commercial surface

Grahic Jump Location
Figure 16

Comparison of experimental and predicted bubble shapes during nucleate boiling on a simulated commercial surface for (Fluid: saturated water, ΔTw=7°C, g=1.0ge)

Grahic Jump Location
Figure 17

The variation of wall heat flux with wall superheat for water at g=1.0ge

Grahic Jump Location
Figure 18

Comparison of experimental and predicted bubble shapes during nucleate boiling on a simulated commercial surface (Fluid: saturated water, ΔTw=7°C, g=0.01ge)

Grahic Jump Location
Figure 19

The variation of wall heat flux with wall superheat for water at g=0.01ge




Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In