Lattice Boltzmann simulation of bubble evolution at boiling on surfaces with different wettability

The hybrid Lattice Boltzmann method is adopted for a detailed study of surface wettability effects on evolution of vapor bubbles and temperature field of heat exchange surface at boiling. The simulation results show that the bubble departure diameter increases with an increase of static contact angle, and its value normalized by the capillary length approaches 3 for superlyophobic surface. In the range of contact angles of 110° – 129° the size of the dry area bounded by triple contact line greatly increases compared to the bubble departure diameter. At contact angles θ ⩾ 153° the dry area does not shrink even at the bubble departure stage and the typical film boiling regime is observed. It has been shown that the deterioration of wettability affects the evolution of the temperature field beneath the bubble, which leads to a significant change in the local heat transfer rate.


Introduction
Boiling is a common phenomenon widely used in practical applications, including energy generation and transformation technologies. However, some questions related to the theoretical description of this process, including the influence of wettability of the heated wall on the local and integral characteristics of boiling heat transfer are still open. Experimental studies [1][2][3][4][5][6][7][8] demonstrate that usage of hydrophobic surfaces reduces superheat required for the boiling onset, significantly influences the evolution of vapor bubbles and leads to substantial decrease of the critical heat flux at water boiling in comparison to hydrophilic surfaces. A number of studies [3,4] have reported dramatic reduction of the boiling heat transfer intensity for water on superhydrophobic surfaces, while other reports [5][6][7] indicate heat transfer enhancement at low heat fluxes. Recent papers [9][10][11] present experimental studies of the pool boiling processes on the surfaces with mixed wettability (biphilic) fabricated to achieve the best heat transfer performance. In this regard, numerical simulations provide a unique opportunity to study the effect of the surface wettability on the boiling process with changing the free surface energy in a wide range.
A classic approach to two-phase flows simulation, including boiling, comes down to a numerical solution of Navier-Stokes equations with an additional interface tracking technique that requires introduction of external parameters. The Lattice Boltzmann Method (LBM) provides an alternative approach to fluid dynamics simulations based on the solution of discrete Boltzmann equation. The kinetic nature of this method allows simulating complex systems, including the emergence and transformation of interphase boundaries and interaction with solid surfaces [12]. The LB simulations of bubble behavior at liquid boiling and the effect of gravity, wetting contact angle, and Jakob number were successfully performed in [13,14]. In [15,16] the effect of heater surface wettability on boiling heat transfer and development of boiling crisis were studied using a hybrid thermal LBM. The present study is devoted to LB simulation of evolution of a single vapor bubble and triple contact line, and local heat transfer rate at nucleate boiling on surfaces with different wettability in a wide range of wetting contact angles.

Model
A detailed description of the numerical model used in the work is presented in the recent paper [16]. The hybrid model is based on the lattice Boltzmann method and the heat conduction equation in a two-phase medium. Peng-Robinson equation of state is used. The interaction of the liquid-vapor medium with a solid surface is taken in the form [15]. The surface wettability is usually quantified by the static contact angle θ, which is obtained at the equilibrium between the interfacial tensions acting at liquid-solid-vapor contact interfaces (often measured on a sessile drop on the surface). Fig. 1 shows the LBM simulation results for liquid droplets on the surfaces with the wetting contact angles of 43°, 90° and 129°, respectively. For simulation of bubble evolution on lyophilic surfaces, 160×250 mesh was used with a heater of 80×20 size placed in the middle of bottom boundary. For lyophobic surfaces we used computational domains with the same height and width of up to 1600 (the corresponding heater width of up to 1440). We applied periodic boundary conditions on the left and the right borders, bounceback condition on the bottom and the heater, and constant pressure on the top. The computational domain was initially filled with liquid, and the initial temperature was set to T0 = 0.9 Tc. The temperature at the bottom border of the heater was specified to be constant T = T0+ΔT and on the top border was set to T0. An insulating boundary condition was applied in the remaining part of the bottom border. Spatial and temporal steps in our simulations were h = 5×10 -5 m and Δt = 1.5×10 -5 s, respectively. Thermal conductivity and heat capacity of the heater were set to 5 W/(m·K) and 10 6 J/m 3 , respectively.

Results and discussion
At the first stage, the dynamics of a single bubble growth and departure at boiling on lyophilic surfaces (θ < 90°) was numerically simulated. The process of the boiling on lyophilic surfaces has been well studied, and a number of experimental results, available in literature for water boiling on surfaces from different materials with θ ranging within 30° -70°, provide a reliable benchmark for the model performance evaluation.   Fig. 2 illustrates LBM simulation results for bubble evolution at Ja = 1.0 and θ = 43°, where Ja = cpΔTρl/(rρv) is the Jakob number. The simulation demonstrates that after enabling the heat source, the liquid layer adjoint to the surface is heated up to vapor bubble appearance. Recent experimental study [17] and numerical modeling [18] of boiling on lyophilic surfaces distinguish two following stages: bubble growth and departure. The first stage is characterized by an increase of both the bubble diameter and the area bounded by the triple contact line (Fig. 2a,b). When buoyancy force exceeds surface tension force, which keeps the bubble on the substrate, the stage of bubble departure begins. At this stage the triple contact line, outlining the dry spot, shrinks (Fig. 2c,d) while the diameter of bubble base remains nearly constant. This process is accompanied by elongation of the bubble in the direction opposite to the gravity force. Collapse of the contact line (dry spot rewetting) indicates a complete bubble detachment. After it the liquid adjoint to the surface is heated up the appearance of a new bubble, and the abovedescribed process is repeated. The simulation results illustrated in Fig. 2 are in a qualitative agreement with experimental data obtained by using high speed visualization of single bubble evolution at boiling [17,19]. The evolution of bubbles at boiling on substrates with θ equal to 90° and 129° is illustrated in Fig. 3. Time is measured from the moment of bubble separation from the surface. The first frame in Fig. 3a depicts the density distribution at the moment after bubble detachment from the surface with θ = 90°. As shown in Fig. 3a, the dry area bounded by the triple contact line does not collapse completely after bubble departure, in contrast to lyophilic surfaces (Fig.2), and vapor remains on the surface. Therefore, the waiting time of the next bubble appearance on the surface, tw, reduces to zero that is in agreement with experimental observations [1,2]. As depicted in the following frames in Fig.3a, in the first ~100 ms after bubble detachment the dry area under the bubble on the surface rapidly expands, while the height of the bubble remains approximately the same, which suggests the increase of the apparent contact angle θa. After the stabilization of the contact line, the increase of the dry area is accompanied with the increase of the bubble height. At some moment, the motion of the triple contact line stops, and the bubble growth proceeds through an increase of its height along the direction of buoyancy force (t = 645 ms, Fig. 3a). When buoyancy force exceeds the surface tension, the process of bubble detachment begins. This stage is characterized by gradual shrinkage of the dry area (t = 975 ms and t = 1035 ms, Fig. 3a). The main difference of these results from ones obtained for lyophilic surface (Fig. 2) is that the 4 apparent contact angle remains larger than 90°, and the interface boundary close to the triple contact line has a distinct concave shape. As a result, during the further bubble evolution a thin neck is formed, and its rapid collapse leads to the detachment of a part of vapor as a distinct bubble. The simulation results for a contact angle of 110° show that, in general, the evolution of the bubbles is similar to the data presented in Fig. 3a. However, with an increase in the degree of lyophobicity, there is a noticeable increase in the size of the region bounded by the triple contact line. Fig. 3b shows the results obtained for boiling on the surface with contact angle θ = 129°. The size of the heater increases from 4 to 8 mm for more accurate measurements of the area bounded by the triple contact line. The first frame corresponds to the moment of time after the bubble separation. Within the next ~30 ms the vapor bubble in the center coalesces with neighboring bubbles formed at the heater periphery. As a result, a thin vapor film covering almost the entire heater surface is formed (t = 90 ms, Fig. 3b). During the following ~330 ms the bubble grows in the direction of buoyancy force, while the dry area under the bubble remains nearly the same (t = 420 ms, Fig. 3b). Analysis of the observed process suggests that the evolution of the vapor bubble on the surface with such wettability reminds the classic hydrodynamics problem of the development of Taylor instability. In this case, the hydrodynamic instability initiates the bubble departure stage, during which the dry area begins to rapidly decrease. During the process the shape of the bubble reminds an elongating "bell" (t = 765 ms, Fig. 3b). It is noticeable that the bubble dynamics is similar to evolution and detachment of a droplet from a liquid film covering the horizontal tube in falling film evaporators [20,21].   Fig. 4 the bubble departure diameter for lyophilic surface is close to the capillary length  and increases with an increase of the contact angle which are in agreement with experimental and numerical results [8,14,17,22]. For superlyophobic surface the value of Ddep/Λ approaches 3 which is consistent with the result of numerical simulation [14]. The size of drop detachment from a liquid film covering a horizontal tube at vapor condensation normalized by capillary length is also equal to 3 [20,21] that may be typical for two-phase systems whose evolution is determined by the development of Taylor instability. As shown in Fig. 4, the slope of the dependence of dry area size on the contact angle in the range of θ < 129° is higher compared to the slope of dependence of bubble departure diameter, and two curves intersect at θ = 60° -90°. In the range of contact angle of 110° < θ ≤ 129°, a sharp increase in the maximum size of the dry spot during bubble growth is observed. Moreover, while the maximum size of the dry region for surface with θ = 129° does not reach the size of the heater and the bubble departure stage is characterized by partial rewetting, for However, it should be noted that the value of Dds/Λ = 6.4 for superlyophobic heat exchange surface is not a limit and is defined by the heater size. The sharp increase in the size of the dry area relative to the bubble departure diameter, observed in simulations at 110°< θ ≤ 129°, may lead to a significant change of the heat transfer rate at transition from single-bubble to quasi-film boiling regime.
The simulation results on the evolution of temperature profile of the heating surfaces with different wettability at boiling are shown in Fig. 5. Surface superheat is normalized as follows T* = (Th -Tsat)/Tsat). The evolution of temperature profile of the heated surface with θ = 43° illustrated in Fig. 5a shows the existence of regions with different heat transfer rate during bubble growth. The central part with higher temperature corresponds to the area of dry spot formed under the bubble after its appearance. A region with lower temperature at the periphery corresponds to the area of triple contact line including liquid microlayer, where the maximum heat transfer rate is achieved. The simulation results are in a qualitative agreement with experimental observations [19,23] of the temperature field evolution under a single bubble at boiling on thin film heaters obtained using high speed IRthermography. For a surface with θ = 90° (Fig.5b) at the moment of bubble detachment the integral temperature of the heater is minimal. There is a remaining overheated area -dry spot in the region of nucleation site, which results from the presence of vapor remaining at the surface after bubble detachment, while near the triple contact line an apparent temperature drop is observed. At the stage of bubble growth, the surface temperature in the area of dry spot and the size of dry area increase. At the stage of bubble detachment, the temperature drops both in the area of triple contact line and dry spot due to cold liquid inflow and partial rewetting of the heated surface. The analysis of temperature fields of heat exchange surfaces with different lyophobicity at boiling shows that during the most part of bubble evolution cycle, the integral surface temperature beneath a growing bubble is noticeably higher compared with lyophilic surface. The experimentally observed enhancement of heat transfer at boiling on weakly-hydrophobic surfaces [5][6][7], therefore, is not associated with an increase of the local heat transfer rate in the area of a single vapor bubble, as it follows from the comparison of the surface temperature evolution at boiling on surfaces with different wettability obtained by LB simulation. In accordance with experimental results [5][6][7], heat transfer enhancement is mainly attributed to a significant increase of the nucleation site density and a decrease of the temperature required for onset of nucleation, which also follows from the results of LBM.

Conclusions
A detailed study of surface wettability effect on vapor bubble and temperature field evolution at boiling on lyophilic and lyophobic surfaces has been performed with using hybrid lattice Boltzmann simulation. The results show that at boiling on lyophobic surfaces the stage between bubble departure and new bubble nucleation disappears which is confirmed by experimental observations. Bubble departure diameter and size of dry spots beneath the bubbles increase with deterioration of surface wettability. In the range of contact angle of 110°<θ ≤129° there is a significant increase in the maximum size of the dry area bounded by the triple contact line relative to the bubble departure diameter.
The results of LBM simulations of the evolution of the heating surface temperature field demonstrate that enhancement of heat transfer observed experimentally for boiling on weakly-hydrophobic surfaces due to a decrease in the temperature of onset of boiling, and the corresponding increase of the nucleation site density, rather than an increase of the local heat transfer rate by a single vapor bubble. Moreover, the discrepancy of experimental results on the heat transfer rate for boiling on weakly-hydrophobic and superhydrophobic surfaces is associated with transition from single-bubble to quasi-film and film boiling regime with an increase in lyophobicity level, which is confirmed by the lattice Boltzmann simulation.