On the Intrinsic Properties Effects of the Heater Surface in Nucleate Pool Boiling: A Simulation Study

1. Introduction

In contemporary research, power dissipation has emerged as a critical factor across various domains, spanning from the production of miniaturized electronics to nuclear power plants, as highlighted in reviews by [1,2,3]. Within the electronic industry, as MEMS processing technology advances towards high-speed, large buffer memory in smaller devices, there is a growing focus on the extensive study of pool boiling to achieve efficient cooling.

Nucleate boiling proves to be effective in dissipating a substantial heat flux with a minimal temperature difference, facilitating efficient component cooling. In pool boiling, heat is transferred from a surface to liquid in a macroscopic state of rest. The initiation of boiling occurs when the local temperature is sufficiently high to permit the formation and growth of vapor bubbles on surface imperfections. These imperfections or cavities act as nucleation sites for the growth of bubbles, arising from vapor initially trapped within them [4]. Surface roughness, porosity, and wettability are believed to influence boiling behavior, although bubble nucleation can also take place on a smooth surface without imperfections. In such cases, nucleation energy is contingent upon the contact angle, defined as the angle between a tangent to the liquid surface and the solid surface. Typically, a right angle is considered neutrally wetting, with lesser angles indicating hydrophilic properties and greater angles indicating hydrophobic characteristics. Both hydrophilic and hydrophobic surfaces contribute positively to improving boiling performance, with hydrophilic surfaces aiding in rewetting following bubble departure, and hydrophobic surfaces promoting bubble nucleation seeded by a stored vapor in cavities [1,2]. Additionally, the contact angle between the liquid-vapor interface at the bubble base and the surface varies during the stages of bubble growth and departure [5,6,7].

Recent research in nucleate pool boiling has emphasized the development of novel surface manipulation methods to enhance phase-change heat transfer. These methods include mechanical machining, chemical treatments, nanoparticle coatings, and micro-/nanoelectromechanical systems techniques, such as photolithography and reactive ion etching, as well as fiber-laser texturing that avoids additional layers [8,9,10,11,12]. In the context of ongoing efforts to improve heat transfer, this study investigates the impact of thermophysical properties on nucleate boiling for two materials with distinct characteristics: copper and silicon oxide. Notably, copper exhibits thermal diffusivity and conductivity two orders of magnitude greater than silicon oxide.

While numerous numerical models and experiments have explored the bubble cycle and precise bubble cycle data [13,14,15], predicting boiling heat transfer remains intricate due to the necessity of considering phenomena occurring over multiple scales—from the adsorbed vapor at the nanometer scale [16,17] to the bubble diameter at the millimeter scale. Timescales also vary widely, from microlayer formation at the microsecond scale ($10-100\mu s$) to evaporation at the millisecond scale ($1-10ms$), significantly contributing to bubble growth [18,19,20,21,22]. Our proposed model is highly versatile, spanning multiple scales from nano to milli, and timescales from tenths of milliseconds to minutes, allowing the independent variation of parameters through computational simulations.

Reproducing boiling data is challenging due to contamination and cavity reactivation. Thus, our model analyzes bubble departure frequency in small copper and silicon oxide heaters, maintaining a fixed contact angle and comparing the results to established expressions and experimental data. Focusing on bubble departure frequency, which is simpler to observe during experimentation than other bubble characteristics, serves as a valuable metric for assessing the effectiveness of boiling on a given surface. This frequency is linked to the heater’s dissipation potency, depending on both the growth time of the bubbles and the recovery time of the microlayer of fluid, which rewets and cools the heater upon bubble departure [23].

Over the past half-century, significant efforts have been directed towards developing heat transfer models to determine these times. A theory presented in [24] outlines the rate of growth of a vapour bubble from a heated wall in a liquid near its saturation temperature. Predictions from this analysis bear some relation ( 25 per cent) to experimental results under conditions where assumptions are approximately valid. A corrected version is established in [25], substantially reducing the growth constant of bubbles and approaching experimental results more closely. However, recovery or wait times notably vary in each experiment. An expression for wait time was analytically derived by C.-H. Han and P. Griffith [26], considering transient heat conduction through the heater at a constant temperature to the microlayer of fluid, but significantly underestimate to experimental values of wait times. In an attempt to bridge this gap, Podowski et al. [27] consider the wall surface temperature as a time-dependent, fluctuating parameter, solved through a one-dimensional heat conduction equation through the wall. Although this consideration introduces a dependence on wall material properties in wait times, it still deviates from experimental results. In [28], an empirical relationship is proposed, fitting well with experimental data, and is the one utilized in our model.

Through numerical simulations, this study systematically investigates the effects of the thermophysical properties of the heater on nucleate pool boiling for two distinctly different heater materials, testing the proposed different mechanisms with semi-empirical relationships known in the literature.

2. Model

Consider a scenario of pool boiling, where a small heater plate is positioned horizontally in a container filled with stagnant refrigerant at atmospheric pressure (see Figure 1). The model consists of two distinct domains operating at different scales: the heat conduction process within the heater material and the growth and departure of bubbles from the surface. These domains are interconnected through the heat transfer mechanisms involved in surface cooling.

2.1. Heat conduction in the heater

Starting with the Fourier equation given by

$$\frac{1}{{\alpha }_m}\frac{dT}{d\tau }={\overline{\nabla }}^{2}T+\frac{\dot{q}}{{\kappa }_m},$$

(1)

the temperature field in the heater is estimated using the alternating-direction implicit (ADI) finite difference method, which is modified by an f factor ($0<f<1$) [30]. In Eq. (1), ${\alpha }_m$ and ${\kappa }_m$ represent the heater’s thermal diffusivity and conductivity, respectively, and $\dot{q}$ denotes the heat flux sources and/or sinks.

The heater is divided into grid cells, and in each grid cell, the heat fluxes due to the first term in Eq. (1) are affected by the factor f, remaining the total heat flux according to Eq. (1) in each direction over a complete time step. The advantage of the modified ADI method lies in the fact that not only tridiagonal matrices need to be solved, but also the time steps can be significantly larger without compromising the convergence and stability of the solution. The time step limit for the conventional ADI method can now be increased by a factor of $1/f$ through the use of this modified ADI method. For example, the time step limit can be increased by one order of magnitude with $f=0.1$, and the solutions still remain stable with high accuracy (see [30] for more details).

The heat generated by the heater is extracted from its top into the liquid. The extraction mechanisms include radiation, microconvection, microlayer evaporation due to bubble growth, and natural convection generated by the fraction of the free cell surface in contact with the liquid without bubbles. Radiation is always present, whereas the other extraction mechanisms are mutually exclusive. However, radiation heat is negligible in nucleate boiling. In essence, microconvection and microlayer evaporation are the relevant processes for heat extraction, despite being intermittent (refer to Section 2.2). As these mechanisms undergo relatively rapid changes over time, solving Eq. (1) using the ADI method is advantageous, especially since the precision is on the order of 2 in time ($O\left(t^2\right)$), allowing for better tracking of temporal changes compared to explicit finite difference methods.

2.2. Heat transfer mechanisms

The mechanisms for extracting thermal energy from the heater and transferring it to the liquid are illustrated in Figure 2:

• Natural convection ($q_{nc}$): The density difference between the hot and cold fluid induces convective flows, extracting heat from the heater across its free area, $A_{n}c$.
• Microlayer evaporation, ($q_{me}$): It is assumed that the vapor mass forming the bubble originates from the evaporation of the superheated microlayer between the heater surface and the bubble itself.
• Microconvection, ($q_{mc}$): Detaching bubbles remove the superheated liquid layer in their vicinity. Colder liquid from the bulk of the pool quenches an area $A_{mc}$, and heat is required to recover the microlayer subcooling.
• Radiation, ($q_{rad}$): Heat is transferred in the form of electromagnetic waves from the heater to the liquid.

The total heat rate extracted from a cell to the liquid, ${\dot{q}}_{out}$, is given by the sum of individual contributions:

$${\dot{q}}_{out}={\dot{q}}_{nc}+{\dot{q}}_{mc}+{\dot{q}}_{me}+{\dot{q}}_{rad}.$$

(2)

The net heat rate extracted from a cell, $\dot{q}$ (see Eq. (1)), represents the difference between the power generated in the cell and the heat rate extracted from the upper cells through the aforementioned mechanisms.

Nucleation sites (cavities) are randomly placed at the cell centers of the heater, with radii chosen from normal distributions given by $6\pm 3;\mu m$ [32]. At the active sites, bubbles grow attached to the cells until a critical detachment radius, $r_d$, at which point bubbles depart, initiating the growth process again with a new bubble if the site becomes reactivated.

If a growing bubble touches other bubbles on the surface, they coalesce to form a single bubble at the nucleation site of the largest one, preserving the vapor mass.

2.2.1. Effect of the contact angle and natural convection

The current model incorporates the influence of the contact angle $\varphi$ on the site activation process through the activation rule. Accordingly, a given site is active if its radius r satisfies [25]:

$$r_{min}<r<r_{max},$$

(3)

where

$$r_{min,max}=\frac{\delta }{2C_1}\left\{1-\frac{{\theta }_s}{{\theta }_w}\pm \sqrt{{\left[1-\frac{{\theta }_s}{{\theta }_w}\right]}^2-\frac{4\zeta C_3}{\delta {\theta }_w}}\right\}$$

(4)

with

$$\zeta =\frac{2\sigma T_{sat}}{{\rho }_v{h}_{fg}},C_1=\frac{1+cos\varphi }{sin\varphi }and{C}_3=1+cos\varphi ,$$

(5)

where

$\delta$: thermal boundary layer thickness

$\varphi$: contact angle between a tangent to the liquid surface and the solid surface

${\theta }_s$: liquid subcooling, i.e., $T_{sat}-T_{\infty }$

${\theta }_w$: wall superheat of the cell containing the nucleation site, i.e., $T_w-T_{sat}$

${\rho }_v$: vapor density

$h_{fg}$: latent heat of phase change

$T_{sat}$: saturation temperature

$T_{\infty }$: liquid bulk temperature

$\sigma$: surface tension.

The thermal boundary layer thickness is estimated as $\delta =\kappa /\overline{h}$, where $\kappa$ is the liquid heat conduction coefficient, and $\overline{h}$ is the average heat transfer coefficient assessed using the average temperature of the heater. The average heat transfer coefficient is given by [33,34,35]:

$$\overline{h}=\frac{\kappa }{L}\overline{Nu},$$

(6)

where L is the heater side characteristic length, $L=4\frac{A_L}{P}$, with $A_L$ and P being the heater area and perimeter, respectively. $\overline{Nu}$ is the Nusselt number calculated as:

$$\overline{Nu}=\left\{\begin{array}{cc}\hfill 0.54c_{nc}R{a}^{1/4},& \mathrm{if}{10}^4\le R{a}_L<{10}^7\\ \hfill 0.15c_{nc}R{a}^{1/3},& \mathrm{if}{10}^7\le R{a}_L<{10}^{11}\end{array}\right.$$

(7)

where the Rayleigh number is defined as:

$$Ra=\frac{g{\beta }_l\left(T_w-T_{\infty }\right)L^3}{\nu {\alpha }_l},$$

(8)

with g being the gravitational constant, $\beta$ the thermal expansion coefficient of the liquid, $\nu$ the kinematic viscosity, and ${\alpha }_l$ the thermal diffusivity of the liquid. The coefficient $c_{nc}$ is included to account for border effects in small heaters and extra agitation due to the detachment of nearby bubbles.

Finally, natural convection heat transfer acts on the portion of the heater surface in direct contact with the liquid, as shown in Figure 2. For every top cell of the heater, the convection heat transfer rate is given by:

$$\frac{d{q}_{nc}}{dt}=c_{nc}\overline{h}{A}_i\left(T_w-T_{\infty }\right),$$

(9)

where $A_i$ is the area of cell i in contact with the liquid.

2.2.2. Microlayer evaporation and microconvection

In the microlayer evaporation model initially proposed by N. W. Snyder [36], a thin liquid microlayer is present beneath the hemispherical bubble during a specific stage of growth, occurring shortly after the initial formation of the equilibrium radius. Heat is transmitted through the microlayer, and its thickness diminishes due to evaporation in accordance with

$$q_{me}=\frac{4}{3}\pi \left(r{(t+\Delta t)}^3-r{\left(t\right)}^3\right){\rho }_v{h}_{lg},$$

(10)

where $r\left(t\right)$ is radius of the bubble at time t, and $\Delta t$ is the Monte Carlo step. The asymptotic growth of the radius of the attached bubble is provided by [25]:

$$r\left(t\right)=0.470P{r}^{-1/6}Ja{\left({\alpha }_{l}t\right)}^{1/2},$$

(11)

where, $Pr$ is the Prandtl number, $Ja$ is the Jakob number. Eq. (11) is valid until the growing bubble radius reaches the detaching critical value $r_d$, given by the Fritz criterion:

$$r_d=\frac{1}{2}0.0148\varphi \sqrt{\frac{2\sigma }{g\left({\rho }_l-{\rho }_v\right)}},$$

(12)

where g is the gravity and ${\rho }_l$ and ${\rho }_v$ are the liquid and vapor densities, respectively. When a bubble departs from the hot surface, it carries away a fraction of the superheated thermal layer. Colder liquid from the bulk of the pool quenches the surface. If the site on the surface is active, it implies that a new superheated thermal layer was formed (see Figure 2). According to this, the heat removal by the microconvection mechanism can be calculated using [37]

$$q_{mc}=c_l{\rho }_l\frac{2}{3}\pi r_d{\left(c_{mc}{r}_d\right)}^2\left(\frac{T_w-T_{\infty }}{2}-T_{\infty }\right),$$

(13)

where $c_l$ is the specific heat of the liquid, and $c_{mc}$ is the microconvection coefficient representing the liquid volume around a bubble that is carried with the bubble after detachment. In the present model, it is taken to be equal to $1.2-1.3$, as suggested by [37], page 404.

The time required from the beginning of the generation of the thermal layer to the beginning of bubble growth is defined as the waiting period of a bubble, $t_w$. Basu et al. [28] correlated the wait time and found that it was a strong function of wall superheat according to

$$t_w=\frac{139.1}{{(T_{cavity}-T_{sat})}^{4.1}},$$

(14)

where $T_{cavity}$ is the temperature of the cavity when the bubble departs. Thus, the heat $q_{mc}$ is removed by microconvection during the wait time $t_w$. During that time of recovery of the microlayer, non-bubble growth occurs, and the cavity remains inactive. (The microlayer affected occupies an area double that of the bubble).

2.2.3. Radiation

The heat per unit of time removed through the radiation process is accounted for in the model by

$$\frac{d{q}_{rad}}{dt}=A_Tϵ{\sigma }_{sb}\left(T_w^4-T_{\infty }^4\right),$$

(15)

where $A_T$ is the total area of the cell, $ϵ$ is the emissivity of the upper surface of the heater (assumed to be equal to 0.8 in the present model), and ${\sigma }_{sb}$ is the Stefan-Boltzmann constant (${\sigma }_{sb}=5.67\times {10}^{-8};W/m^2{K}^4$)

3. Results

3.1. Comparative Cooling Performance of $Si{O}_2$ and $Cu$ Heaters

In this section, we examine how the intrinsic properties of the heater affect its refrigeration. The transfer of heat flux from the heater to the liquid relies on both the mean frequency of bubbles, denoted as $\nu$, and their number, denoted as $n^{\prime }$, on the heater surface. In the literature [38], it is noted that this transfer is proportional to $\nu$ for the mean number of active sites. Figure 3 illustrates the product $\nu \times n^{\prime }$ in the steady state as a function of the input power per area $q"$ for both copper ($Cu$) and silicon dioxide ($Si{O}_2$) heaters using water as the refrigerant, with a fixed contact angle of $\varphi ={20}^{\circ }$. The heaters share identical dimensions of $1\times 1\times 0.25c{m}^3$, with 20 cavities randomly distributed, sized according to Gaussian distributions given by $6\pm 3;\mu m$. As $q"$ increases, $\nu \times n^{\prime }$ rises for both heater types. Notably, for a given input power per area, the refrigeration, i.e., $\nu \times n^{\prime }$, is higher for copper heaters compared to silicon dioxide heaters. Thus, materials with greater thermal conductivity and diffusivity achieve better refrigeration through nucleate boiling.

It is important to note that these results are applicable under the assumption of a fixed contact angle and a constant cavity densities with Gaussian size distributions. In practical experiments, maintaining constant surface properties can be challenging. In essence, the model facilitates the systematic study of this phenomenon by varying one parameter at a time.

Upon cooling less efficiently, the $Si{O}_2$ heater attains a higher mean temperature than the $Cu$ heater, as illustrated in Figure 4 Temp-Bubbles (a-d) and (e-h) for $Si{O}_2$ and $Cu$ heaters, respectively, under a constant input power per area $q"$ at different time intervals. Furthermore, upon comparing these figures, it becomes evident that temperature fluctuations are more pronounced in the $Si{O}_2$ heater than in the $Cu$ heater. These variations are particularly noticeable in the active cavities, where bubbles grow and subsequently depart from the surface.

3.2. Mean Bubble Frequency and Associated Times: A Comparison between $Si{O}_2$ and $Cu$ Heaters

In Figure 5(a), the mean bubble frequency $\nu$ in the steady state is depicted against the input power per area $q"$ for both $Cu$ and $Si{O}_2$ heaters, with a fixed contact angle of $\varphi ={20}^{\circ }$ and 20 cavities randomly distributed according to Gaussian distributions ($6\pm 3;\mu m$). The mean frequency is calculated by averaging over time windows when the steady state is reached. Transient states, lasting up to ten seconds depending on the heater material, the refrigerant liquid, and $q"$, are illustrated in Figure 5(b). In this figure, the power per area extracted $q{"}_{mc}$ by microconvection is plotted for randomly chosen active cavities where the microlayer is recovering. This mechanism extracts significant heat volumes, especially for the $Cu$ heater, reaching around 100 $W/c{m}^2$ when the microlayer is recovering, as reported in experiments [39]. Notably, Figure 5(b) reveals substantial variations in $q{"}_{mc}$ for the $Si{O}_2$ heater, consistent with the temperature variations observed in the active cavities (refer to Figure 4 (e-h)).

Figure 5(c) displays the bubble radius growth $r\left(t\right)$ on an active site as a function of time for both $Cu$ and $Si{O}_2$ heaters under fixed input power per area. The radii asymptotically grow until reaching the critical radius $r_D$ at which point the bubble departs from the heater surface (see Eqs. (11) and (12)). During bubble growth, the cavity temperature tends to decrease, as depicted in Figure 5(d), where the asymptotic T of the cavity versus time is plotted for the same data as in part (c) of the same figure. For the $Si{O}_2$ heater, the variations in T of the cavity are greater than those for the $Cu$ heater. Moreover, quiescent states are observed in bubbles on the $Si{O}_2$ heater, where the growth of their radii halts due to a collective effect with neighboring cavities. This collective behavior leads to a decrease in T of the cavity following the eventual growth and departure of bubbles (refer to parts (c) and (d)). Preliminary results from spatial correlations support this assertion, although they are beyond the scope of this study.

Finally, for saturated pool boiling, a recommended empirical correlation for the bubble frequency times the departure diameter $D_d=2r_d$ is the following [40]:

$$\nu D_d=1.18\frac{t_g}{t_g+t_d}{\left[\frac{\sigma g({\rho }_l-{\rho }_g)}{{\rho }_l^2}\right]}^{1/4}.$$

(16)

Here, $t_g$ represents the average growth time of a bubble, and $t_d$ is the average time lag between consecutive bubbles (Figure 5(c)). It is important to note that this latter time can be greater than the wait time $t_w$ defined by Eq. (14), as the microlayer may recover, and the cavity is influenced by the temperature of neighboring cavities. The values of $t_g$ and $t_d$ can be calculated using the evolution of the bubble radius at each active nucleation site as predicted by the model. The solid curves in Figure 6 are fitted to Eq. (16) for $Cu$ (violet line) and $Si{O}_2$ (orange line) heaters, and they are comparable to $\nu \times D_d$ for each type of heater (violet circle and orange circle for $Cu$ and $Si{O}_2$ heater, repctively).

4. Conclusions

We analyze the effect of the intrinsic properties of the heater material on nucleate boiling through a model introduced in this work. This model incorporates many ingredients known in the literature, describing the various mechanisms involved in this phenomenon. In experiments, it is impossible to change one parameter without modifying the others. In other words, these experiments are challenging to reproduce even under the same conditions. For example, the vapor trapped in surface imperfections changes significantly from one trial to another, making it difficult to study which mechanism is predominant when comparing results. Simulations allow us to overcome this problem and study changes by modifying one parameter at a time.

We compared two materials, $Cu$ and $Si{O}_2$, and found a substantial difference in nucleate boiling refrigeration between a good thermal conductor like $Cu$ and a poor conductor like $Si{O}_2$. Even assuming that the properties of both surfaces remain the same, i.e., they have the same density of cavities that can be activated with a certain Gaussian size distribution and that both surfaces maintain the contact angle with the refrigerant liquid, we found that the average bubble frequency varies from one material to another. This is attributed to both the bubble growth time and the time lapse between bubbles generated in the same active cavity, both of which exhibit a notable dependence on the material. We observed an influence of neighboring bubbling sites, especially in $Si{O}_2$ due to its very low thermal conductivity and diffusivity. In fact, temperature variations become highly visible in such materials due to their significant inertia in recovery. Finally, we contrasted the simulation results with a semi-empirical correlation known for the bubble frequency times the departure diameter, obtaining good agreement between them (see Eq. (16) and Figure 6).