Modelling of pool boiling on the structured surfaces using Lattice Boltzmann method

. The process of boiling on spatially structured surfaces is simulated by a hybrid model based on Lattice Boltzmann Method and heat transfer equation. The model permits to study heat transfer at boiling in a wide range of surface superheats for different surface structural characteristics. The regimes of natural convection, nucleate boiling, and transition to film boiling are studied. The boiling curves for the surfaces with different structural and wetting properties are obtained. It was shown that the onset of nucleate boiling occurs at lower wall superheat on the structured surfaces than that on the smooth surface. However, at high wall superheats the heat flux and the critical heat flux at modified surfaces are lower. It was obtained that special modification of both structural and wetting properties of heat exchange surface permits to obtain higher removal heat flux as well as higher critical heat flux.


Introduction
Boiling is one of the most efficient mechanisms for heat transfer from a heated surface.Thus, it is widely used in various technological applications related to heat and mass transfer, e.g. in cooling systems.The main parameters characterizing the efficiency of heat transfer at boiling are the removal heat flux and the critical heat flux (CHF) at which the transition to film boiling occurs.In order to simultaneously enhance the removal heat flux and to increase the CHF, various methods of modifying the heat exchange surface are used which is connected with changing the surface wettability properties and introducing structural spatial modifications on it [1].
The surface can be modified by structuring using mechanical processing or laser texturing, and in some works, micro-and nanostructured or porous coatings are used to increase the efficiency of heat transfer [1,2].The addition of protrusions and depressions increases the number of vaporization centres.It leads to an earlier onset of nucleate boiling (ONB), permits to increase the CHF due to preventing vapour babbles merging, as well as to enhance the efficiency of heat transfer.
However, the search for the optimal configuration of a structured surface using experimental studies is a rather laborious task.First, in order to fully provide the desired surface characteristics, specialized equipment is required.There are works that consider all the morphological parameters of the surface at once, including porosity, roughness, and wettability.It is difficult to analyse the influence of each effect separately [3] due to a large number of parameters that need to be considered and a fast speed of the process occurred.
Along with experimental studies, methods for numerical simulation of the boiling process are currently being actively developed, which make it possible to more effectively determine the optimal surface parameters without carrying out expensive experimental studies.It has been proven that methods of computational fluid dynamics (CFD) are good alternative to experiment for studying two-phase flows [4,5].However, additional models should be used for tracing the gas-liquid interface.Moreover, they do not allow modelling of the nucleation process, which is important in obtaining the temperature of the onset of boiling, and the density of nucleation sites.Thus, it is impossible to calculate the boiling curves for an ensemble of bubbles that form on extended surfaces.
One of alternatives to CFD for gases and liquids flows is Lattice Boltzmann method (LBM).It is based on the calculation of the distribution function of pseudoparticles over discrete velocities and space.In recent years, this method has become popular in modelling liquid boiling [6][7][8][9][10].One of the features of the method is the possibility of modelling the process of vapour phase generation without specifying additional initial conditions based on empirical relationships.In addition, in LBM it is easy to set a solid boundary of an arbitrary shape with a different contact wetting angles.
The effect of contrast wettability on the smooth surface was considered in previous papers [9,10].The aim of the work is to determine the optimal configuration of a structured surface in terms of heat transfer enhancement at boiling using a hybrid model based on the Lattice Boltzmann method and solving the heat transfer equation in a two-phase medium in a wide range of wall superheats.The following model is used that was thoroughly described elsewhere [9,10].Only the main features of this model shall be presented below.In presented model the time evolution of the boiling process is described through discrete velocity distribution functions fi.These functions are calculated with the LBM approach by the following set of linear equations:

E3S Web of
where x is the position of the cell located on a straight regular grid, t is the dimensionless time, Ωi is the collision operator, Si source term for forces.In the model, two dimensional D2Q9 lattice is used, hence, discrete set of velocities is chosen: |ci| = (0, 1, √2), i = {1, 2 ... 5, 6 ... 9}.
where eq i f is the equilibrium distribution function, τ is the relaxation time.
The source term is calculated by the exact difference method [11,12]: where Δu = τF/ρ.The total force F acting on the particles consists of the following components: where FSC is the fluid-gas interaction force, Fsolid describes the fluid-gas interaction with the solid surface, and Fgravity is the gravity force.Said forces are determined as follows: where G,  are the interaction coefficients and ψ is a pseudopotential: where pressure is determined by the Peng-Robinson equation of state and s c x t is an acentric factor.In the calculations, the dimensionless variables are used for the values of temperature T, pressure P and density ρ.These dimensionless variables are presented in units Tc, Pc and ρc.These units correspond to the fluid parameters at the critical point.
Evolution of the temperature spatial distribution T(x,t) in the computational domain is calculated by the heat conduction equation taking into account diffusion, convection, work of pressure forces, as well as phase transition: where uf = u+0.5Δu is a physical fluid velocity.The cv and λ are the heat capacity and thermal conductivity coefficient.
The model algorithm is as follows: 0. The initial spatial distributions of density and velocity are set.In accordance with these distributions, discrete velocity distribution functions eq i f are calculated.
1.According to the forces spatial distribution (6) and the source term (4), the equations of collision and streaming (1) are calculated.2. The boundary conditions are enforced, which in this model are given by the Bounce Back approximation.3.According to the current distribution of the fluid physical velocity, the spatial distribution of temperature is calculated.4. Transition to step 1.This calculation is carried out until the desired boiling evolution is obtained.
Computational domain is shown in Fig. 1.Liquid/vapor medium are presented in red/blue colors while metal heater is presented in grey.Spatial and time steps are equal to Δx = 30×10 -6 m and Δt = 5×10 sec, respectively.In this work, 800 and 500 spatial cells in horizontal (axis x) and vertical direction (axis y) are used, correspondingly.The metal heater with thickness of 30 spatial cells is placed on the bottom boundary.Different rectangular steps made of the same metal were placed on its surface.Thus, we obtain the structured heat exchange surface with a number of rectangular caverns.In calculations, the height h and width l of the caverns as well as their number N on the heater surface can be modified.For example, five caverns, N = 5, with 40 spatial cells size, h = l = 1.2 mm, are presented in Fig. 1.At the left and the right boundaries of the solution region the periodic boundary conditions are applied.Constant temperature T0 = 0.9 Tc and pressure P0 in accordance with equation of state are specified on the top boundary.Inside of the metal heater only heat diffusion equation ( 6) is solved with constant temperature.Constant temperature Th is applied at the bottom of the metal heater.The calculations were performed for different values of wall superheat ΔT = Th -T0, which is equal to the temperature difference between the top and the bottom boundaries of the solution region.Thermal conductivity and heat capacity of the metal heater are set for copper.Heat capacity, thermal conductivity, and viscosity of the fluid are taken for water.

Results
The numerical simulations were performed for different configuration of the structured surface.The key parameters of the structures are the height and the width of the caverns, h and l, and their total number on the surface N. In this paper the width of the surface was nonviable parameter and was equal to 24 mm.Depending on the number N of the caverns on the surface, the pitch distance p between the caverns varies.At the first stage the simulation of pool boiling at the surface with square caverns uniformly distributed on it were performed.It was obtained that the gas phase nucleation occurs even at the caverns with small sizes, i.e. h = l = 0.12 mm corresponding to 4 spatial cells.In the paper, the regimes with h = 0.12, 0.24 and 0.48 mm were considered.The number of the caverns on the surface for these sizes are N = 100, 50 and 25, respectively.These regimes were compared with the boiling process on the smooth surface, i.e. h = 0 mm.All these calculations were made for the surface with neutral wetting properties, i.e. for 90º contact angle.Finally, the regime of complex modification of the surface was studied.The wetting properties of the bottom of the caverns was taken to be hydrophobic (110º contact angle) and the rest of the heating surface was hydrophilic (67º contact angle).

Boiling curves
At the initial stage, the boiling curves <q>(ΔT΄) were calculated, where ΔT΄ = <Tw>-T0 is the surface superheat and T0 is the saturation temperature (Fig. 2).Average heater surface temperature <Tw> was obtained by averaging the temperature Tw(x,t) over time and the heater surface.To obtain each point of the dependence <q>(ΔT΄), temperature value was set on the lower wall of the heater Th = const.In turn, the time and the surface averaged heat flux <q> through the metal heater at the height Hh = 0.9 mm was determined as <q> = -λh (<Tw>-Th)/Hh.In all considered regimes Hh = 0.9 mm was the bottom of the caverns.
As it was described previously [10], at low wall superheat the heat transfer occurs in the mode of singlephase natural convection.That is why the boiling curves <q>(ΔT΄) approximately do not depend on the surface configuration.After the onset of boiling, the slope of the <q>(ΔT΄) curves sharply increases, which corresponds to an increase in the intensity of heat transfer with the transition from natural convection to nucleate boiling.It can be seen that the temperature corresponding to ONB is lower for structured surfaces than that on the smooth surface.In the region of low surface superheat, the removal heat flux increases with increasing of size of the caverns However, the situation changes dramatically in the area of high wall superheats.The value of the CHF as well as the removal heat flux decrease significantly with an increase in the caverns size.It should be noted that after the exceeding the CHF all the curves fall on to the line corresponding to film boiling regime.
The obtained absolute values of the removal heat flux <q> are smaller than in experiments on water boiling on a smooth surface.It is connected both with the limitations of LBM (liquid to vapor ratio 30 in present work) and with the fact that in the calculations absolutely smooth surfaces are considered meanwhile in real experiments some roughness are always presented on the heater surface that could increase the heat flux due to some stabilization of the vapor bubbles and the presence of vaporization sites.

Boiling on the surfaces with evenly distributed square caverns
To describe such a behaviour of the boiling curves let us consider the density contour plots for described regimes.Let us compare the boiling process on the smooth surface and on the structured surface with N = 100 caverns for the same wall superheat ΔT = 0.09 Tc.The corresponding regimes are denoted in Fig. 2 by points A and B, respectively.In Fig. 3 the boiling process on the smooth neutral surface is presented.At the same wall superheat, the boiling process on the structured surface is completely different (see Fig. 4).The babbles are merged together possibly due to the presence of the caverns in which the vapour phase is easily formed.Large areas of the heater surface are covered with the vapour films preventing the heat transfer from the surface due to low coefficient of thermal conductivity of the vapour.It should be noted that the degradation of heat transfer on the modified surfaces compared to the unmodified lyophilic surface was also observed experimentally [13][14].

Boiling on specially modified surface
At the next stage, the simulation of pool boiling was performed for specially modified surface.Previous results have shown that at low superheats the enhancement of the ONB and intensifying of the removal heat flux occurs due to the presence of the caverns on the heat exchange surface.However, the number of caverns slightly influences these parameters.Moreover, at high surface superheats, the presence of caverns leads to a formation of the areas with vapour films.To avoid bubbles merging it was decided to increase the pitch size, i.e. the distance between neighbour caverns.The regime with N = 5 caverns with size h = l = 0.3 mm was considered.As it was shown and discussed in previous paper [10], the pitch size between the hydrophobic partens should be of the order of the bubble departure diameter to achieve heat transfer enhancement.In this section p = 4.8 mm was used, which is slightly exceeds mean value of the bubble departure diameter Dd ≈ 3 mm for our conditions.Moreover, in order to prevent the bubble merging on the top boundary of the heater, these areas were chosen to be hydrophilic.To reduce the temperature threshold of ONB, the bottom parts of the caverns were set to be hydrophobic.The obtained boiling curve for the modified surface is presented in Fig. 2. The density contour plot for point C in Fig. 2 is presented in Fig. 5.It is seen that the heat removal flux is higher than that on all previously considered regimes in the whole range of surface superheats.Moreover, the CHF on the specially modified surface is the highest.

Conclusion
With the help of the hybrid Lattice Boltzmann method the process of boiling on the structured surfaces was studied.Different spatial structure and wetting parameters of the heat exchange surface were investigated.Heat transfer enhancement and an increase in CHF at boiling on the structured surface with alternative wetting properties were demonstrated.To obtain optimal conditions for heat transfer enhancement further investigations will be performed.

Fig. 1 .
Fig. 1.Scheme of solution region.Red colour is liquid, blue colour is vapour, grey colour is metal heater.h and l are the height and the width of the caverns.N is number of caverns on the surface of the heater, Dd is bubble departure diameter, p is distance between caverns.

Fig. 5 .
Fig. 5.The density contour plot illustrating the boiling process at surface with specially distributed caverns.Point C in Fig. 2. ΔT = 0.09 Tc.