A numerical investigation of the effect of surface wettability on the boiling curve

Surface wettability is recognized as playing an important role in pool boiling and the corresponding heat transfer curve. In this work, a systematic study of pool boiling heat transfer on smooth surfaces of varying wettability (contact angle range of 5° − 180°) has been conducted and reported. Based on numerical simulations, boiling curves are calculated and boiling dynamics in each regime are studied using a volume-of-fluid method with contact angle model. The calculated trends in critical heat flux and Leidenfrost point as functions of surface wettability are obtained and compared with prior experimental and theoretical predictions, giving good agreement. For the first time, the effect of contact angle on the complete boiling curve is shown. It is demonstrated that the simulation methodology can be used for studying pool boiling and related dynamics and providing more physical insights.


Introduction
Boiling occurs in a variety of industrial applications such as high heat flux electronic devices [1], chemical processes [2], power plants [3], etc. During boiling, a heated surface is adjacent to a liquid, which vaporizes. The large latent heat of vaporization makes it an efficient mode of heat transfer in the nucleate boiling mode. Boiling is quantified in terms of a plot of heat flux versus the wall superheat defined as the temperature of the wall minus the saturation temperature (i.e. the boiling point) at the pressure of the liquid. The heat flux divided by the wall superheat is the heat transfer coefficient (HTC). Large HTC is an indication of efficient heat transfer.
There are various modes of boiling. During nucleate boiling, vapor bubbles form at a superheated surface which rise up in the liquid. Eventually, the heat flux reaches a maximum, which is called the critical heat flux (CHF). Increasing the superheat of the heated surface beyond the CHF value leads to drastically reduced HTC. This is because more vapor with lower thermal conductivity accumulates near the surface, eventually forming a stable film in the film boiling mode. In pool boiling applications, such as electronic equipment cooling, the drastic reduction in boiling heat flux after CHF may lead to devastating results.
Significant enhancement of HTC during boiling has been reported for porous surfaces [4], and surfaces with carbon nanowires [5,6], silicon [7], and copper [7,8]. In all cases, it was a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 found that there is significantly more vapor bubble nucleation in the presence of nanoscale roughness compared to smooth surfaces. Another factor that affects boiling heat transfer is wettability, which is the ability of a liquid to wet (or be spread over) a solid surface [9,10]. It is known that wettability can change the temperature at which CHF occurs or the temperature at which the transition of film boiling to nucleate boiling occurs, i.e. Leidenfrost point. Surface characteristics are critical to determine the efficiency of heat transfer. Surface properties such as roughness and wettability can alter the transition temperature from the nucleate boiling to the film-boiling regime. Experiment has revealed that hydrophobic surface brings out high critical heat flux (CHF), and hydrophobic surface benefits bubble nucleation. Increasing the energy efficiency by tailoring optimal surfaces could profoundly impact many industrial applications. However, contact angle effect always comes with surface roughness effect and it is difficult to distinguish from one to the other.
Boiling process involves complicated physics of phase change. With the present technology and measuring equipment, it is still difficult to observe the real dynamics of interface phase change. However, experiment has captured the temperature distribution along the liquidvapor phase boundary [11,12]. As for numerical parts, there are only limited numbers of numerical studies on phase change since it is a multi-scale physics involved problem which still remains a great challenge using existing tools and the drastic density change between liquid and vapor within less than microscale interface causes numerical difficulties.
The objective of this study is to perform numerical simulations of boiling phenomena to investigate the effect of contact angle (CA) or wettability on boiling heat transfer using a volume of fluid (VOF) method combined with a static contact angle model. The primary focus is on reproducing qualitative trends in the whole boiling curve using one single model and gaining more physical insights into the underlying mechanisms. In this work, numerical simulations were performed for three-dimensional cylindrical and two-dimensional planar unsteady laminar flow of incompressible liquid and vapor. This was also the first time the effect of contact angle is shown in the complete boiling heat transfer curve.

Numerical method and modeling
In this work, the VOF method is chosen since the fluid mass can be conserved appropriately and it can be applied on a larger scale with any grid compared with the LBM and MD method. The continuing work can be extended to a multi-bubble problem to predict real industrial phase change applications.

Governing equations
In the VOF method, the two phases are represented by phase volume fractions such that where α is the volume fraction, subscripts l and v represent the liquid and vapor phases, respectively. The governing equations consist of the continuity equations for the two phases and a one-fluid model for the momentum and energy conservation equations.
@ @t ðrũÞ þ r Á ðrũũÞ ¼ À rp þ r Á ðmrũÞ þ rg þF; @ @t ðrC p TÞ þ r Á ðrũC p TÞ ¼ À r Á ðkrTÞ þ S h ; where ρ v and ρ l are the vapor and liquid densities, respectively, t is the time,ũ is the average fluid velocity, and _ m lv is the mass source due to liquid to vapor phase change. In the momentum conservation (third) equation, p is the pressure,g is the gravitational force, and ρ and μ are the density and viscosity of the mixture of liquid and vapor, respectively. In the energy conservation (fourth) equation, T is the temperature, C p is the mixture specific heat, k is the mixture thermal conductivity, and S h is the heat source. ρ, μ, C p , and k are: where subscripts l and v represent the liquid and vapor phases, respectively. Along the liquidvapor interface, surface tension results from the greater attraction force between liquid molecules than to the molecules in the vapor. Brackbill porposed a continuum surface force (CSF) model to include the surface tension effect [27], which is commonly used in the continuum VOF model. The origin of this source term can be considered from the specific case where the surface tension is constant along the interface, and where only the forces normal to the interface are considered. The pressure drops across the interface can be estimated in terms of the surface tension coefficient, σ, and the surface curvature as measured by two radii in orthogonal directions, R 1 and R 2 : Hence, in the momentum equation,F is the surface tension force between the two phases that is expressed as a volume force density: where σ is the surface tension and the interface curvature is given by: The effect of contact angle at fluid interface in contact with solid boundary then can be estimated within the CSF model in terms of θ w which is the equillibrium contact angle between the solid and fluid. It is a static contact angle which is measured when the fluid is at rest. If θ w is the contact angle at the wall, then the surface unit normaln at the calculation cell iŝ n ¼n w cos y w þn t sin y w ð7Þ wheren w is the unit vector normal to the wall andn t is the unit vector tangential to the wall.

Phase-change model
The mass and energy exchange between the two phases can be governed by some phase change model coupled with the governing equations [28]. In the model used in our simulations, the thermal conductivity of the interfacial grid cells is equal to that of the unsaturated phase [29]. The interfacial heat flux Q pc causing the liquid to vapor phase change is then calculated as [29] Q pc ¼ where C I 's are interfacial grid cells (see Fig 1), f denote faces of C I , A f is the face area, rT f is the temperature gradient at the cell face, andñ f ;in is the unit normal vector at the cell face which points into C I . The mass source due to phase change at C I can be calculated by: where V CI is the volume corresponding to grid cell C I , and h lv is the latent heat of vaporization. The corresponding heat source term in the energy equation is calculated as:

Problem definition
The primary goal is to qualitatively investigate how the contact angle influences the boiling curve, i.e., the heat transfer during various phases of boiling. Fig 2 shows the configuration used in this study. Computations were performed for three and two-dimensional, unsteady, incompressible flow with heat transfer. For the VOF method, an interface tracking method, small amount of vapor phase has to be specified at the initial of the computation [16,25]. This is the reason why the VOF method cannot be used to study the nucleation process. Vapor phase was initially placed at the bottom wall. If the film boiling mode is not stable, then the vapor wraps into a bubble and the dynamics proceeds in the nucleate or transition boiling regime. The governing equations were solved using the finite volume method. The governing equations are solved using the software ANSYS 1 Fluent Academic Research, Release 16.2.

Geometry and boundary conditions
Two dimensional simulation domain is with a width of λ 0 and a height of 3λ 0 , and three dimensional simulation uses cylinder computational domain. λ 0 is the Taylor-Rayleigh instability wavelength, which is calculated using the working fluid properties listed in Table 1: The boundary conditions are as follows. The upper boundary condition is set as a pressure outlet with temperature T = T sat . Symmetric boundary conditions are used for both sidewalls of the simulation domain. The bottom boundary is a no-slip superheated wall with a constant contact angle along the surface. The bottom wall has a constant temperature ranging from 2K to 100K above saturation temperature T sat = 373K. To consider the effect of contact angle, the bottom wall has a contact angle ranging from 5 o À 180 o . For simplicity, static contact angles were considered as the key trends were successfully resolved by static contact angle simulations although studying the effect of dynamic contact angles was feasible [25]. Initially, a linear temperature profile from T w at the bottom wall to T sat at the liquid-vapor interface is specified in the vapor domain. The liquid domain has an initial temperature equal to T sat . A "microlayer" is deemed relevant for hydrophilic cases [15,16]. In these cases, the microlayer contribution to heat transfer has been reported to be around 20% [23]. Hence following a prior work [25], and for simplicity, the microlayer was not modeled in this work [30].
The initial shape of the vapor-liquid interface is perturbed to initiate the bubble growth. The initial interface position, which is the height of vapor film from the bottom wall, is given by: Gravity is pointed in the vertically downward direction.

Fluid properties
For the purposes of obtaining better numerical stability, the density ratio of liquid and vapor cannot be large [31], so the working fluid is chosen to be an artificial one with properties listed in Table 1. For comparison, the properties of water (not the fluid in our simulations) are listed in Table 2. Further comparison between the properties of working fluid and water can be done based on non-dimensional parameters. The primary non-dimensional parameters are listed below: where [Á] denotes scale of the corresponding physical variable. In the above list, Ja is the Jakob number, Pe is the Peclet number, Fr is the Froude number, Re is the Reynolds number, We is the Weber number, and Pr is the Prandtl number. Liquid properties are used in defining these parameters. Additional non-dimensional parameters are the ratios of density, viscosity, thermal conductivity, and specific heat of the two phases (liquid and vapor). These ratios can be deduced from Tables 1 and 2. The non-dimensional parameters for the working fluid in our simulations are listed below: Pe ¼ 3:849; Re ¼ 3:85; We ¼ 1:023;  The Jakob number (Ja) represents the ratio of sensible heat to latent heat absorbed during the phase change process. The Jakob number of working fluid in simulations is larger than that of water which represents a better ability to phase change. The typical critical heat flux (CHF) for water in experiments is around 10 6 W/m 2 ; a lower CHF for the working fluid may be due to different Ja. The Reynolds numbers (Re) in both fluids are in the laminar regime. The Weber numbers (We), which depend on surface tension, are close to unity in both fluids. The Prandtl numbers (Pr) in both fluids are on the order of unity, implying that thermal diffusivity is close to momentum diffusivity.

Results and discussion
The numerical results are presented in three-dimensional cylindrical and two-dimensional planar simulations, respectively.

Two-dimensional planar simulation
The two-dimensional planar numerical results are presented in four subsections: transient bubble dynamics, boiling curve, critical heat flux (CHF), and Leidenfrost point (LFP).

Transient bubble dynamics.
Vapor is initially placed next to the bottom boundary which is superheated with a static contact angle imposed. If the film boiling mode is not stable, then the vapor wraps into a bubble and the dynamics proceeds in the nucleate or transition boiling regime. The continuing bubble generated depends upon the amount of remaining vapor accumulated at the heated wall. In the transition and film boiling regimes, a new vapor bubble is formed naturally from the interfacial instability. In the nucleate boiling regime, while the surface can become fully wetted without any vapor, one can numerically introduce repeated nucleation bubbles artificially. The rate of introduction of new nucleate would depend in nucleation models that take into account the nature of the surface. This gives no difference from prior models in literature. To validate this work, this simulation is compared with the result of Lattice Boltzmann method [16].   critical heat flux, transient boiling, and LFP, and contact angles of 10˚, 60˚, 120˚, and 160˚, respectively. Note that contact angles less than 90˚are regarded as hydrophilic, whereas those above 90˚are regarded as hydrophobic. The red fluid area represents the vapor phase state while the blue area represents the liquid phase state.
(i) Hydrophilic surface: For the contact angles of 5˚− 90˚, it is observed that liquid tends to wet the bottom wall more. For example, at contact angle = 10˚ (Fig 6(a)-1), the initial vapor layer forms a bubble and the liquid wets the wall. The vapor layer is increasingly unstable for more hydrophilic walls. In the nucleate boiling phase, bubbles are formed at the bottom wall. With an increase in the wall superheat, the vapor-production rate increases and the growth period decreases. It is seen that the bubble diameter at departure increases with wall superheat. For a fixed contact angle, the departure bubble diameter depends on the growth rate, which increases with wall superheat.
For a contact angle = 10˚, ΔT = 45K critical heat flux will result in Fig 6(a)-2. This is deduced based on a plot for heat flux (Fig 7). For wall temperatures around the critical heat flux condition, there are significant fluctuations in local surface heat flux due to bubble dynamics. The dry and wet regions change continuously. More vapor regions are formed during transient boiling (Fig 6(a)-3) with increasing wall superheat. In the case of contact angle = 10˚, the vapor film is stable at the bottom boundary beyond ΔT = 60K. Thus, in this case, ΔT = 60K is the Leidenfrost point (LFP). The LFP is also the minimum heat flux point (Fig 7). The vapor film covers the entire bottom boundary leading to the film boiling regime (Fig 6(a)-4) at superheats beyond the LFP.
(ii) Hydrophobic surface: For contact angles of 90 − 160˚, it is observed that original vapor layer is still unstable at low superheats but it is increasingly stable for greater contact angles. As a result, the nucleate and transition boiling phases occur much earlier. For a contact angle of 160˚, there is practically no transition boiling regime with a vapor film present in nearly all cases (Fig 6(d)-1 to 6(d)-4). The hydrophobic surface repels the liquid and stabilize the vapor film so that the small superheat causes Leidenfrost behavior. This is consistent with prior results [32].    Many models of CHF have been developed [30,[36][37][38][39][40][41]. A commonly used model by Zuber [41] to calculate CHF is given by Using properties of the working fluid in the simulations (ρ l = 200 kg/m 3 , ρ v = 5 kg/m 3 , h lv = 1 × 10 4 J/kg, and surface tension σ = 0.1 N/m), the critical heat flux is calculated to be 1.0891 × 10 5 W/m 2 . This estimate is at the same order of magnitude as the maximum critical heat flux calculated in simulations. However, this correlation does not capture the effect of the contact angle on the CHF. Kandlikar [30] proposed a model for CHF that accounts for the effect of surface wettability: where κ is a surface-dependent parameter, which is large for a poor wetting surface but small for a strong wetting surface. It is given by: where θ is the contact angle and f is the heater orientation angle relative to the horizontal. It has been found that our simulation results are quantitatively smaller than the values from the model proposed by Kandlikar [30], but qualitatively the trends are similar and comparable favorably for hydrophobic surfaces (Fig 13). For hydrophilic surfaces, the agreement between the model and numerical simulations is less favorable, the numerical data are noisy due to significant fluctuations in transient heat flux. Actually, similar non-monotonic data are obtained from experiments for hydrophilic surfaces as seen in Fig 13. In comparison, greater fluctuation in CHF data is found to be related to liquid-vapor dynamics. CHF denotes the onset of transition boiling where there is greater tendency to form and break vapor film next to the wall. Since liquid prefers to remain in contact with a hydrophilic surface, there is greater tendency to destabilize the vapor film formation process. As a result, sometimes an asymmetric vapor film is formed (see Fig 6( Experimental data in the transition region of the boiling curve have been speculated to be "noisy" due to dynamic advancing and receding contact angles [9]. In this work we have static contact angles. From pooling boiling experiment [42], Ramanujapu and Dhir have shown that advancing and receding contact angles are within only ±5 deg of the static contact angle. For the first step of this study, the static contact angle model throughout the bubble evolution process is feasible. Yet, the fluctuation in the transition regime persists due to the presence of different unstable modes as discussed above. Some prior numerical results have shown that it can still captures key features of boiling without considering the pinning effect [43]. Thus, at the first step in this work, the pinning effect is not included in the simulations while being able to capture the primary mechanism.

Leidenfrost point.
The Leidenfrost point characterizes the onset of stable film boiling and minimum heat flux. Fig 12 shows the plot of minimum heat flux at different contact angles. Berenson's minimum heat flux model is given by [44]: From the fluid properties used in simulations, minimum heat flux is found to be 1.168 × 10 3 W/m 2 , which is close to our simulation values. Fig 14 shows the plot of the Leidenfrost point (LFP) at different contact angles determined from our simulation results. The trend from simulations is compared to experimental data by plotting normalized values of the Leidenfrost temperature [32,45]. In both cases, the LFP decreases with increasing contact angle; however, experimental data show a sharper decrease of LFP with contact angles in the hydrophilic regime when compared to simulations. This may be due to the effect of surface roughness. Consider two experimental cases-"textured surface" and "smooth surface" from a prior work [32], as shown in Fig 14. Smoother experimental surfaces (experimental surfaces are not perfectly smooth) are found to show a shallower plot for Leidenfrost temperature due to fewer nucleation sites. Simulation data in this work are for smooth surfaces without any nucleation models. In this case, it is expected that the vapor film will be less stable for hydrophilic walls as indicated by a shallower plot for Leidenfrost temperature in the hydrophilic regime. This suggests that smooth wall simulations may provide a baseline limiting case for the development of suitable models for LFP. 4.2.5 Temperature field. Following experimental results and Boltzmann simulation work [11,12,15,16,22], in this research the liquid-vapor interface is not constrained to be at the saturation temperature. Instead, an interfacial heat flux exchange model is adapted . Fig 15a-15d show temperature profiles with a gradient of temperatures across the liquid-vapor interface which is similar to prior LBM work [15,16]. Fig 15e is the temperature contour result of vapor bubble within the Lattice Boltzmann method which alo shows a higher temperature inside the bubble and gradient interface temperature. From experiment and other related research https://doi.org/10.1371/journal.pone.0187175.g015 [11,12,22], the temperature at the interface as well as inside the bubble is higher than T sat for growing bubble on the heating surface which is in agreement with our results. The corresponding boiling curves are found to have similar trends as expected. It was found that the imposition of saturation temperatures on the liquid-vapor interface destabilized the formation of the vapor film. Consequently, the transition to film boiling in the boiling curve was significantly delayed or not observed for the parameters we tested. A detailed investigation is not within the scope of this work but is warranted in future.

Conclusion
In this work, numerical simulations of evolving liquid-vapor interfaces during pool boiling on a horizontal smooth surface have been performed to study the surface wettability effect and related dynamics. Instead of a saturated temperature constrained interface, an interfacial heat flux exchange model at the interface is adapted. The simulation results based on the volumeof-fluid method and static contact angle model have been carried out and compared with some prior theoretical and experimental predictions, demonstrating good agreement. The effect of surface contact angle and superheat on the complete boiling heat transfer curve is obtained for the first time and the corresponding dynamics has been qualitatively captured. It is verified this approach can be used for investigating boiling phenomena and providing more physical insights into the corresponding dynamics. In addition, it can provide some guidance for more time consuming three dimensional cylindrical numerical simulations. In a near future, specific boiling regime will be focused and more physics behind the dynamics can be assured.