Numerical simulation and modeling of liquid film evaporation inside axisymmetric reentrant cavities

Evaporation of thin liquid films inside reentrant cavities occurs in several boiling processes where enhanced surfaces are utilized. In this work, evaporation from a single reentrant cavity with an additional thin channel is studied. The channel allows the backflow of liquid from the pool into the cavity during bubble growth. Direct numerical simulations were performed, showing a strong relation between flow to the film inside the cavity and bubble growth at the pore. Additionally, a model was created with a novel modeling approach which is based on solving the Young-Laplace equation. From the model characteristic nondimensional parameters can be obtained and the influence of geometry variations on hydrodynamics can be


Introduction
Surfaces with reentrant-type structures are well known to strongly enhance heat transfer in pool boiling.For reentrant-structures consisting of subsurface tunnels with pores, Nakayama et al. [1] suggested three boiling modes: the flooded mode, the suction evaporation mode, and the dried up mode.In the flooded mode, the cavities are completely filled with liquid.The suction evaporation mode is the desired mode of operation and is characterized by thin liquid films evaporating inside the tunnels.In the dried up mode, the tunnels are completely filled with vapor.Since then, many experimental studies confirmed Nakayama's observations but still there are many open questions concerning the dependence of the boiling modes on geometric, process, and fluid parameters.
Experiments with tubes with reentrant-type surface structures are performed since almost 30 years.A comprehensive review is given by Webb and Kim [2].Ayub and Bergles [3] varied the gap width of a GEWA-T b tubular surface and found an optimal gap width with highest heat transfer coefficients.Several authors confirmed that there should exist an optimal gap width or pore size for each boiling condition.Kim and Choi [4] introduced the concept of a "crossover characteristic", stating that surfaces with a small pore area perform best at low heat fluxes and surfaces with a large pore area perform best at high heat fluxes.In general, this trend can be observed for most reentranttype surface structures so far.
Single reentrant-cavities were manufactured for having stable nucleation sites in single bubble experiments or to create surfaces for chip cooling.Kubo et al. [5] manufactured single circular reentrant-cavities with pore diameters in the range from 1.6 ȝm to 3.1 ȝm.They studied the influence of degassing and subcooling on the boiling hysteresis and heat transfer coefficient with FC-72.Shoji and Takagi [6] manufactured cavities with a pore diameter of 100 ȝm in copper and performed measurements with water as working fluid.The cavities worked as stable nucleation sites and wall superheats were lower compared to conical cavities.Phadke et al. [7] created surfaces for chip cooling with an array of square reentrant cavities.The surfaces were mounted vertically allowing interaction of the cavities through rising bubbles.
As the review of Wörner [8] shows, numerical simulations have become a powerful tool to study boiling phenomena at single bubbles with good accuracy.In order to capture the evaporation of thin films at the bubble foot, subgrid scale models based on the work of Wayner et al. [9] were created.Numerical simulations of bubble growth at structured surfaces were performed by Lee and Son [10] as well as Lee et al. [11].The studied geometries include circular cylindrical cavities, multi step cavities and a combination of circular fins with cylindrical cavities.A level set method was employed to model the two phase flow with a fixed contact angle and 5 °C wall superheat.They found that the creation of thin films inside the corners of the cavities leads to a substantial increase in heat transfer coefficient.
First attempts in modeling boiling from surfaces with reentrant-type structures were undertaken by Nakayama et al. [12] , who split the total heat transfer into internal heat transfer inside the structure and external heat transfer on the outer surface.The internal heat flux is modelled as latent heat transfer of evaporating thin films and the external heat flux as sensible heat transfer caused by convection of departing bubbles.Similar approaches were chosen by Chien and Webb [13], Ramaswamy et al. [14] as well as Chen et al. [15] reducing the number of empirical constants and refining the calculation of heat transfer inside and outside of the structure.Focus of these models is heat transfer while modeling of the flow from and to the cavities is more empirical in nature with input from experimental measurements required.
As a first step towards improving the understanding of the boiling modes observed with subsurface tunnels, direct numerical simulations of single reentrant-cavities are performed.The geometry is chosen such that liquid films can exits in the cavities and be fed by liquid from the pool.From the observations with those simulations, a simplified numerical model with a novel modeling approach, focusing on hydrodynamics, is created.The model enables the investigation of a wide parameter range and analyzing characteristic parameters.

Direct numerical simulation
The model employed for numerical simulation is based on the work of Kunkelmann and Stephan [16] and is implemented in the open source CFD software package OpenFOAM.The basic incompressible conservation equations for mass, momentum and energy are solved with mass source terms accounting for phase change according to the method of Hardt and Wondra [17].Convection of the liquid-vapor interface is modeled with the Volume-of-Fluid method based on the work of Hirt and Nichols [18].Kunkelmann and Stephan [16] added an interface reconstruction such that calculation of curvature of the liquid-vapor interface and of the temperature gradient at the interface can be improved.A smoothing operation on the curvature reduces spurious currents by an order of magnitude compared to the original implementation in OpenFOAM.The temperature gradient at the interface is used to calculate the evaporation rate.Evaporation at the three phase contact line is calculated with a subgrid scale model, based on the work of Stephan and Busse [19].From the contact line model, the contact angle is calculated locally and transient for each cell containing a piece of contact line and employed in the CFD simulation.The temperature equation is solved in the fluid region and in the solid wall until convergence of the temperature at the fluid-solid interface is reached at each time step.
The full model was validated by Kunkelmann and Stephan [16] with experimental boiling data and by Herbert et al. [20] with single droplet impingement data.

Computational domain
Figure 1 shows the axisymmetric computational domain for the numerical simulation.The fluid region is coupled to a solid region at the wall.There are two small pieces of wall where an adiabatic wall boundary condition is applied, as the geometries differ at this location.In the fluid region, there is a small channel at this location, allowing liquid flow from the bulk into the cavity.The solid region extends across this channel to allow heat flow to the pore.A heat flux boundary condition is applied at the bottom of the solid region.The sides are adiabatic.The fluid region above the cavity extends 2 mm to the top where the pressure is fixed and 0.7 mm to the side, which has the properties of an adiabatic slip wall.CFD simulations were performed for refrigerant R-134a at T sat = 20 °C for heat fluxes of 10 kW/m 2 , 40 kW/m 2 , 80 kW/m 2 , 120 kW/m 2 and 240 kW/m 2 for one geometry.Properties are taken from the VDI heat atlas [21].For the solid wall, properties of pure copper are applied.Five geometries were studied with the dimensions given in Table 1.Calculations are performed until a quasi-steady state is reached.A grid dependence study was performed for the case R75W3 with a heat flux of 120 kW/m 2 .Dividing cell sizes in half leads to deviations less than 8 % for the wall superheat and less than 4 % for bubble frequency.

Simplified quasi-steady model
The simplified model of flow to the cavity is based on parameterizing and solving the Young-Laplace equation (1) in MATLAB.

K p p
The modeling approach is chosen based on the observation that inertia effects are minor even at high heat fluxes and that more than 95 % of the total heat is removed as latent heat with the bubble from the surface.
By solving the Young-Laplace equation, the volume of the bubble and the volume of the liquid film can be related with a pressure jump across the liquid-vapor interface of the bubble and of the film.From the pressure differences, the flow to the film can be calculated and from the heat flux together with volume conservation inside the cavity the bubble growth.The procedure is illustrated in Figure 2.
The volume change rate of the bubble is given by equation In here it is assumed that the heat flow Q is completely transferred as latent heat.The heat flow to the liquid film film Q can be any portion of Q .In this work the heat flow to the film is taken to be proportional to the contact line length of the liquid film to the total contact line length.This is just a rough estimate but the exact value of the distribution of heat flow has only very little influence as long as the density ratio l v / ȡ ȡ is small.With the volume change rates and the timestep Ĳ Δ the bubble and film volume can be updated explicitly.The bubble volume is reset once the bubble departure diameter is reached.If the film volume exceeds the maximal value possible, bubble departure is induced as well.If the film volume becomes negative, the calculation is aborted.The flow resistance ȟ is calculated in this work by solving the Navier-Stokes equation for the flow of liquid through a thin channel.

Calculation of bubble shape, pressure, and departure diameter
Parameterizing the Laplace equation in cylindrical coordinates leads to the system of equations ( 4) or in nondimensional form to equations (5).For each pore radius there is a maximal volume, the equation can be solved for.In this case, the local angle of the path as well as the contact angle reach 90 degrees at the pore.During bubble departure, the interface typically contracts above the pore, leaving some vapor at the pore.The size of this vapor rest is not known from the solution, in this work the bubble with the volume leading to the maximal pressure is taken to remain at the pore.This corresponds to a bubble shape close to half of a sphere with the radius of the pore.
The obtained departure volumes are nondimensionalized with the capillary length and follow the trend of the fit to experimental data by Bari and Robinson [22].Deviations are less than 5 % as shown in Figure 3.

Calculation of film shape and pressure
The interface shape inside the cavity can be calculated with a similar procedure as calculating the bubble shape.
For the film, the function ) , (  The contact angle is enforced on the vertical wall and if the upper end of the film is not pinned also on the horizontal part.For the upper contact line being pinned at the pore, the contact angle can be increased until again the vertical part inside the pore limits the contact angle.Once the film reaches this point, it is assumed that bubble departure is induced as was observed in the numerical simulations.
As shown in Figure 4, the pressure jump increases with decreasing film volume, allowing the system to be stable.For the pinned contact line, dependence of pressure on volume is lower.Consequently, if operated in this range, variations in film volume become larger.The contact angle has to be taken from the direct numerical simulation, as the relation between wall superheat and heat flux is not modeled currently.Modeling of wall superheat should be possible for cases being completely dominated by thin film evaporation inside the cavity and evaporation at the three phase contact line.As will be shown below, contact line evaporation plays an important role in this work, but other effects can not be neglected completely.There is only one liquid film inside the cavity increasing the relevance of forced convection on the outside and there is pumping action of cold liquid through the thin channel, which should be considered when calculating wall superheat.

Results from direct numerical simulation
Numerical simulations show that bubble growth at the pore and liquid flow through the thin channels to the cavity are strongly related.Liquid films become larger with lower heat flux and larger channel width W. Smaller pore radii lead to smaller departure diameters and consequently the pressure in the vapor stays higher and liquid films thinner.The evolution of liquid in the film follows the bubble growth.The variation in film volume during one bubble cycle is large for low heat fluxes and large channel widths.If the channel is very thin, the film volume stays almost constant during the bubble cycle.
Depending on heat flux and cavity dimensions three different modes of bubble departure were observed, which are illustrated in Figure 6.On the very left, the bubble departs after reaching the maximum volume possible for the given pore radius.The bubble contracts at the bubble foot and the contact line at the bubble foot moves to some extend into the pore.The liquid film is hardly affected by this procedure.This type of bubble departure can be observed for higher heat fluxes or narrow channels for liquid backflow.
The central picture shows a bubble departure induced by the liquid film.In this case liquid is flowing to the film until it extends upwards through the pore and reaches the bubble foot.The contraction appears inside the cavity leading to a large residual of liquid in the cavity after departure.As film and bubble connect during departure, liquid can flow not only through the thin channel but also through the pore into the cavity.Due to the early departure, bubble volumes are smaller than the maximum volume possible for the given pore radius.
On the very right, bubble departure with a shallow cavity is shown.The liquid film grows and reaches the bottom of the geometry, leading to immediate flooding of the cavity.
Contact angles smaller than 90° have no effect on bubble departure diameters but the pressure difference across the liquid-vapor interface of the liquid film is strongly affected.For a contact angle of 45° the interface of the liquid film becomes flat in the corner.As the pressure jump across the liquid-vapor interface at the bubble will be larger than the hydrostatic head from the liquid pool to the liquid film for most cases, no liquid film can exist.This demonstrates that contact angles lower than 45° are required for liquid flow into the cavity.In the numerical simulation a contact angle of 42.5° is obtained for the case of a heat flux of 240 kW/m 2 , which is large enough to keep the film very thin in the corner at the channel.
The simulations show that 95 % to 98 % of the heat is transferred as latent heat from the surface.Up to a distance of 0.5 ȝm from the three phase contact line 53 % to 69 % of the total heat are transferred.If the full numeric cell containing the three phase contact line is considered, this value increases up to 93 %.Generally, the portion of heat transferred in the vicinity of the three phase contact line increases with decreasing heat flux.As the length of the contact line increases with increasing pore radius and cavity radius, the heat transfer coefficient increases as well.Larger cavities require more space and consequently fewer cavities can be placed on a surface.If the heat transfer coefficient is related to the area required by the cavity, small cavity diameters perform better than large cavity diameters.
Figure 7 shows the heat transfer coefficient for different heat fluxes with the geometry R100W10.It shall be noted that the absolute value of the heat transfer coefficient h depends on the size of the boundary where the heat flux boundary condition is applied.As this value can be chosen rather arbitrary with a single cavity, conclusions should be drawn only from qualitative observations.In the simulations, the heat transfer coefficient decreases with increasing heat flux.This can be explained with the strong influence, the three phase contact line exerts on the heat transfer.The local heat transfer coefficient in the vicinity of the three phase contact line decreases with increasing heat flux, as the contact angle increases.This affects the heat transfer from the full cavity as the contact line plays such a dominant role.

Results from simplified quasi-steady model
Without heat transfer being modeled, the simplified model allows to study the influence of cavity dimensions, fluid properties and contact angle on hydrodynamics.
The nondimensionalization of the system gives characteristic parameters which govern the flow.The parameters are given in Table 2.
The nondimensional parameters together with the equations show some interesting characteristics of the system.Flow resistance to the liquid film and heat flow HEAT 2014 01005-p.5 are interchangeable.An increase in flow resistance has the same effect on the flow behavior as increasing the heat flow.The interchangeability of flow resistance and heat flow corresponds to the experimentally observed "crossover characteristic".It shall be noted that in a real boiling process an increase in heat flow will lead to an increase in wall superheat which will affect the contact angle.In the simplified model, the contact angle can be controlled independently from the heat flow.Furthermore, the heat flow appears in the nondimensional time scale, which means that an increase in heat flow will accelerate the processes, which in return will affect convective heat transfer.
The influence of the geometric parameters pore radius, cavity radius and depth of the cavity are related to the capillary length of the fluid.From the equations it can be taken that changing the density ratio has a similar effect as changing the flow resistance.
The relation of pore radius with capillary length indicates that the optimal pore diameter should decrease with increasing pressure and that pore diameters for water must be much larger than for other refrigerants, which was also suggested by Nakayama et al. [12].
A parametric study allows identifying the influence of geometry and fluid properties on the behavior of the liquid film.Figure 8 shows the variation of the maximal film depth in terms of pore radius, heat flow and contact angle.As was shown with the direct numerical simulations, if the liquid film reaches the bottom of the cavity, flooding might occur.A film extending far down furthermore shows that pressure in the cavity becomes low and liquid can easily flow to the film, making the system more prone to flooding.
An increase of the parameter * ) ( ȟ Q either stands for an increase in heat flow or an increase in flow resistance of the liquid film.An increase in contact angle can also be interpreted as an increase in wall superheat.Filled markers in red color indicate a condition, where bubble departure was induced by the liquid film connecting with the bubble interface through the pore, leading to early departure.
As expected, liquid films become smaller with an increase in heat flux or contact angle.While an increase in heat flux increases the pressure drop for flow through the channel, an increase in contact angle decreases the pressure drop across the liquid-vapor interface and thus the pressure difference between pool and liquid film.
For the lowest heat flux, the liquid film reaches its maximum size for all pore dimensions.In those cases bubble growth is slow enough for the pressure drop across the channel to become negligible.Consequently the curvature of the liquid film follows closely the curvature of the bubble.
For 1 ) ( = * ȟ Q and all contact angles, the maximal depth of the film first increases with the pore radius and than decreases again.This is caused by the counteracting of two effects influencing the maximal depth.Generally, a small pore allows larger film volumes because there is more space between the pore and the cavity wall.But with a small pore radius bubble departure diameters are smaller, leading to higher pressures inside the cavity.At small pore radii the flow to the film is limited by the large pressure in the vapor, at large pore radii, the maximal film size is limited by the space between pore corner and cavity corner.This effect diminishes with rising heat flux and the pore radius corresponding to the largest liquid films increases.At the highest heat flux, the liquid film stays small over the full bubble cycle.
The influence of contact angle is similar to the influence of heat flux.Maximal film depths increase with decreasing contact angle.With increasing contact angle liquid films become smaller, until they disappear completely.

Conclusion
Direct numerical simulations were performed for a circular reentrant cavity, with a channel connecting a liquid film in the cavity with the liquid pool.The simulations show that concerning the flow, dynamic effects are small and bubble growth and liquid flow to the cavity are strongly related.
Three types of bubble departure have been observed depending on geometry and heat flux.For liquid films being small, the bubble contracts outside the cavity and departs with the largest diameter possible.Bubble departure can also be induced by the liquid film growing through the pore and connecting with the bubble liquidvapor interface.If the cavity is too shallow and the liquid film reaches the bottom, the cavity will be flooded completely.

01005-p.6
Heat transfer is governed by processes in the vicinity of the three phase contact line.The heat transfer coefficient therefore decreases with increasing heat flux.
For the same system a simplified quasi-steady model was created.The new modeling approach is able to depict flow to the cavity and bubble growth with similar accuracy as the direct numerical simulation.The model allows studying flooding conditions based on nondimensional characteristic parameters.From the dimensional analysis and parametric study, effects can be explained that were observed experimentally for reentrant-type surfaces.

Figure 2 .
Figure 2. Solution procedure.For an initial volume of the bubble bubble V and of liquid inside the cavity film V , the solution of the Young-Laplace equation gives the pressure difference between vapor and pool pool v p p − as well as between vapor and

t
is the coordinate along the path of the interface, y the vertical and r the radial coordinate.ϑ is the local angle of the path and can be converted into the current contact angle with method.Once the contact angle int Ĭ reaches a certain value Ĭ or r is equal to the pore radius, the bubble volume bubble V is the contact line is pinned.For the bubble foot radius being larger than the pore radius, are obtained for bubbles of all sizes and by interpolating it is possible to obtain a function

Figure 4 .
Figure 4. Pressure in the liquid film.

Figure 5
compares direct numerical simulation results (left) with the calculated shape from the full simplified model (right) for three different instances of time at a heat flux of 120 kW/m 2 .

Figure 6 .
Figure 6.Evolution of the liquid-vapor interface with time during bubble departure.Different modes of bubble departure.

Figure 7 .
Figure 7. Heat transfer coefficient over heat flux.

Figure 8 .
Figure 8. Maximal size of liquid film for variations of pore radius, channel flow resistance and contact angle.