Numerical investigation of boiling process in micropores

. The model for direct numerical simulation of the boiling process on porous surfaces was proposed and verified in the present paper. The VOF method using the CSF method to simulate surface tension was applied to solve this problem. The model was verified on the problem of pool boiling on a flat plate, where boiling patterns were obtained, and the heat transfer coefficient was compared with the experiment and analytical solution. Images of vapour bubble formation and dependence of heat transfer coefficient on temperature difference were obtained for a two-dimensional problem with uniform packing at various wall temperatures.


Introduction
Air or water cooling was suitable for the heat dissipation needs of many industries.Due to the increase in heat flux density, it became necessary to resort to various tricks to further using of single-phase cooling systems, such as increasing the thermal conductivity of the coolant [1], changes in the geometry of cooling systems and their overall optimization [2].Such improvements were primarily needed for power microelectronics, computer microelectronics, laser systems, and other industries.Also, the optimization of two-phase cooling leads to the reduction of already used devices such as air conditioner condensers [3] or refrigeration devices [4], leading to a reduction in their material capacity and price.
In the microelectronics industry heat pipes and evaporation chambers are mainly used for heat dissipation.These devices use a porous medium to transfer heat from the hot region to the cold region, which acts as a wick to transport working fluid from the condenser to the evaporator due to the capillary effect.The thermal resistance of the porous medium is an important obstacle to heat transfer into the evaporator section [5].Evaporation and condensation at the wickvapour interface are the main phase transition mechanisms in a heat pipe.At a certain heat flux, a bubbling regime of liquid boiling is formed in a wick, which results in a significant decrease in thermal resistance [6,7].Previous studies have shown that the coating of the surfaces with porous media causes an increase in the heat transfer coefficient for bubble boiling mode on a flat plate.However, the porous coating increases thermal resistance between the heated surface and vapour.Therefore, thermal resistance in the wick-substrate system should be reduced to improve wick performance.Also, an increase in heat transfer can be significantly affected by the next: an increase in the number of vaporization centres due to using the porous surface and an increase in heat transfer rate of phase transition due to an increase in the number and rate of bubble detachment from surface.To better explain the mechanisms involved, the mechanism of vapour bubble growth in a porous structure needs to be better understood.In several studies, for example [5], experimental investigations have been carried out for porous structures during boiling.
The backfill method and its porosity (anisotropic or isotropic) are also important.Thus, in paper [8] was demonstrated that anisotropic backfill with porosity in the range of 50-71% showed the best heat dissipation results.This is because vapour bubbles leave the vapour formation zone faster than for lower porosity.Isotropic backfill results in lateral growth of vapour bubbles, which subsequently form a continuous vapour region, thus negatively affecting heat transfer.There was also considered a variant of filling where individual particles formed a structure similar to a wick.With this formulation of the problem, the lateral interval of individual particles increased, and the vertical interval decreased.It was shown that using particles of 250 µm in diameter improved heat transfer compared to one with a diameter of 150 µm.Further increasing the diameter, however, did not improve the results.Variation of particle diameter for separately formed wicks was also considered.It was shown that wicks assembled from particles with a diameter of 50 µm had the best ability to conduct heat.Also in paper [8], the studies of the optimization of wick height of porous media have been carried out.As a result of calculations, the best results were shown for the value of height to width ratio of 4-5.
Many factors influence to behaviour of vapour bubbles: drag force, surface tension force, bubble size at detachment from the surface, pore size, porosity coefficient, and many others.All these parameters must be considered when modelling the boiling on a porous surface.To solve this problem, the VOF method with the CSF scheme for surface tension modelling was applied.The method has also been adapted for its use at a microscale.

Mathematical model
To date, the most widespread methods are those implementing the idea of continuous markers.Due to their efficiency and ease of implementation, among algorithms of continuous volume markers, the Volume of Fluid (VOF) method has gained the most popularity today, which has proven itself well for calculating free surface flows.The idea of this method is that liquid and gas are considered as a single two-component medium and the spatial distribution of phases within the computational domain is determined using a special marker function F(x,y,z,t).The volume fraction of the liquid phase in the computational cell is taken as F(x,y,z,t) = 0 if the cell is empty, F(x,y,z,t) = 1 if the cell is completely filled with liquid, and 0 < F(x,y,z,t) < 1 if phase boundary passes through the cell.
Since the free surface moves with the liquid, tracking the movement of the free boundary in space is carried out by solving the transport equation for the volume fraction of the liquid phase in the cell: where v is the velocity vector of the two-phase medium, found from the solution of a hydrodynamic equation system consisting of a mass conservation equation or continuity equation: and motion equations: ( ) where τ is tensor of viscous stresses, Fs is volume force vector, p is static pressure, and ρ is the density of twophase medium.The components of the viscous stress tensor τ ij are: where µ is the dynamic viscosity coefficient of twophase medium.
The density and the molecular viscosity of considered two-component medium are determined through volume fraction of liquid phase in cell using the mixing rule: (1 ) where ρ l , μ l are the density and the viscosity of the liquid, defined below, ρ v , μ v are the density and the viscosity of vapour.
Obtained values of density ρ and viscosity μ are used in motion equations and determine the physical properties of the two-phase medium.Special attention is paid to surface tension when considering liquid flows with a phase boundary.Therefore, one of the advantages of the VOF method is that it allows for relatively easy modelling of the effect of surface tension forces.
The CSF (Continuum Surface Force) algorithm, which involves introducing an additional volume force Fs into motion equations, was used to model surface tension in this problem.The magnitude of this force is determined from the relation: where σ is the coefficient of surface tension, k is the curvature of the free surface, which is determined as the divergence of normal vector: The normal to free surface is calculated as the gradient of the volume fraction of the liquid phase in the cell: The magnitude of the normal vector is determined by contact angle θ on the solid wall: cos sin where n w and τ w are normal and tangential to the wall vectors.
The energy equation is given by: where the enthalpy h is defined as: The effective thermal conductivity k eff is common to both phases: The term S h in this equation represents the volumetric heat source due to phase change at the interface: where lv m is the mass transfer rate from liquid to vapour, vl m is the mass transfer rate from vapour to liquid, and h lv is the latent heat of vaporization from liquid to vapour.
The mass transfer rate is implemented using the Lee model: where constant γ has dimension of 1/s and is defined as follows: where  is the accommodation coefficient, b d is the bubble diameter, � is the molar mass, R is the universal gas constant, sat T is the saturation temperature, and L is the latent heat.
After analyzing the works [9], [10], and [11], the constant γ in this work was assumed to be equal to 100 for both evaporation and condensation.

Model Verification
As a test, a numerical study of pool boiling on a flat surface with dimensions of 26.5×26.5 mm 2 was conducted.A constant heat flux density of 25 kW/m 2 , 50 kW/m 2 , 75 kW/m 2 , and 100 kW/m 2 was set on this surface.Lateral walls were adiabatic, and upper wall was an outlet.As a result of calculation, values of heat transfer coefficient at given heat flux densities were obtained and compared with results of experiment [4] and theoretical dependence.In case of bubble boiling of a stationary liquid in a large volume, heat transfer coefficient can be calculated using the formula: Here ν is kinematic viscosity coefficient, C p is heat capacity, r is heat of vaporization, k is thermal conductivity, σ is surface tension of liquid at saturation temperature T s ; ρ', ρ" is densities of liquid and vapor at T s , T s is saturation temperature in Kelvin.
The comparison of obtained results is presented in Figure 1, which shows the dependence of the heat transfer coefficient during boiling of pure water on heat flux density on the heater.As can be seen from this figure, the obtained results using the considered methodology are quite consistent with the calculation formula and experimental data [12].
Moreover, this methodology allows the tracking of nucleation and evolution of vapour bubbles dynamically.Isosurfaces at a vapour concentration of 0.5 at different times are presented in Figure 2, which visualizes the interface between the liquid and vapour.As can be seen from Figure 2, the used approach allows for modelling both the boiling process and dynamics of vapour bubble development, as well as obtaining integral and instantaneous values of various quantities and characteristics of the boiling process.Further calculations were carried out for the problem of boiling on the surface of a capillary-porous coating, which was modeled as a packing of spheres with specified sizes in a two-dimensional setting.

2D setting of problem
For a two-dimensional setting, a geometry was constructed with sphere diameters of 150 microns and a distance of 10 microns between them, forming a uniform packing, as shown in figure 3. The total size of the computational domain is 0.61×1 mm 2 .Periodic boundary conditions were set on lateral walls, a free outlet was set on the upper boundary, temperature and adhesion conditions were set on spheres, and adiabatic boundary conditions and adhesion conditions were set on the lower wall.Gravity force is directed from top to bottom.The contact angle was chosen to be 160º.The second-order QUICK scheme was used to approximate convective terms of hydrodynamics equations.An implicit first-order scheme was used to approximate unsteady terms of hydrodynamics equations.The second-order upwind scheme was used to approximate the energy equation.The VOF equations were solved using the Geo-Reconstruct scheme.Diffusion fluxes and source terms were approximated using a second-order scheme.The relationship between velocity and pressure fields was implemented using the Piso procedure.
Liquid and vapour considered as incompressible fluids, thermophysical properties were constant (at a temperature of 373K).
An unstructured mesh with a step of 1 micron and a total of 386,000 cells was constructed.At initial time, the domain was filled with the water at a temperature of 373 K. To correctly resolve water-vapour interface and ensure numerical stability, the Courant number was taken to be 0.2.
Figures 4 and 5 show the process of water boiling in sphere packing.At wall temperatures of 378 K (figure 4), small bubbles are only on the surface of the packing after the time of 2•10 -4 s, but for temperatures of 383 K (figure 5), vapour bubbles form throughout the entire packing and have a larger size.At a calculation time of 2•10 -3 s, the vapour almost fills the entire volume of packing for all wall temperatures.Figure 6 also shows the dependence of the heat transfer coefficient on the temperature difference.

Conclusion
A mathematical model was developed for direct numerical simulation of the boiling process on the surface of a capillary-porous coating.Such coating was modelled using a packing of spheres with specified sizes.The model was based on the VOF method for tracking the water-vapour interface and takes into account capillary effects.The Lee model was also used to describe phase transition and constants of this model were selected for the considered tasks.Debugging and test calculations were carried out to determine the minimum size of the computational domain, the detail of computational grids, time step, boundary conditions, and other numerical parameters to obtain correct results of boiling process modelling in the porous coating.
The model was verified on the task of pool boiling on a flat plate, where results obtained using the considered methodology are quite consistent with the calculation formula and experimental data.This methodology also allows the tracking of nucleation and evolution of vapour bubbles dynamically.
Images of vapour bubble formation were obtained for a two-dimensional problem with uniform packing at various wall temperatures.At a wall temperature of 378K, small bubbles form only on the surface of the packing at the initial time, while at a wall temperature of 383K, vapour bubbles form throughout the packing and have a larger size.When thermal balance is reached, almost the entire internal volume of packing is occupied by vapour.
The study was supported by a grant from the Russian Science Foundation No. 22-49-08018.

Fig. 1 .
Fig. 1.Dependence of heat transfer coefficient at boiling on heat flux density at the heater.