- Email: [email protected]
Bubble nucleation, growth, and departure: A new, dynamic understanding
International Journal of Heat and Mass Transfer 145 (2019) 118803
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Bubble nucleation, growth, and departure: A new, dynamic understanding H. Jeremy Cho a,b, Evelyn N. Wang a,⇑ a b
Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA, United States University of Nevada, Las Vegas, 4505 S Maryland Pkwy, Las Vegas, NV, United States
a r t i c l e
i n f o
Article history: Received 11 February 2019 Received in revised form 23 September 2019 Accepted 27 September 2019
Keywords: Boiling Bubble Nucleation Growth Departure Wetting
a b s t r a c t The bubble nucleation, growth, and departure cycle is a fundamental aspect of nucleate pool boiling. While much research on this subject has been performed in the previous century, new correlations and models have not been developed in light of important results in simulations and experiments in recent decades. In this work, we provide an updated understanding of nucleation, growth, and departure with analytical models that are validated by recent work. We found that nucleation at incipience is correlated to departure size, suggesting that the momentum-induced bubble snap-off at departure may be responsible for the size of trapped vapors in cavities. We also developed a bubble growth model that takes into account a thermal boundary layer whereas previous works typically considered a uniform superheating. In addition, we developed a bubble departure description, wherein the velocity boundary layer is responsible for closing the base of the bubble. Finally, we provide a methodology to calculate bubble growth and departure for contact angle changing with time, which is relevant in studying surfactantenhanced boiling. Ó 2019 Published by Elsevier Ltd.
1. Introduction Boiling is utilized in a variety of industrial and domestic applications. Despite more than half a century of research, many uncertainties about the process still remain. Boiling is a highly complex process involving such phenomena as turbulent flows, large temperature gradients, liquid–vapor phase change, and contact line dynamics. To gain a simplified understanding, much of the research on boiling fundamentals has focused on the bubble cycle: the subsequent events of bubble nucleation on a boiling surface, bubble growth, and bubble departure from the surface. A single bubble cycle transfers an amount of heat per time qsingle ; knowing how to calculate qsingle allows one to calculate heat transfer for the entire boiling surface. The heat transfer during boiling is nqsingle where n is the nucleation site density. While many theories, models, and correlations were developed from the 1950s to 1970s that allow for the calculation of qsingle [1–3], the advent of computer simulations and more sophisticated experiments in recent decades, have revealed a richer physical picture of the bubble cycle and accurate bubble cycle data [4–6]. In particular, it ⇑ Corresponding author. E-mail address: [email protected] (E.N. Wang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118803 0017-9310/Ó 2019 Published by Elsevier Ltd.
has been shown that heat transfer primarily occurs through the liquid–vapor interface either through the microlayer or from the bulk surrounding liquid [7]. In light of these new experiments and simulations, we present a way to determine qsingle that is based on new physical insights and validated by the latest research. We do not discuss how to calculate nucleation site density since it is highly surface-dependent [8] and is often a feature that is designed for in modern boiling surfaces [9].
1.1. Background on nucleation The first stage of the bubble cycle is bubble nucleation. According to trapped vapor theory (also known as entrapped vapor theory), cavities on the surface trap air bubbles; these trapping sites serve as nucleation sites in boiling [10]. Because of the convex curvature of trapped vapor bubbles at cavities, the vapor pressure inside the bubble is higher than the surrounding liquid. As such, if the surrounding liquid is at the saturation temperature, its equilibrium vapor pressure would not match the vapor pressure inside the trapped bubble. However, if the liquid were sufficiently superheated, its equilibrium vapor pressure would be greater than the bubble vapor pressure, causing the bubble to grow. The particular condition when this bubble growth—or nucleation—occurs is when
2
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
Nomenclature
al
liquid thermal diffusivity molar gas constant velocity boundary layer clv liquid-vapor surface tension l chemical potential ql liquid density qv vapor density qc T Ja modified Jakob number ( lqp;lh sat ) v fg ReL Reynolds number h contact angle (radians) A area abase acceleration of base contact line Aq surface area where heat transfer occurs on growing bubble C trans transition time scaling factor Cb growth rate scaling factor cp;l liquid specific heat base diameter Dbase Dequiv;depart bubble equivalent diameter at departure Dequiv bubble equivalent diameter g acceleration due to gravity hboil boiling heat transfer coefficient R dv
the chemical potential of the liquid matches the chemical potential of the vapor. This is often approximated by combining the Clausius–Clapeyron relation with the Young–Laplace equation [11]:
DT
2clv T sat v fg : hfg Rb
ð1Þ
Thus, if the radius of curvature, Rb , of a nucleation site is known, then the superheat at which bubble nucleation occurs can be calculated using Eq. (1), which we denote as the Clausius–Clapeyron-Y oung–Laplace (CCYL) equation. Two common ways to obtain the bubble radius of curvature, Rb , is through the Lorenz and Hsu models [12,13]. However, we have found that these models do not agree with experiments of nucleation activation at various pressures. Thus, we provide a simple correlation for Rb based on bubble departure size that more accurately reflects these experiments. 1.2. Background on growth The next stage of the bubble cycle is bubble growth. The growth rate depends on wetting behavior and superheat. Plesset and Zwick [14] modeled this growth assuming that heat transfer from a superheated liquid to a bubble is diffusion-limited. In this diffusion-limited case, Rb / t 1=2 . Mikic et al. further extended this model to include an inertia-dominated growth regime where Rb / t at very short timescales before transitioning to the diffusion-limited growth regime where Rb / t 1=2 [2]. van Stralen et al. considered growth with a thin evaporation microlayer underneath the bubble and also found that Rb / t 1=2 [15]; thus, it is difficult to differentiate between microlayer and diffusion-limited growth [16]. Nonetheless, many others have investigated the relative heat transfer contribution from the microlayer [7,17,18]. Regardless of the microlayer contribution, heat transfer involves conduction through a liquid layer in both microlayer and bulk DT ffi nature is regions. Thus, heat transfer of a semi-infinite, q00 pklffiffiffiffiffiffi pal t
present around the bubble, leading to a Rb / t 1=2 dependence. Recently, experimentally-validated computer simulations of single-bubble growth that take into account microlayer effects
hfg kl LBo LGr n n P Psat qsingle R2 T t td T sat t trans t wait V
v fg
z
specific enthalpy of liquid-vapor phase change liquid conductivity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi characteristic Bond length scale ( gðqclv q Þ) l v characteristic Grashof length moles of ideal gas nucleation density pressure saturation pressure single bubble heat transfer coefficient of determination temperature time time at which bubble growth transitions to t 1=4 growth saturation temperature transition time when constant base closing force starts bubble wait time volume specific volume of liquid-vapor phase change dimension normal to surface
and macro-scale effects indicate that the Rb / t1=2 dependence breaks down at longer timescales closer to departure [4–6]. In this work, we took into account the thermal boundary layer to develop an analytical model of bubble growth that agrees very well with experimentally-validated computer simulations. We consider semi-infinite conduction through the liquid around the bubble— assumed to take a pendant bubble geometry—as a means to lump the effects of bulk and microlayer heat transfer together since they lead to similar growth signatures. Our model shows an initial diffusion-limited Rb / t1=2 growth characteristic of bulk and microlayer heat transfer, and an eventual Rb / t 1=4 at longer timescales.
1.3. Background on departure The final stage of the bubble cycle is bubble departure. The departure of bubbles from the boiling surface is the result of an interaction of buoyancy, momentum, and surface tension forces. Most correlations are based off of a static force balance between surface tension and buoyancy where the departure diameter is approximately on the order of the characteristic Bond length: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi LBo clv =ððql qv Þg Þ. Fritz developed a simple correlation [19] where
Dequiv;depart ¼ 1:19hLBo :
ð2Þ
Later correlations incorporated wall superheat effects but did not incorporate contact angle effects [1]. Currently, there are no departure models that consider the time-dependent contact angle behavior that are present due to dynamic contact line movement, surfactants [20], or microlayer evaporation [21]. To address this limitation, rather than develop a generic bubble departure size correlation, we provide a method to calculate the time of departure based on the movement of the contact line taking into account both superheat and contact angle—including dynamic contact angles. Our departure description incorporates effects of the velocity boundary layer and agrees well with experimentally validated numerical simulations of single bubbles [4–6].
3
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
2. Bubble nucleation
equations are solved simultaneously to determine Rb ; ql , and qv from a known DT; T sat , and P sat :
2.1. Limitations of existing models
lIAPWS ðT sat þ DT; ql Þ ¼ lIAPWS ðT sat þ DT; qv Þ; PIAPWS ðT sat þ DT; ql Þ ¼ Psat 2 c ð T þ DT Þ PIAPWS ðT sat þ DT; qv Þ ¼ Psat þ lv :
To evaluate the limitations of existing nucleation models, we compared the models’ predictions of radii of curvature, Rb , to experimentally derived values as a function of pressure. Boiling studies at various pressures show that activation superheats increase with decreasing pressure [22–25] (Fig. 1). From boiling curve data, it is possible to determine Rb from the superheat at which nucleation occurs. This can be done using a ClausiusClapeyron-Young-Laplace (CCYL) equation (Eq. (1)). However, since the CCYL approximation is only valid at small superheats [26], we chose to perform an exact calculation of equivalent potentials. This calculation is possible with steam tables. We have performed our analysis with water properties as given by the International Association for the Properties of Water and Steam (IAPWS) [27]. Here, we equate chemical potentials (or Gibbs free energies) between the vapor and liquid states, and relate pressure to bubble curvature, Rb , using the Young–Laplace equation in an approach we denote as the chemical-equilibrium-Young–Laplace (CEYL) method. To apply the CEYL approach, the following
Fig. 1. Calculated bubble radii of curvature as a function of pressure from experiments [22–25,36] show that Rb scales with the bubble departure size (blue curve) instead of advancing-front models—red dashed curve for constant volume (Eq. S3) [12] and red curve for constant Helmholtz energy (Eq. S12). We have taken data from the literature and estimated the temperature at incipient boiling from the boiling curve using a least-squares fitting, then applied Eqs. (3)–(5) to find Rb;incipient . It is evident that Rb;incipient decreases with increasing pressure. This is contrary to the advancing-front description, which predicts either an invariance of Rb;incipient with pressure (constant volume model) or a decrease in Rb;incipient with decreasing pressure (constant Helmholtz model). Thus, there must be some physical mechanism for vapor trapping outside of simple wetting in an advancing liquid front. The advancing-front models plotted uses a conical cavity with / ¼ 25 ; h ¼ 50 , and 1=3 Rcav ¼ 5 lm. The Hsu model [13] plotted has C Hsu ¼ 1:6 and dT ¼ 0:5 gbmaDlT and predicts the opposite trend with pressure. The idea of trapped gases compressing due to increased saturation pressure is an insufficient predictor of incipience as modeled using the ideal gas law with capillary pressure (Eq. (6), purple curve). The critical pressure of water, P crit , is 22.1 MPa. Our model (blue curve) is likely invalid very close to the critical point since the departure correlation (Eq. (9)) diverges due to an infinite specific heat [37]; however, we show agreement with experimental results up to 150 MPa. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Rb
ð3Þ ð4Þ ð5Þ
The incipient superheat, DT, used for plotting literature data in Fig. 1 was estimated by performing a linear regression to the region of the boiling curve data where a steep rise occurs; the intercept with the superheat axis serves as an approximate incipient superheat. Eqs. (3)–(5) can also be used by any other fluid as long as the equation of state (steam table equations) are known. Since solving these equations is computationally inexpensive, we recommend using the CEYL method whenever possible, especially at higher superheats. 2.1.1. Advancing-front models Lorenz developed a geometric model of the trapped vapor geometry assuming a quasi-static advancing liquid front over a conical cavity [12]. This description depends entirely on contact angle and cavity geometry. Tong et al. improved upon Lorenz’ model by recognizing that the dynamic contact angle should be involved in this cavity-wetting behavior [26]. Some relaxation assumption about the bubble geometry must be made as the bubble relaxes from the moment it is trapped (state 1 in Fig. S1) to the embryo state where a trapped bubble of curvature Rb is formed (state 2 in Fig. S1). In Lorenz’ model (Section S1), a constant volume between state 1 and state 2 is assumed (V 1 ¼ V 2 ); thus, the radius of curvature can be related to cavity geometry and contact angle, h (Eq. S3). Assuming that the surrounding liquid is uniformly at a superheated temperature, the CCYL relation (Eq. (1)) can be used to relate the superheat of nucleation activation, DT, to cavity geometry and contact angle (Eq. S4). The result is that nucleation activation superheat decreases with increasing contact angle—a typical observation in experiments [28]. However, Lorenz’ description does not capture any dependence on pressure (Fig. 1, red dashed line)—we did not include pressure-dependent contact angle since this effect is likely minimal [29]. Even using a different relaxation assumption that is more physically reasonable, where V 1 – V 2 and the Helmholtz free energy is conserved (isothermal condition, which is derived in Section S1), results in a pressure dependence that is completely counter to the experimental data (Fig. 1, red solid line). Therefore, there is little support for the quasi-static advancing front description with uniformly superheated liquid from literature data. One possible explanation for the lowering of incipient radii at higher pressures could be the compression of vapor regions that were formed during cavity wetting at lower pressures. To investigate this possible mechanism, we modeled the compression of an ideal gas with elevated bubble pressure from capillarity:
2c ðT sat ; P sat Þ 4 3 Psat þ lv pRb ¼ nRT: Rb 3
ð6Þ
After solving Eq. (6) for Rb using IAPWS properties across different saturation conditions, it is clear that simple compression of trapped vapor regions underpredicts the range of incipient radii from experimental data (Fig. 1, purple curve). Another possible explanation for this decrease in incipient radii with pressure could be the dissolution of gases into water at elevated pressures. Dissolution of gases in water follows Henry’s Law, where solubility is proportional to pressure [30]. At saturation conditions, this proportionality still holds since the Henry’s Law constant does not vary much with temperature [30]. Thus, at elevated pressures, it is likely that bubbles would become completely dissolved since solubilities would be
4
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
increased many times over. Regardless, the dissolution of bubbles should be highly system-specific since the relative amounts of water and entrapped vapor can vary, and hence the amount of gas dissolution can vary. From the experimental data plotted in Fig. 1, however, it is clear that bubble incipience follows a general trend across independent data sources; thus, the effects of solubility, if any, appear to be minimal. 2.1.2. The Hsu bubble viability model Hsu [13] and others [31] considered a growing bubble within a quiescent superheated thermal boundary layer. After every departure event, a one-dimensional, transient conduction occurs within the thermal boundary layer with thickness dT . Hsu showed that if the liquid within this layer does not sufficiently heat the trapped bubble embryo according to the CCYL relation, then the site remains inactive. Mathematically, there is a discriminant (quantity within a square-root) that must be positive for a solution of a bubble actively growing to exist. This condition, which we denote as bubble viability, is
DT ¼
8clv T sat ð1 þ cos ðc þ hÞÞ hfg qv dT
ð7Þ
when at saturation conditions. Combining the Hsu superheat criterion (Eq. (7)) with the CEYL equations (Eqs. (3)–(5)) yields a relationship between DT and Rb as a function of pressure (Fig. 1, green curve). It is clear that the Hsu criterion is not consistent with experimental results as it predicts a decrease in Rb with a decrease in pressure. Thus, there is little experimental support for the description of quiescent, one-dimensional, transient conduction in the thermal boundary layer driving bubble nucleation. 2.2. New correlation for the bubble radius of curvature at nucleation In contrast to the advancing-front and Hsu bubble viability models, we hypothesized that convective motion may play a role in bubble nucleation. More convective motion implies that bubbles depart earlier (smaller bubble departure diameter, Ddepart ) and that the velocity and thermal boundary layer are thinner. With thinner thermal boundary layers, trapped vapor embryos would be exposed to higher superheats, leading to smaller Rb at nucleation. Thus, to test the link between convection and nucleation, we examined whether Rb decreases with decreasing bubble departure size. Using the departure diameter correlation by Cole and Rohsenow [32] where
Ddepart ¼ 1:5 104
clv
ðql qv Þg
1=2
Ja5=4
ð8Þ
and
Ja
ql cp;l T sat ; qv hfg
ð9Þ
is the modified Jakob number to describe the state of saturation, it was found that Rb / Ddepart : 3
Rb 1:6 10 Ddepart :
ð10Þ
The agreement with experimental data over several decades of pressure measurement (Fig. 1, blue line) confirms our hypothesis that increased convection causes smaller bubble embryos to nucleate, perhaps through a snap-off mechanism. Thus, using Eq. (10), a correlation for bubble departure (a new departure correlation is presented in Section 4), and the CCYL equation (Eq. (1)) or CEYL equations (Eqs. (3)–(5)), it is possible to calculate the superheat at which a nucleation site will activate.
3. Bubble growth Given the importance of convection as suggested by our nucleation description, we modeled bubble growth in the presence of velocity and thermal boundary layers. In our simplified modeling approach, we assume that growing bubbles can be approximated by the static pendant bubble shape in an unchanging thermal and boundary environment, the properties of which are determined by the superheat. We note that these assumptions and simplifications are not completely physical; however, we validate this simplified description through agreement with available data on bubble growth from literature sources. To quantify bubble size, we use equivalent radii and diameters, 1=3 Requiv and Dequiv , respectively: Requiv ¼ Dequiv =2 ¼ 43p V . As Mikic et al. have found [2], bubble growth can be described by two regimes. At short timescales, there is inertially dominated growth in which Requiv / t and is described by the Rayleigh equation. At longer timescales, growth is heat diffusion limited in which Requiv / t1=2 , consistent with Plesset and Zwick’s description [14]. Since the timescale at which this transition occurs is less than one millisecond whereas bubble departure times are typically in the tens of milliseconds, we assumed that growth is only diffusion-limited. While Plesset and Zwick and Mikic et al. did not consider contributions from microlayer evaporation, van Stralen et al. showed that microlayer growth also has a t 1=2 dependence [3]. Thus, growth from the microlayer should mathematically resemble the diffusion-limited growth as described by Plesset and Zwick [16]. As a result, in our model, any microlayer effects are effectively taken into account by a prefactor, C b . Plesset and Zwick stated that the conductive heat transfer into and around a growing spherical bubble is equivalent to the latent heat times the mass change rate of the bubble:
kl DT dV : Alv pffiffiffiffiffiffi qv hfg dt al t
ð11Þ
From this relationship, we can derive an expression for the bubble growth rate as a function of fluid properties and contact angle for a bubble with a spherical cap geometry:
Requiv ðt Þ h 1=3 ¼ Rðt Þ 4 2 cos 2 cos4 2h cos ðhÞ
pffiffi 8kl DT pffiffiffiffiffi t : ð3 þ 2 cos ðhÞ cos ð2hÞÞqv hfg al
ð12Þ
Though the spherical cap geometry is simple, it is unphysical as it neglects the effects of gravity. Bubbles adhered to a horizontal wall in the presence of gravity will have a pendant geometry (Section S3). Furthermore, the pendant bubbles would be in the presence of a thermal boundary layer. Thus, we developed a growth model that assumes a pendant bubble in the presence of a thermal boundary layer. As shown in Fig. 2, we have used a simple thermal boundary layer model where the liquid temperature linearly decreases from T ¼ T wall at z ¼ 0 (at the solid-liquid interface) to T ¼ T sat at z ¼ dT so
z DT ¼ DT wall 1 dT
ð13Þ
where dT is the thermal boundary length. We also assume that the temperature of the vapor bubble is at T sat . Use a scaling analogous to Eq. (11), but with the incorporation of the thermal boundary layer,
Z Aq
kl ðT ðzÞ T sat Þ p dDequiv pffiffiffiffiffiffi dA D2equiv qv hfg 2 dt al t
ð14Þ
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
5
more closely matches this flattening due to its transition to a Dequiv / t 1=4 dependence. Furthermore, our growth model can be extended to cases when surface tension is time-dependent (Section S3). 4. Bubble departure In our treatment of bubble departure, we track the movement of the three-phase contact line with time based on our growth model. To do so, we relate the diameter of the three-phase contact circle, Dbase , to the bubble size:
Dbase ¼ 0:8hDequiv
Fig. 2. Heat transfer around a growing bubble within a thermal boundary layer. A bubble with height lb having a pendant shape grows as a result of lateral heat conduction, q00 , at its bottom. The heat flux is larger near the surface and decreases to zero at the thermal boundary layer height, dT . The temperature profile in the thermal boundary layer is linear.
where dA is a differential bubble surface area at a height z and Aq is the surface area where heat transfer occurs (any bubble surface that is within the thermal boundary layer). By incorporating the pendant bubble geometry (Section S3), the growth rate can be solved analytically:
8 pffi 2 > < 2dc T 1 exp C b bc DdT wall t 0 6 t < td T Requiv ðt Þ ¼ pffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi > : 2 C bdT DT pffiffit d2 ðln 41Þ t P t d b wall 2c2
ð15Þ
where
h c ¼ 2 cos ; 2 td ¼
dT ln 2 2
C b bc DT wall
ð16Þ !2 :
ð17Þ
Here, we assume the contact angle, h, to be the static contact angle. Using the dynamic contact angle could lead to more accurate modeling. However, the use of the static contact angle in this case is reasonable since we are considering the growth of isolated, single-bubbles, which tends to occur at lower superheats where dynamic contact angle effects are small [33]. We also provide a methodology for incorporating dynamic contact angles in our Supplementary Information (S4). When t < td , the height of the bubble, lb is less than the thermal boundary thickness, dT , and when t > t d ; lb > dT . Here, Eq. (15) has the behavior of t1=2 growth when the bubble is much smaller than the boundary length (same result as Plesset and Zwick where a uniform superheating is assumed) and t1=4 growth when the bubble is larger than the boundary length. By performing a least-squares fit of Eq. (15) to nine data sets of experimentally verified single-bubble simulations from several studies [4–6] (Fig. 4), we have found that the scaling factor C b ¼ 0:534 is appropriate for saturated water at atmospheric conditions. For comparison, we also fit the bubble growth model for a spherical cap under uniform superheating (Eq. (12)) to the same data sets. We found that the agreement of our pendant bubble in
ð18Þ
where h is radians (applicable for a variety of Bond numbers). This approximation is close to the actual Young-Laplace solution (Fig. S2). As suggested by our description of bubble nucleation, we hypothesized that convective forces determine the point at which bubble departure occurs. Closer observation of the base diameter velocity as a function of time for a bubble simulated by Mukherjee and Kandlikar [5] as shown in Fig. 3 suggests that a momentum force may be responsible for bubble departure. Initially, the base diameter velocity is high and positive as the bubble grows rapidly. However, the rate at which this base radial velocity changes is not constant (inward radial acceleration decreases). At some transition time, ttrans , the radial velocity of the contact line decreases linearly with time (constant inward radial acceleration). This constant inward radial acceleration is evidence of a convective force, consistent with our hypothesis. Thus, we modeled the evolution of Dbase with time such that, initially, Dbase increases according to pendant bubble geometry (Eq. (18)). After the transition time, t trans , the acceleration of the base closing, abase , is maintained at its value at ttrans (second derivative of Dbase with respect to time). This causes Dbase to decrease until the point of departure when Dbase ¼ 0. Until the point of departure, growth proceeds according to Eq. (15). We can determine ttrans from a velocity boundary layer scaling analysis. According to the boundary layer thickness equation, pffiffiffiffiffiffiffi dv L= ReL where dv is the velocity boundary layer thickness and L is the boundary layer length. Using ReL ¼ L2 =ðt ml Þ, the characteristic timescale associated with boundary layer development is d2v =m. Assuming ttrans is on the order of this time scale, ttrans d2v =m. Assuming the velocity boundary layer thickness
boundary layer model (Eq. (15)) was superior (R2 ¼ 0:994) to that of spherical cap bubble in uniform superheat model (Eq. (12)) (R2 ¼ 0:985). The Dequiv / t1=2 growth associated with uniform superheating does not closely follow the flattening of the growth curve at longer times (Fig. 4). On the other hand, our growth model
Fig. 3. Base diameter velocity shows a transition to a constant acceleration zone. Data of bubble growth and departure with DT wall ¼ 10 C and h ¼ 54 from a numerical simulation by Mukherjee and Dhir [5] shows the base diameter velocity transitioning to a constant acceleration zone.
6
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
(dictated by the Grashof number) can be modified by the contact angle (smaller contact angles should have thinner layers), then 1=3 2 . The resulting ttrans is then dv h gbm DT l
t trans ¼ C trans h2
m1=3
ð19Þ
ðgbl DT Þ2=3
where
C trans ¼ 0:261
ð20Þ
was found to be an appropriate scaling factor for saturated water at atmospheric pressure. However, the analysis above is in the limit where convection is important. In the limit of small bubbles (very hydrophilic surfaces; low h), bubble departure may simply result from a static force balance between buoyancy and surface tension:
Dequiv;depart
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 24 sin h clv ¼ : 3 2 þ 3 cos h g ðql qv Þ
ð21Þ
One way to determine the relative importance of convection versus the static force balance is to compare the characteristic length scales of the contact angle-dependent Bond and Grashof numbers. As
hLBo =LGr ¼
1=2 1=3 h lv bl DT 1=3 g 1=6 ð l v Þ1=2 2=3
c q
q
m
ð22Þ
approaches unity, then the effect of the velocity boundary layer is reduced and bubble departure should be dictated by the static force balance between surface tension and buoyancy forces. For the nine different data sets in Fig. 4, we show values of hLBo =LGr . It can be seen that the lowest value of hLBo =LGr is for Nam et al., which happens to be where our departure model is least successful. Thus, our model should be more accurate for situations where hLBo =LGr > 5 (approximate). Results of our growth and departure model calculated for contact angles up to 90° and DT ¼ 20 C are shown in Fig. 5 along with the boundary where hLBo =LGr ¼ 5 which delineates the limit of our model’s applicability. When hLBo =LGr < 5, we recommend using a static departure description such as the spherical cap departure criterion (Eq. (21)) of Fritz correlation (Eq. (2)). According to Fig. 5b,c, a static departure description is valid for very hydrophilic surfaces (h < 20 ), which is associated with small bubble departure sizes. As shown in Fig. 4, modeled values for Dequiv;depart agree with previous studies within 4:1 to 5:6% (2r). Modeled values for tdepart agree with previous studies within 16 to 28% (2r). With the ability to determine bubble departure and model growth, it is possible to calculate the heat transfer rate associated with a single nucleation site (Fig. 5c): 4
qsingle ¼ qv hfg 3
p Dequiv;depart =2 t depart þ twait
3 ð23Þ
where twait is the waiting time between bubble departure and bubble growth.
Fig. 4. Our bubble growth and departure model for saturated water at atmospheric conditions is in good agreement with experimentally-verified simulation results from the literature. Bubble growth is modeled according to Eq. (15) and departure is determined as described in Section 4. Comparisons are shown for (a) ðh ¼ 54 ; DT ¼ 10:0 CÞ [4], (b) ðh ¼ 38 ; DT ¼ 12:5 CÞ [4], (c) ðh ¼ 38 ; DT ¼ 6:2 CÞ [4], (d) ðh ¼ 50 ; DT ¼ 7:0 CÞ [4], (e) ðh ¼ 25 ; DT ¼ 8:5 CÞ [4], and (f) ðh ¼ 30 ; DT ¼ 8:5 CÞ [4], (g) ðh ¼ 38 ; DT ¼ 8:5 CÞ [4], (h) ðh ¼ 50 ; DT ¼ 8:5 CÞ [5], and (i) ðh ¼ 10 ; DT ¼ 5:3 CÞ [6].
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
7
Fig. 5. Bubble departure results from our model for saturated water at atmospheric pressure. (a) Departure equivalent diameter is shown as a function of superheat (colors) and contact angle. Departure equivalent diameter is determined from bubble size (Eq. (15)) at the moment of bubble departure. How the bubble departs depends on hLBo =LGr . If hLBo =LGr > 5, the momentum boundary layer is important and departure is as described in Section 4—a simple approximation of Dequiv;depart is given by Eq. (21). Conversely, if hLBo =LGr < 5, then departure is described by a static departure description (surface tension vs buoyancy) such as the spherical cap departure criterion (Eq. (21)) or the Fritz correlation (Eq. (2)). The static descriptions of departure size have no dependence on superheat (colors). For both momentum and static regimes (spherical cap criterion), (b) departure time and (c) single bubble heat transfer are plotted as a function of superheat and contact angle. Single bubble heat transfer takes into account latent heat, departure size, departure time, and waiting time (Eqs. (23) and (24)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
To calculate the waiting time accurately, one should take into account superheat reduction effects during the bubble growth stage, which could include microlayer evaporation effects. For instance, microlayer evaporation could lower superheat within the growth period, effectively increasing the waiting time for a bubble. This effect could also be highly dependent on the Jakob number. For simplicity, we have used an approximate twait assuming semi-infinite conduction into the thermal boundary layer [31,34] and a thermal boundary layer as dictated by the Rayleigh number:
t wait ¼
d2T
pal
¼
m2=3
: 1=3
ð24Þ
g 2=3 bl2=3 DT 2=3 pal
From Fig. 5c, it is apparent that lower contact angles have lower heat transfer. Thus, if nucleation density were constant for a given surface, higher contact angles would be favorable for better heat transfer. Using the results of our growth and departure models, the departure diameter has been correlated to h and DT as
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dequiv;depart ¼ 0:652h0:581 DT 0:313
clv
g ðql qv Þ
ð25Þ
with R2 ¼ 0:999 for both our model results as well as simulation results [4,5].
It should be noted that our departure model assumes natural convection boundary layers, which is completely reasonable for an isolated bubble. However, at high nucleation densities, bubbles will start to interact with each other and the natural convection assumption will break down. This could be addressed using a different boundary layer correlation such as the one suggested by Tien [35]:
dT ¼ 3:22
kl hboil
ð26Þ
where hboil is the heat transfer coefficient in boiling and kl is the liquid conductivity. It should also be noted that we have assumed that wetting properties are constant. However, our departure model, like the growth model, can also be extended to cases where wetting properties are time-dependent (Section S4, Fig. S6). 5. Conclusions Studying the bubble cycle of nucleation, growth, and departure is fundamental and necessary to gain a more complete understanding of nucleate pool boiling. We have found that nucleation is linked to bubble departure, which suggests that nucleation criteria is dependent on wettability and convection. We have developed a growth model that assumes a bubble grows within a thermal boundary layer as opposed to the uniform superheating assumed
8
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
in previous theories. We theorized a new mechanism of bubble departure wherein a velocity boundary layer is responsible for the rewetting underneath the bubble. Both the growth and departure models agree well with experimentally verified numerical simulation studies of single bubbles. By tracking the movement of the bubble base, we have also developed a methodology to calculate growth and departure for contact angle changing with time, which may be due to changes in surface tension and/or contact angle. Through the use of a simplified model, the results of this work provide a more physical understanding of the boiling process and simple quantitative tools to model boiling heat transfer. To improve upon this understanding further, future work could incorporate more detailed modeling of dynamic contact angle and microlayer evaporation, as well as explore non-aqueous fluids, mixtures, and nanofluids. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This work has been funded by the Singapore-MIT Alliance for Research and Technology (SMART). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer. 2019.118803. References [1] R. Cole, Bubble frequencies and departure volumes at subatmospheric pressures, AIChE J. 13 (4) (1967) 779–783, https://doi.org/10.1002/ aic.690130434. http://doi.wiley.com/10.1002/aic.690130434 . [2] B. Mikic, W. Rohsenow, P. Griffith, On bubble growth rates, Int. J. Heat Mass Transf. 13 (4) (1970) 657–666, https://doi.org/10.1016/0017-9310(70)900402. http://linkinghub.elsevier.com/retrieve/pii/0017931070900402 . [3] S. van Stralen, M. Sohal, R. Cole, W. Sluyter, Bubble growth rates in pure and binary systems: combined effect of relaxation and evaporation microlayers, Int. J. Heat Mass Transf. 18 (3) (1975) 453–467, https://doi.org/10.1016/0017-9310 (75)90033-2. http://linkinghub.elsevier.com/retrieve/pii/0017931075900332 . [4] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface, J. Heat Transf. 121 (3) (1999) 623, https://doi.org/10.1115/1.2826025. http://heattransfer. asmedigitalcollection.asme.org/article.aspx?articleid=1443768 . [5] A. Mukherjee, S.G. Kandlikar, Numerical study of single bubbles with dynamic contact angle during nucleate pool boiling, Int. J. Heat Mass Transf. 50 (1–2) (2007) 127–138, https://doi.org/10.1016/j.ijheatmasstransfer.2006.06.037. http://linkinghub.elsevier.com/retrieve/pii/S0017931006004157 . [6] Y. Nam, E. Aktinol, V.K. Dhir, Y.S. Ju, Single bubble dynamics on a superhydrophilic surface with artificial nucleation sites, Int. J. Heat Mass Transf. 54 (7–8) (2011) 1572–1577, https://doi.org/10.1016/j.ijheatmasstransfer.2010.11.031. http:// linkinghub.elsevier.com/retrieve/pii/S0017931010006459 . [7] J. Kim, Review of nucleate pool boiling bubble heat transfer mechanisms, Int. J. Multiph. Flow 35 (12) (2009) 1067–1076, https://doi.org/10.1016/j. ijmultiphaseflow.2009.07.008. http://linkinghub.elsevier.com/retrieve/pii/ S0301932209001311 . [8] N.I. Kolev, How Accurately Can We Predict Nucleate Boiling?, Springer, Berlin, Heidelberg, 2012, pp 143–178, https://doi.org/10.1007/978-3-642-21372-4_6. [9] H.-J.J. Cho, D.J. Preston, Y. Zhu, E.N. Wang, Nanoengineered materials for liquid–vapour phase-change heat transfer, Nat. Rev. Mater. 2 (2) (2016) 16092, https://doi.org/10.1038/natrevmats.2016.92. http://www.nature.com/ articles/ natrevmats201692 . [10] S.G. Bankoff, Entrapment of gas in the spreading of a liquid over a rough surface, AIChE J. 4 (1) (1958) 24–26, https://doi.org/10.1002/aic.690040105. http://doi.wiley.com/10.1002/aic.690040105 . [11] P. Griffith, J.D. Wallis, The Role of Surface Conditions in Nucleate Boiling, Tech. rep., Massachusetts Institute of Technology, Cambridge, MA, 1958.
[12] J.J. Lorenz, The Effects of Surface Conditions on Boiling Characteristics (Ph.D. thesis), Massachusetts Institute of Technology, 1972. [13] Y.Y. Hsu, On the size range of active nucleation cavities on a heating surface, J. Heat Transf. 84 (3) (1962) 207, https://doi.org/10.1115/1.3684339. http:// heattransfer.asmedigitalcollection.asme.org/article.aspx?articleid=1432587 . [14] M.S. Plesset, S.A. Zwick, The growth of vapor bubbles in superheated liquids, J. Appl. Phys. 25 (4) (1954) 493–500, https://doi.org/10.1063/1.1721668. http:// aip.scitation.org/doi/10.1063/1.1721668 . [15] S.J.D. van Stralen, R. Cole, W.M. Sluyter, M.S. Sohal, Bubble growth rates in nucleate boiling of water at subatmospheric pressures, Int. J. Heat Mass Transf. 18 (5) (1975) 655–669, https://doi.org/10.1016/0017-9310(75)90277-X. [16] V.P. Carey, Liquid Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation Processes in Heat Transfer Equipment, second ed., Taylor & Francis, 2007. [17] P. Stephan, T. Fuchs, E. Wagner, N. Schweizer, Transient local heat fluxes during the entire vapor bubble life time, in: ECI International Conference on Boiling Heat Transfer, Brazil, 2009. [18] T. Yabuki, O. Nakabeppu, Heat transfer mechanisms in isolated bubble boiling of water observed with MEMS sensor, Int. J. Heat Mass Transf. 76 (C) (2014) 286–297. [19] W. Fritz, Maximum volume of vapor bubbles, Phys. Z. 36 (11) (1935) 379–384. [20] L. Cheng, D. Mewes, A. Luke, Boiling phenomena with surfactants and polymeric additives: a state-of-the-art review, Int. J. Heat Mass Transf. 50 (13–14) (2007) 2744–2771. [21] A. Zou, A. Chanana, A. Agrawal, P.C. Wayner, S.C. Maroo, Steady state vapor bubble in pool boiling, Sci. Rep. 6 (1) (2016) 20240, https://doi.org/10.1038/ srep20240. http://www.nature.com/articles/srep20240 . [22] J.N. Addoms, Heat Transfer at High Rates to Water Boiling Outside Cylinders (Ph.D. thesis), Massachusetts Institute of Technology, 1948. [23] A. Sloan, S. Penley, R.A. Wirtz, Sub-atmospheric pressure pool boiling of water on a screen-laminate enhanced surface, 2009 25th Annual IEEE Semiconductor Thermal Measurement and Management Symposium, vol. 50, IEEE, 2009, pp. 246–253, https://doi.org/10.1109/STHERM.2009.4810771. http://ieeexplore. ieee.org/document/4810771/ . [24] S.M. Kwark, M. Amaya, R. Kumar, G. Moreno, S.M. You, Effects of pressure, orientation, and heater size on pool boiling of water with nanocoated heaters, Int. J. Heat Mass Transf. 53 (23–24) (2010) 5199–5208, https://doi.org/ 10.1016/j.ijheatmasstransfer.2010.07.040. http://linkinghub.elsevier.com/ retrieve/pii/S0017931010004187 . [25] X.-F. Yang, Z.-H. Liu, Pool boiling heat transfer of functionalized nanofluid under sub-atmospheric pressures, Int. J. Therm. Sci. 50 (12) (2011) 2402–2412, https://doi.org/10.1016/j.ijthermalsci.2011.07.009. http://linkinghub. elsevier.com/retrieve/pii/S1290072911002213 . [26] W. Tong, A. Bar-Cohen, T. Simon, S. You, Contact angle effects on boiling incipience of highly-wetting liquids, Int. J. Heat Mass Transf. 33 (1) (1990) 91– 103, https://doi.org/10.1016/0017-9310(90)90144-J. http://www.sciencedirect. com/science/article/pii/001793109090144J . [27] W. Wagner, A. Pruß, The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use, J. Phys. Chem. Ref. Data 31 (2) (2002) 387–535, https://doi.org/10.1063/1.1461829. http://aip.scitation.org/doi/10.1063/1.1461829 . [28] C.H. Wang, V.K. Dhir, Effect of surface wettability on active nucleation site density during pool boiling of water on a vertical surface, J. Heat Transf. 115 (1993) 659, https://doi.org/10.1115/1.2910737. [29] N. Siemons, H. Bruining, H. Castelijns, K.-H. Wolf, Pressure dependence of the contact angle in a CO2–H2O–coal system, J. Colloid Interface Sci. 297 (2) (2006) 755–761. [30] J.B. Goodman, N.W. Krase, Solubility of nitrogen in water at high pressures and temperatures, Ind. Eng. Chem. 23 (4) (1931) 401–404. [31] H. Chi-Yeh, P. Griffith, The mechanism of heat transfer in nucleate pool boiling— Part I, Int. J. Heat Mass Transf. 8 (6) (1965) 887–904, https://doi.org/10.1016/00179310(65)90073-6. http://linkinghub.elsevier.com/retrieve/pii/0017931065900748 http://linkinghub.elsevier.com/retrieve/pii/0017931065900736 . [32] R. Cole, W.M. Rohsenow, Correlation of bubble departure diameters for boiling of saturated liquids, Chem. Eng. Prog. Symp. Ser. (1969). [33] R. Raj, C. Kunkelmann, P. Stephan, J. Plawsky, J. Kim, Contact line behavior for a highly wetting fluid under superheated conditions, Int. J. Heat Mass Transf. 55 (9–10) (2012) 2664–2675. [34] M.Z. Podowski, R.M. Podowski, Mechanistic Multidimensional modeling of forced convection boiling heat transfer, Sci. Technol. Nucl. Instal. 2009 (2009) 1–10. [35] C.L. Tien, A hydrodynamic model for nucleate pool boiling, Int. J. Heat Mass Transf. 5 (6) (1962) 533–540. [36] H. Fedders, Messung des Wärmeüberganges beim Blasensieden von Wasser an metallischen Rohren, no. Jül-740 RB, Kernforschungsanlage Jülich, 1971. [37] J.V. Sengers, J.M.H.L. Sengers, Thermodynamic behavior of fluids near the critical point, Annu. Rev. Phys. Chem. 37 (1) (1986) 189–222, https://doi.org/ 10.1146/annurev.pc.37.100186.001201.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Bubble nucleation, growth, and departure: A new, dynamic understanding H. Jeremy Cho a,b, Evelyn N. Wang a,⇑ a b
Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA, United States University of Nevada, Las Vegas, 4505 S Maryland Pkwy, Las Vegas, NV, United States
a r t i c l e
i n f o
Article history: Received 11 February 2019 Received in revised form 23 September 2019 Accepted 27 September 2019
Keywords: Boiling Bubble Nucleation Growth Departure Wetting
a b s t r a c t The bubble nucleation, growth, and departure cycle is a fundamental aspect of nucleate pool boiling. While much research on this subject has been performed in the previous century, new correlations and models have not been developed in light of important results in simulations and experiments in recent decades. In this work, we provide an updated understanding of nucleation, growth, and departure with analytical models that are validated by recent work. We found that nucleation at incipience is correlated to departure size, suggesting that the momentum-induced bubble snap-off at departure may be responsible for the size of trapped vapors in cavities. We also developed a bubble growth model that takes into account a thermal boundary layer whereas previous works typically considered a uniform superheating. In addition, we developed a bubble departure description, wherein the velocity boundary layer is responsible for closing the base of the bubble. Finally, we provide a methodology to calculate bubble growth and departure for contact angle changing with time, which is relevant in studying surfactantenhanced boiling. Ó 2019 Published by Elsevier Ltd.
1. Introduction Boiling is utilized in a variety of industrial and domestic applications. Despite more than half a century of research, many uncertainties about the process still remain. Boiling is a highly complex process involving such phenomena as turbulent flows, large temperature gradients, liquid–vapor phase change, and contact line dynamics. To gain a simplified understanding, much of the research on boiling fundamentals has focused on the bubble cycle: the subsequent events of bubble nucleation on a boiling surface, bubble growth, and bubble departure from the surface. A single bubble cycle transfers an amount of heat per time qsingle ; knowing how to calculate qsingle allows one to calculate heat transfer for the entire boiling surface. The heat transfer during boiling is nqsingle where n is the nucleation site density. While many theories, models, and correlations were developed from the 1950s to 1970s that allow for the calculation of qsingle [1–3], the advent of computer simulations and more sophisticated experiments in recent decades, have revealed a richer physical picture of the bubble cycle and accurate bubble cycle data [4–6]. In particular, it ⇑ Corresponding author. E-mail address: [email protected] (E.N. Wang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118803 0017-9310/Ó 2019 Published by Elsevier Ltd.
has been shown that heat transfer primarily occurs through the liquid–vapor interface either through the microlayer or from the bulk surrounding liquid [7]. In light of these new experiments and simulations, we present a way to determine qsingle that is based on new physical insights and validated by the latest research. We do not discuss how to calculate nucleation site density since it is highly surface-dependent [8] and is often a feature that is designed for in modern boiling surfaces [9].
1.1. Background on nucleation The first stage of the bubble cycle is bubble nucleation. According to trapped vapor theory (also known as entrapped vapor theory), cavities on the surface trap air bubbles; these trapping sites serve as nucleation sites in boiling [10]. Because of the convex curvature of trapped vapor bubbles at cavities, the vapor pressure inside the bubble is higher than the surrounding liquid. As such, if the surrounding liquid is at the saturation temperature, its equilibrium vapor pressure would not match the vapor pressure inside the trapped bubble. However, if the liquid were sufficiently superheated, its equilibrium vapor pressure would be greater than the bubble vapor pressure, causing the bubble to grow. The particular condition when this bubble growth—or nucleation—occurs is when
2
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
Nomenclature
al
liquid thermal diffusivity molar gas constant velocity boundary layer clv liquid-vapor surface tension l chemical potential ql liquid density qv vapor density qc T Ja modified Jakob number ( lqp;lh sat ) v fg ReL Reynolds number h contact angle (radians) A area abase acceleration of base contact line Aq surface area where heat transfer occurs on growing bubble C trans transition time scaling factor Cb growth rate scaling factor cp;l liquid specific heat base diameter Dbase Dequiv;depart bubble equivalent diameter at departure Dequiv bubble equivalent diameter g acceleration due to gravity hboil boiling heat transfer coefficient R dv
the chemical potential of the liquid matches the chemical potential of the vapor. This is often approximated by combining the Clausius–Clapeyron relation with the Young–Laplace equation [11]:
DT
2clv T sat v fg : hfg Rb
ð1Þ
Thus, if the radius of curvature, Rb , of a nucleation site is known, then the superheat at which bubble nucleation occurs can be calculated using Eq. (1), which we denote as the Clausius–Clapeyron-Y oung–Laplace (CCYL) equation. Two common ways to obtain the bubble radius of curvature, Rb , is through the Lorenz and Hsu models [12,13]. However, we have found that these models do not agree with experiments of nucleation activation at various pressures. Thus, we provide a simple correlation for Rb based on bubble departure size that more accurately reflects these experiments. 1.2. Background on growth The next stage of the bubble cycle is bubble growth. The growth rate depends on wetting behavior and superheat. Plesset and Zwick [14] modeled this growth assuming that heat transfer from a superheated liquid to a bubble is diffusion-limited. In this diffusion-limited case, Rb / t 1=2 . Mikic et al. further extended this model to include an inertia-dominated growth regime where Rb / t at very short timescales before transitioning to the diffusion-limited growth regime where Rb / t 1=2 [2]. van Stralen et al. considered growth with a thin evaporation microlayer underneath the bubble and also found that Rb / t 1=2 [15]; thus, it is difficult to differentiate between microlayer and diffusion-limited growth [16]. Nonetheless, many others have investigated the relative heat transfer contribution from the microlayer [7,17,18]. Regardless of the microlayer contribution, heat transfer involves conduction through a liquid layer in both microlayer and bulk DT ffi nature is regions. Thus, heat transfer of a semi-infinite, q00 pklffiffiffiffiffiffi pal t
present around the bubble, leading to a Rb / t 1=2 dependence. Recently, experimentally-validated computer simulations of single-bubble growth that take into account microlayer effects
hfg kl LBo LGr n n P Psat qsingle R2 T t td T sat t trans t wait V
v fg
z
specific enthalpy of liquid-vapor phase change liquid conductivity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi characteristic Bond length scale ( gðqclv q Þ) l v characteristic Grashof length moles of ideal gas nucleation density pressure saturation pressure single bubble heat transfer coefficient of determination temperature time time at which bubble growth transitions to t 1=4 growth saturation temperature transition time when constant base closing force starts bubble wait time volume specific volume of liquid-vapor phase change dimension normal to surface
and macro-scale effects indicate that the Rb / t1=2 dependence breaks down at longer timescales closer to departure [4–6]. In this work, we took into account the thermal boundary layer to develop an analytical model of bubble growth that agrees very well with experimentally-validated computer simulations. We consider semi-infinite conduction through the liquid around the bubble— assumed to take a pendant bubble geometry—as a means to lump the effects of bulk and microlayer heat transfer together since they lead to similar growth signatures. Our model shows an initial diffusion-limited Rb / t1=2 growth characteristic of bulk and microlayer heat transfer, and an eventual Rb / t 1=4 at longer timescales.
1.3. Background on departure The final stage of the bubble cycle is bubble departure. The departure of bubbles from the boiling surface is the result of an interaction of buoyancy, momentum, and surface tension forces. Most correlations are based off of a static force balance between surface tension and buoyancy where the departure diameter is approximately on the order of the characteristic Bond length: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi LBo clv =ððql qv Þg Þ. Fritz developed a simple correlation [19] where
Dequiv;depart ¼ 1:19hLBo :
ð2Þ
Later correlations incorporated wall superheat effects but did not incorporate contact angle effects [1]. Currently, there are no departure models that consider the time-dependent contact angle behavior that are present due to dynamic contact line movement, surfactants [20], or microlayer evaporation [21]. To address this limitation, rather than develop a generic bubble departure size correlation, we provide a method to calculate the time of departure based on the movement of the contact line taking into account both superheat and contact angle—including dynamic contact angles. Our departure description incorporates effects of the velocity boundary layer and agrees well with experimentally validated numerical simulations of single bubbles [4–6].
3
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
2. Bubble nucleation
equations are solved simultaneously to determine Rb ; ql , and qv from a known DT; T sat , and P sat :
2.1. Limitations of existing models
lIAPWS ðT sat þ DT; ql Þ ¼ lIAPWS ðT sat þ DT; qv Þ; PIAPWS ðT sat þ DT; ql Þ ¼ Psat 2 c ð T þ DT Þ PIAPWS ðT sat þ DT; qv Þ ¼ Psat þ lv :
To evaluate the limitations of existing nucleation models, we compared the models’ predictions of radii of curvature, Rb , to experimentally derived values as a function of pressure. Boiling studies at various pressures show that activation superheats increase with decreasing pressure [22–25] (Fig. 1). From boiling curve data, it is possible to determine Rb from the superheat at which nucleation occurs. This can be done using a ClausiusClapeyron-Young-Laplace (CCYL) equation (Eq. (1)). However, since the CCYL approximation is only valid at small superheats [26], we chose to perform an exact calculation of equivalent potentials. This calculation is possible with steam tables. We have performed our analysis with water properties as given by the International Association for the Properties of Water and Steam (IAPWS) [27]. Here, we equate chemical potentials (or Gibbs free energies) between the vapor and liquid states, and relate pressure to bubble curvature, Rb , using the Young–Laplace equation in an approach we denote as the chemical-equilibrium-Young–Laplace (CEYL) method. To apply the CEYL approach, the following
Fig. 1. Calculated bubble radii of curvature as a function of pressure from experiments [22–25,36] show that Rb scales with the bubble departure size (blue curve) instead of advancing-front models—red dashed curve for constant volume (Eq. S3) [12] and red curve for constant Helmholtz energy (Eq. S12). We have taken data from the literature and estimated the temperature at incipient boiling from the boiling curve using a least-squares fitting, then applied Eqs. (3)–(5) to find Rb;incipient . It is evident that Rb;incipient decreases with increasing pressure. This is contrary to the advancing-front description, which predicts either an invariance of Rb;incipient with pressure (constant volume model) or a decrease in Rb;incipient with decreasing pressure (constant Helmholtz model). Thus, there must be some physical mechanism for vapor trapping outside of simple wetting in an advancing liquid front. The advancing-front models plotted uses a conical cavity with / ¼ 25 ; h ¼ 50 , and 1=3 Rcav ¼ 5 lm. The Hsu model [13] plotted has C Hsu ¼ 1:6 and dT ¼ 0:5 gbmaDlT and predicts the opposite trend with pressure. The idea of trapped gases compressing due to increased saturation pressure is an insufficient predictor of incipience as modeled using the ideal gas law with capillary pressure (Eq. (6), purple curve). The critical pressure of water, P crit , is 22.1 MPa. Our model (blue curve) is likely invalid very close to the critical point since the departure correlation (Eq. (9)) diverges due to an infinite specific heat [37]; however, we show agreement with experimental results up to 150 MPa. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Rb
ð3Þ ð4Þ ð5Þ
The incipient superheat, DT, used for plotting literature data in Fig. 1 was estimated by performing a linear regression to the region of the boiling curve data where a steep rise occurs; the intercept with the superheat axis serves as an approximate incipient superheat. Eqs. (3)–(5) can also be used by any other fluid as long as the equation of state (steam table equations) are known. Since solving these equations is computationally inexpensive, we recommend using the CEYL method whenever possible, especially at higher superheats. 2.1.1. Advancing-front models Lorenz developed a geometric model of the trapped vapor geometry assuming a quasi-static advancing liquid front over a conical cavity [12]. This description depends entirely on contact angle and cavity geometry. Tong et al. improved upon Lorenz’ model by recognizing that the dynamic contact angle should be involved in this cavity-wetting behavior [26]. Some relaxation assumption about the bubble geometry must be made as the bubble relaxes from the moment it is trapped (state 1 in Fig. S1) to the embryo state where a trapped bubble of curvature Rb is formed (state 2 in Fig. S1). In Lorenz’ model (Section S1), a constant volume between state 1 and state 2 is assumed (V 1 ¼ V 2 ); thus, the radius of curvature can be related to cavity geometry and contact angle, h (Eq. S3). Assuming that the surrounding liquid is uniformly at a superheated temperature, the CCYL relation (Eq. (1)) can be used to relate the superheat of nucleation activation, DT, to cavity geometry and contact angle (Eq. S4). The result is that nucleation activation superheat decreases with increasing contact angle—a typical observation in experiments [28]. However, Lorenz’ description does not capture any dependence on pressure (Fig. 1, red dashed line)—we did not include pressure-dependent contact angle since this effect is likely minimal [29]. Even using a different relaxation assumption that is more physically reasonable, where V 1 – V 2 and the Helmholtz free energy is conserved (isothermal condition, which is derived in Section S1), results in a pressure dependence that is completely counter to the experimental data (Fig. 1, red solid line). Therefore, there is little support for the quasi-static advancing front description with uniformly superheated liquid from literature data. One possible explanation for the lowering of incipient radii at higher pressures could be the compression of vapor regions that were formed during cavity wetting at lower pressures. To investigate this possible mechanism, we modeled the compression of an ideal gas with elevated bubble pressure from capillarity:
2c ðT sat ; P sat Þ 4 3 Psat þ lv pRb ¼ nRT: Rb 3
ð6Þ
After solving Eq. (6) for Rb using IAPWS properties across different saturation conditions, it is clear that simple compression of trapped vapor regions underpredicts the range of incipient radii from experimental data (Fig. 1, purple curve). Another possible explanation for this decrease in incipient radii with pressure could be the dissolution of gases into water at elevated pressures. Dissolution of gases in water follows Henry’s Law, where solubility is proportional to pressure [30]. At saturation conditions, this proportionality still holds since the Henry’s Law constant does not vary much with temperature [30]. Thus, at elevated pressures, it is likely that bubbles would become completely dissolved since solubilities would be
4
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
increased many times over. Regardless, the dissolution of bubbles should be highly system-specific since the relative amounts of water and entrapped vapor can vary, and hence the amount of gas dissolution can vary. From the experimental data plotted in Fig. 1, however, it is clear that bubble incipience follows a general trend across independent data sources; thus, the effects of solubility, if any, appear to be minimal. 2.1.2. The Hsu bubble viability model Hsu [13] and others [31] considered a growing bubble within a quiescent superheated thermal boundary layer. After every departure event, a one-dimensional, transient conduction occurs within the thermal boundary layer with thickness dT . Hsu showed that if the liquid within this layer does not sufficiently heat the trapped bubble embryo according to the CCYL relation, then the site remains inactive. Mathematically, there is a discriminant (quantity within a square-root) that must be positive for a solution of a bubble actively growing to exist. This condition, which we denote as bubble viability, is
DT ¼
8clv T sat ð1 þ cos ðc þ hÞÞ hfg qv dT
ð7Þ
when at saturation conditions. Combining the Hsu superheat criterion (Eq. (7)) with the CEYL equations (Eqs. (3)–(5)) yields a relationship between DT and Rb as a function of pressure (Fig. 1, green curve). It is clear that the Hsu criterion is not consistent with experimental results as it predicts a decrease in Rb with a decrease in pressure. Thus, there is little experimental support for the description of quiescent, one-dimensional, transient conduction in the thermal boundary layer driving bubble nucleation. 2.2. New correlation for the bubble radius of curvature at nucleation In contrast to the advancing-front and Hsu bubble viability models, we hypothesized that convective motion may play a role in bubble nucleation. More convective motion implies that bubbles depart earlier (smaller bubble departure diameter, Ddepart ) and that the velocity and thermal boundary layer are thinner. With thinner thermal boundary layers, trapped vapor embryos would be exposed to higher superheats, leading to smaller Rb at nucleation. Thus, to test the link between convection and nucleation, we examined whether Rb decreases with decreasing bubble departure size. Using the departure diameter correlation by Cole and Rohsenow [32] where
Ddepart ¼ 1:5 104
clv
ðql qv Þg
1=2
Ja5=4
ð8Þ
and
Ja
ql cp;l T sat ; qv hfg
ð9Þ
is the modified Jakob number to describe the state of saturation, it was found that Rb / Ddepart : 3
Rb 1:6 10 Ddepart :
ð10Þ
The agreement with experimental data over several decades of pressure measurement (Fig. 1, blue line) confirms our hypothesis that increased convection causes smaller bubble embryos to nucleate, perhaps through a snap-off mechanism. Thus, using Eq. (10), a correlation for bubble departure (a new departure correlation is presented in Section 4), and the CCYL equation (Eq. (1)) or CEYL equations (Eqs. (3)–(5)), it is possible to calculate the superheat at which a nucleation site will activate.
3. Bubble growth Given the importance of convection as suggested by our nucleation description, we modeled bubble growth in the presence of velocity and thermal boundary layers. In our simplified modeling approach, we assume that growing bubbles can be approximated by the static pendant bubble shape in an unchanging thermal and boundary environment, the properties of which are determined by the superheat. We note that these assumptions and simplifications are not completely physical; however, we validate this simplified description through agreement with available data on bubble growth from literature sources. To quantify bubble size, we use equivalent radii and diameters, 1=3 Requiv and Dequiv , respectively: Requiv ¼ Dequiv =2 ¼ 43p V . As Mikic et al. have found [2], bubble growth can be described by two regimes. At short timescales, there is inertially dominated growth in which Requiv / t and is described by the Rayleigh equation. At longer timescales, growth is heat diffusion limited in which Requiv / t1=2 , consistent with Plesset and Zwick’s description [14]. Since the timescale at which this transition occurs is less than one millisecond whereas bubble departure times are typically in the tens of milliseconds, we assumed that growth is only diffusion-limited. While Plesset and Zwick and Mikic et al. did not consider contributions from microlayer evaporation, van Stralen et al. showed that microlayer growth also has a t 1=2 dependence [3]. Thus, growth from the microlayer should mathematically resemble the diffusion-limited growth as described by Plesset and Zwick [16]. As a result, in our model, any microlayer effects are effectively taken into account by a prefactor, C b . Plesset and Zwick stated that the conductive heat transfer into and around a growing spherical bubble is equivalent to the latent heat times the mass change rate of the bubble:
kl DT dV : Alv pffiffiffiffiffiffi qv hfg dt al t
ð11Þ
From this relationship, we can derive an expression for the bubble growth rate as a function of fluid properties and contact angle for a bubble with a spherical cap geometry:
Requiv ðt Þ h 1=3 ¼ Rðt Þ 4 2 cos 2 cos4 2h cos ðhÞ
pffiffi 8kl DT pffiffiffiffiffi t : ð3 þ 2 cos ðhÞ cos ð2hÞÞqv hfg al
ð12Þ
Though the spherical cap geometry is simple, it is unphysical as it neglects the effects of gravity. Bubbles adhered to a horizontal wall in the presence of gravity will have a pendant geometry (Section S3). Furthermore, the pendant bubbles would be in the presence of a thermal boundary layer. Thus, we developed a growth model that assumes a pendant bubble in the presence of a thermal boundary layer. As shown in Fig. 2, we have used a simple thermal boundary layer model where the liquid temperature linearly decreases from T ¼ T wall at z ¼ 0 (at the solid-liquid interface) to T ¼ T sat at z ¼ dT so
z DT ¼ DT wall 1 dT
ð13Þ
where dT is the thermal boundary length. We also assume that the temperature of the vapor bubble is at T sat . Use a scaling analogous to Eq. (11), but with the incorporation of the thermal boundary layer,
Z Aq
kl ðT ðzÞ T sat Þ p dDequiv pffiffiffiffiffiffi dA D2equiv qv hfg 2 dt al t
ð14Þ
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
5
more closely matches this flattening due to its transition to a Dequiv / t 1=4 dependence. Furthermore, our growth model can be extended to cases when surface tension is time-dependent (Section S3). 4. Bubble departure In our treatment of bubble departure, we track the movement of the three-phase contact line with time based on our growth model. To do so, we relate the diameter of the three-phase contact circle, Dbase , to the bubble size:
Dbase ¼ 0:8hDequiv
Fig. 2. Heat transfer around a growing bubble within a thermal boundary layer. A bubble with height lb having a pendant shape grows as a result of lateral heat conduction, q00 , at its bottom. The heat flux is larger near the surface and decreases to zero at the thermal boundary layer height, dT . The temperature profile in the thermal boundary layer is linear.
where dA is a differential bubble surface area at a height z and Aq is the surface area where heat transfer occurs (any bubble surface that is within the thermal boundary layer). By incorporating the pendant bubble geometry (Section S3), the growth rate can be solved analytically:
8 pffi 2 > < 2dc T 1 exp C b bc DdT wall t 0 6 t < td T Requiv ðt Þ ¼ pffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi > : 2 C bdT DT pffiffit d2 ðln 41Þ t P t d b wall 2c2
ð15Þ
where
h c ¼ 2 cos ; 2 td ¼
dT ln 2 2
C b bc DT wall
ð16Þ !2 :
ð17Þ
Here, we assume the contact angle, h, to be the static contact angle. Using the dynamic contact angle could lead to more accurate modeling. However, the use of the static contact angle in this case is reasonable since we are considering the growth of isolated, single-bubbles, which tends to occur at lower superheats where dynamic contact angle effects are small [33]. We also provide a methodology for incorporating dynamic contact angles in our Supplementary Information (S4). When t < td , the height of the bubble, lb is less than the thermal boundary thickness, dT , and when t > t d ; lb > dT . Here, Eq. (15) has the behavior of t1=2 growth when the bubble is much smaller than the boundary length (same result as Plesset and Zwick where a uniform superheating is assumed) and t1=4 growth when the bubble is larger than the boundary length. By performing a least-squares fit of Eq. (15) to nine data sets of experimentally verified single-bubble simulations from several studies [4–6] (Fig. 4), we have found that the scaling factor C b ¼ 0:534 is appropriate for saturated water at atmospheric conditions. For comparison, we also fit the bubble growth model for a spherical cap under uniform superheating (Eq. (12)) to the same data sets. We found that the agreement of our pendant bubble in
ð18Þ
where h is radians (applicable for a variety of Bond numbers). This approximation is close to the actual Young-Laplace solution (Fig. S2). As suggested by our description of bubble nucleation, we hypothesized that convective forces determine the point at which bubble departure occurs. Closer observation of the base diameter velocity as a function of time for a bubble simulated by Mukherjee and Kandlikar [5] as shown in Fig. 3 suggests that a momentum force may be responsible for bubble departure. Initially, the base diameter velocity is high and positive as the bubble grows rapidly. However, the rate at which this base radial velocity changes is not constant (inward radial acceleration decreases). At some transition time, ttrans , the radial velocity of the contact line decreases linearly with time (constant inward radial acceleration). This constant inward radial acceleration is evidence of a convective force, consistent with our hypothesis. Thus, we modeled the evolution of Dbase with time such that, initially, Dbase increases according to pendant bubble geometry (Eq. (18)). After the transition time, t trans , the acceleration of the base closing, abase , is maintained at its value at ttrans (second derivative of Dbase with respect to time). This causes Dbase to decrease until the point of departure when Dbase ¼ 0. Until the point of departure, growth proceeds according to Eq. (15). We can determine ttrans from a velocity boundary layer scaling analysis. According to the boundary layer thickness equation, pffiffiffiffiffiffiffi dv L= ReL where dv is the velocity boundary layer thickness and L is the boundary layer length. Using ReL ¼ L2 =ðt ml Þ, the characteristic timescale associated with boundary layer development is d2v =m. Assuming ttrans is on the order of this time scale, ttrans d2v =m. Assuming the velocity boundary layer thickness
boundary layer model (Eq. (15)) was superior (R2 ¼ 0:994) to that of spherical cap bubble in uniform superheat model (Eq. (12)) (R2 ¼ 0:985). The Dequiv / t1=2 growth associated with uniform superheating does not closely follow the flattening of the growth curve at longer times (Fig. 4). On the other hand, our growth model
Fig. 3. Base diameter velocity shows a transition to a constant acceleration zone. Data of bubble growth and departure with DT wall ¼ 10 C and h ¼ 54 from a numerical simulation by Mukherjee and Dhir [5] shows the base diameter velocity transitioning to a constant acceleration zone.
6
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
(dictated by the Grashof number) can be modified by the contact angle (smaller contact angles should have thinner layers), then 1=3 2 . The resulting ttrans is then dv h gbm DT l
t trans ¼ C trans h2
m1=3
ð19Þ
ðgbl DT Þ2=3
where
C trans ¼ 0:261
ð20Þ
was found to be an appropriate scaling factor for saturated water at atmospheric pressure. However, the analysis above is in the limit where convection is important. In the limit of small bubbles (very hydrophilic surfaces; low h), bubble departure may simply result from a static force balance between buoyancy and surface tension:
Dequiv;depart
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 24 sin h clv ¼ : 3 2 þ 3 cos h g ðql qv Þ
ð21Þ
One way to determine the relative importance of convection versus the static force balance is to compare the characteristic length scales of the contact angle-dependent Bond and Grashof numbers. As
hLBo =LGr ¼
1=2 1=3 h lv bl DT 1=3 g 1=6 ð l v Þ1=2 2=3
c q
q
m
ð22Þ
approaches unity, then the effect of the velocity boundary layer is reduced and bubble departure should be dictated by the static force balance between surface tension and buoyancy forces. For the nine different data sets in Fig. 4, we show values of hLBo =LGr . It can be seen that the lowest value of hLBo =LGr is for Nam et al., which happens to be where our departure model is least successful. Thus, our model should be more accurate for situations where hLBo =LGr > 5 (approximate). Results of our growth and departure model calculated for contact angles up to 90° and DT ¼ 20 C are shown in Fig. 5 along with the boundary where hLBo =LGr ¼ 5 which delineates the limit of our model’s applicability. When hLBo =LGr < 5, we recommend using a static departure description such as the spherical cap departure criterion (Eq. (21)) of Fritz correlation (Eq. (2)). According to Fig. 5b,c, a static departure description is valid for very hydrophilic surfaces (h < 20 ), which is associated with small bubble departure sizes. As shown in Fig. 4, modeled values for Dequiv;depart agree with previous studies within 4:1 to 5:6% (2r). Modeled values for tdepart agree with previous studies within 16 to 28% (2r). With the ability to determine bubble departure and model growth, it is possible to calculate the heat transfer rate associated with a single nucleation site (Fig. 5c): 4
qsingle ¼ qv hfg 3
p Dequiv;depart =2 t depart þ twait
3 ð23Þ
where twait is the waiting time between bubble departure and bubble growth.
Fig. 4. Our bubble growth and departure model for saturated water at atmospheric conditions is in good agreement with experimentally-verified simulation results from the literature. Bubble growth is modeled according to Eq. (15) and departure is determined as described in Section 4. Comparisons are shown for (a) ðh ¼ 54 ; DT ¼ 10:0 CÞ [4], (b) ðh ¼ 38 ; DT ¼ 12:5 CÞ [4], (c) ðh ¼ 38 ; DT ¼ 6:2 CÞ [4], (d) ðh ¼ 50 ; DT ¼ 7:0 CÞ [4], (e) ðh ¼ 25 ; DT ¼ 8:5 CÞ [4], and (f) ðh ¼ 30 ; DT ¼ 8:5 CÞ [4], (g) ðh ¼ 38 ; DT ¼ 8:5 CÞ [4], (h) ðh ¼ 50 ; DT ¼ 8:5 CÞ [5], and (i) ðh ¼ 10 ; DT ¼ 5:3 CÞ [6].
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
7
Fig. 5. Bubble departure results from our model for saturated water at atmospheric pressure. (a) Departure equivalent diameter is shown as a function of superheat (colors) and contact angle. Departure equivalent diameter is determined from bubble size (Eq. (15)) at the moment of bubble departure. How the bubble departs depends on hLBo =LGr . If hLBo =LGr > 5, the momentum boundary layer is important and departure is as described in Section 4—a simple approximation of Dequiv;depart is given by Eq. (21). Conversely, if hLBo =LGr < 5, then departure is described by a static departure description (surface tension vs buoyancy) such as the spherical cap departure criterion (Eq. (21)) or the Fritz correlation (Eq. (2)). The static descriptions of departure size have no dependence on superheat (colors). For both momentum and static regimes (spherical cap criterion), (b) departure time and (c) single bubble heat transfer are plotted as a function of superheat and contact angle. Single bubble heat transfer takes into account latent heat, departure size, departure time, and waiting time (Eqs. (23) and (24)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
To calculate the waiting time accurately, one should take into account superheat reduction effects during the bubble growth stage, which could include microlayer evaporation effects. For instance, microlayer evaporation could lower superheat within the growth period, effectively increasing the waiting time for a bubble. This effect could also be highly dependent on the Jakob number. For simplicity, we have used an approximate twait assuming semi-infinite conduction into the thermal boundary layer [31,34] and a thermal boundary layer as dictated by the Rayleigh number:
t wait ¼
d2T
pal
¼
m2=3
: 1=3
ð24Þ
g 2=3 bl2=3 DT 2=3 pal
From Fig. 5c, it is apparent that lower contact angles have lower heat transfer. Thus, if nucleation density were constant for a given surface, higher contact angles would be favorable for better heat transfer. Using the results of our growth and departure models, the departure diameter has been correlated to h and DT as
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dequiv;depart ¼ 0:652h0:581 DT 0:313
clv
g ðql qv Þ
ð25Þ
with R2 ¼ 0:999 for both our model results as well as simulation results [4,5].
It should be noted that our departure model assumes natural convection boundary layers, which is completely reasonable for an isolated bubble. However, at high nucleation densities, bubbles will start to interact with each other and the natural convection assumption will break down. This could be addressed using a different boundary layer correlation such as the one suggested by Tien [35]:
dT ¼ 3:22
kl hboil
ð26Þ
where hboil is the heat transfer coefficient in boiling and kl is the liquid conductivity. It should also be noted that we have assumed that wetting properties are constant. However, our departure model, like the growth model, can also be extended to cases where wetting properties are time-dependent (Section S4, Fig. S6). 5. Conclusions Studying the bubble cycle of nucleation, growth, and departure is fundamental and necessary to gain a more complete understanding of nucleate pool boiling. We have found that nucleation is linked to bubble departure, which suggests that nucleation criteria is dependent on wettability and convection. We have developed a growth model that assumes a bubble grows within a thermal boundary layer as opposed to the uniform superheating assumed
8
H.J. Cho, E.N. Wang / International Journal of Heat and Mass Transfer 145 (2019) 118803
in previous theories. We theorized a new mechanism of bubble departure wherein a velocity boundary layer is responsible for the rewetting underneath the bubble. Both the growth and departure models agree well with experimentally verified numerical simulation studies of single bubbles. By tracking the movement of the bubble base, we have also developed a methodology to calculate growth and departure for contact angle changing with time, which may be due to changes in surface tension and/or contact angle. Through the use of a simplified model, the results of this work provide a more physical understanding of the boiling process and simple quantitative tools to model boiling heat transfer. To improve upon this understanding further, future work could incorporate more detailed modeling of dynamic contact angle and microlayer evaporation, as well as explore non-aqueous fluids, mixtures, and nanofluids. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This work has been funded by the Singapore-MIT Alliance for Research and Technology (SMART). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer. 2019.118803. References [1] R. Cole, Bubble frequencies and departure volumes at subatmospheric pressures, AIChE J. 13 (4) (1967) 779–783, https://doi.org/10.1002/ aic.690130434. http://doi.wiley.com/10.1002/aic.690130434 . [2] B. Mikic, W. Rohsenow, P. Griffith, On bubble growth rates, Int. J. Heat Mass Transf. 13 (4) (1970) 657–666, https://doi.org/10.1016/0017-9310(70)900402. http://linkinghub.elsevier.com/retrieve/pii/0017931070900402 . [3] S. van Stralen, M. Sohal, R. Cole, W. Sluyter, Bubble growth rates in pure and binary systems: combined effect of relaxation and evaporation microlayers, Int. J. Heat Mass Transf. 18 (3) (1975) 453–467, https://doi.org/10.1016/0017-9310 (75)90033-2. http://linkinghub.elsevier.com/retrieve/pii/0017931075900332 . [4] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface, J. Heat Transf. 121 (3) (1999) 623, https://doi.org/10.1115/1.2826025. http://heattransfer. asmedigitalcollection.asme.org/article.aspx?articleid=1443768 . [5] A. Mukherjee, S.G. Kandlikar, Numerical study of single bubbles with dynamic contact angle during nucleate pool boiling, Int. J. Heat Mass Transf. 50 (1–2) (2007) 127–138, https://doi.org/10.1016/j.ijheatmasstransfer.2006.06.037. http://linkinghub.elsevier.com/retrieve/pii/S0017931006004157 . [6] Y. Nam, E. Aktinol, V.K. Dhir, Y.S. Ju, Single bubble dynamics on a superhydrophilic surface with artificial nucleation sites, Int. J. Heat Mass Transf. 54 (7–8) (2011) 1572–1577, https://doi.org/10.1016/j.ijheatmasstransfer.2010.11.031. http:// linkinghub.elsevier.com/retrieve/pii/S0017931010006459 . [7] J. Kim, Review of nucleate pool boiling bubble heat transfer mechanisms, Int. J. Multiph. Flow 35 (12) (2009) 1067–1076, https://doi.org/10.1016/j. ijmultiphaseflow.2009.07.008. http://linkinghub.elsevier.com/retrieve/pii/ S0301932209001311 . [8] N.I. Kolev, How Accurately Can We Predict Nucleate Boiling?, Springer, Berlin, Heidelberg, 2012, pp 143–178, https://doi.org/10.1007/978-3-642-21372-4_6. [9] H.-J.J. Cho, D.J. Preston, Y. Zhu, E.N. Wang, Nanoengineered materials for liquid–vapour phase-change heat transfer, Nat. Rev. Mater. 2 (2) (2016) 16092, https://doi.org/10.1038/natrevmats.2016.92. http://www.nature.com/ articles/ natrevmats201692 . [10] S.G. Bankoff, Entrapment of gas in the spreading of a liquid over a rough surface, AIChE J. 4 (1) (1958) 24–26, https://doi.org/10.1002/aic.690040105. http://doi.wiley.com/10.1002/aic.690040105 . [11] P. Griffith, J.D. Wallis, The Role of Surface Conditions in Nucleate Boiling, Tech. rep., Massachusetts Institute of Technology, Cambridge, MA, 1958.
[12] J.J. Lorenz, The Effects of Surface Conditions on Boiling Characteristics (Ph.D. thesis), Massachusetts Institute of Technology, 1972. [13] Y.Y. Hsu, On the size range of active nucleation cavities on a heating surface, J. Heat Transf. 84 (3) (1962) 207, https://doi.org/10.1115/1.3684339. http:// heattransfer.asmedigitalcollection.asme.org/article.aspx?articleid=1432587 . [14] M.S. Plesset, S.A. Zwick, The growth of vapor bubbles in superheated liquids, J. Appl. Phys. 25 (4) (1954) 493–500, https://doi.org/10.1063/1.1721668. http:// aip.scitation.org/doi/10.1063/1.1721668 . [15] S.J.D. van Stralen, R. Cole, W.M. Sluyter, M.S. Sohal, Bubble growth rates in nucleate boiling of water at subatmospheric pressures, Int. J. Heat Mass Transf. 18 (5) (1975) 655–669, https://doi.org/10.1016/0017-9310(75)90277-X. [16] V.P. Carey, Liquid Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation Processes in Heat Transfer Equipment, second ed., Taylor & Francis, 2007. [17] P. Stephan, T. Fuchs, E. Wagner, N. Schweizer, Transient local heat fluxes during the entire vapor bubble life time, in: ECI International Conference on Boiling Heat Transfer, Brazil, 2009. [18] T. Yabuki, O. Nakabeppu, Heat transfer mechanisms in isolated bubble boiling of water observed with MEMS sensor, Int. J. Heat Mass Transf. 76 (C) (2014) 286–297. [19] W. Fritz, Maximum volume of vapor bubbles, Phys. Z. 36 (11) (1935) 379–384. [20] L. Cheng, D. Mewes, A. Luke, Boiling phenomena with surfactants and polymeric additives: a state-of-the-art review, Int. J. Heat Mass Transf. 50 (13–14) (2007) 2744–2771. [21] A. Zou, A. Chanana, A. Agrawal, P.C. Wayner, S.C. Maroo, Steady state vapor bubble in pool boiling, Sci. Rep. 6 (1) (2016) 20240, https://doi.org/10.1038/ srep20240. http://www.nature.com/articles/srep20240 . [22] J.N. Addoms, Heat Transfer at High Rates to Water Boiling Outside Cylinders (Ph.D. thesis), Massachusetts Institute of Technology, 1948. [23] A. Sloan, S. Penley, R.A. Wirtz, Sub-atmospheric pressure pool boiling of water on a screen-laminate enhanced surface, 2009 25th Annual IEEE Semiconductor Thermal Measurement and Management Symposium, vol. 50, IEEE, 2009, pp. 246–253, https://doi.org/10.1109/STHERM.2009.4810771. http://ieeexplore. ieee.org/document/4810771/ . [24] S.M. Kwark, M. Amaya, R. Kumar, G. Moreno, S.M. You, Effects of pressure, orientation, and heater size on pool boiling of water with nanocoated heaters, Int. J. Heat Mass Transf. 53 (23–24) (2010) 5199–5208, https://doi.org/ 10.1016/j.ijheatmasstransfer.2010.07.040. http://linkinghub.elsevier.com/ retrieve/pii/S0017931010004187 . [25] X.-F. Yang, Z.-H. Liu, Pool boiling heat transfer of functionalized nanofluid under sub-atmospheric pressures, Int. J. Therm. Sci. 50 (12) (2011) 2402–2412, https://doi.org/10.1016/j.ijthermalsci.2011.07.009. http://linkinghub. elsevier.com/retrieve/pii/S1290072911002213 . [26] W. Tong, A. Bar-Cohen, T. Simon, S. You, Contact angle effects on boiling incipience of highly-wetting liquids, Int. J. Heat Mass Transf. 33 (1) (1990) 91– 103, https://doi.org/10.1016/0017-9310(90)90144-J. http://www.sciencedirect. com/science/article/pii/001793109090144J . [27] W. Wagner, A. Pruß, The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use, J. Phys. Chem. Ref. Data 31 (2) (2002) 387–535, https://doi.org/10.1063/1.1461829. http://aip.scitation.org/doi/10.1063/1.1461829 . [28] C.H. Wang, V.K. Dhir, Effect of surface wettability on active nucleation site density during pool boiling of water on a vertical surface, J. Heat Transf. 115 (1993) 659, https://doi.org/10.1115/1.2910737. [29] N. Siemons, H. Bruining, H. Castelijns, K.-H. Wolf, Pressure dependence of the contact angle in a CO2–H2O–coal system, J. Colloid Interface Sci. 297 (2) (2006) 755–761. [30] J.B. Goodman, N.W. Krase, Solubility of nitrogen in water at high pressures and temperatures, Ind. Eng. Chem. 23 (4) (1931) 401–404. [31] H. Chi-Yeh, P. Griffith, The mechanism of heat transfer in nucleate pool boiling— Part I, Int. J. Heat Mass Transf. 8 (6) (1965) 887–904, https://doi.org/10.1016/00179310(65)90073-6. http://linkinghub.elsevier.com/retrieve/pii/0017931065900748 http://linkinghub.elsevier.com/retrieve/pii/0017931065900736 . [32] R. Cole, W.M. Rohsenow, Correlation of bubble departure diameters for boiling of saturated liquids, Chem. Eng. Prog. Symp. Ser. (1969). [33] R. Raj, C. Kunkelmann, P. Stephan, J. Plawsky, J. Kim, Contact line behavior for a highly wetting fluid under superheated conditions, Int. J. Heat Mass Transf. 55 (9–10) (2012) 2664–2675. [34] M.Z. Podowski, R.M. Podowski, Mechanistic Multidimensional modeling of forced convection boiling heat transfer, Sci. Technol. Nucl. Instal. 2009 (2009) 1–10. [35] C.L. Tien, A hydrodynamic model for nucleate pool boiling, Int. J. Heat Mass Transf. 5 (6) (1962) 533–540. [36] H. Fedders, Messung des Wärmeüberganges beim Blasensieden von Wasser an metallischen Rohren, no. Jül-740 RB, Kernforschungsanlage Jülich, 1971. [37] J.V. Sengers, J.M.H.L. Sengers, Thermodynamic behavior of fluids near the critical point, Annu. Rev. Phys. Chem. 37 (1) (1986) 189–222, https://doi.org/ 10.1146/annurev.pc.37.100186.001201.