An Enhanced VOF Method Coupled with Heat Transfer and Phase Change to Characterise Bubble Detachment in Saturated Pool Boiling

The present numerical investigation identifies quantitative effects of fundamental controlling parameters on the detachment characteristics of isolated bubbles in cases of pool boiling in the nucleate boiling regime. For this purpose, an improved Volume of Fluid (VOF) approach, developed previously in the general framework of OpenFOAM Computational Fluid Dynamics (CFD) Toolbox, is further coupled with heat transfer and phase change. The predictions of the model are quantitatively verified against an existing analytical solution and experimental data in the literature. Following the model validation, four different series of parametric numerical experiments are performed, exploring the effect of the initial thermal boundary layer (ITBL) thickness for the case of saturated pool boiling of R113 as well as the effects of the surface wettability, wall superheat and gravity level for the cases of R113, R22 and R134a refrigerants. It is confirmed that the ITBL is a very important parameter in the bubble growth and detachment process. Furthermore, for all of the examined working fluids the bubble detachment characteristics seem to be significantly affected by the triple-line contact angle (i.e., the wettability of the heated plate) for equilibrium contact angles higher than 45◦. As expected, the simulations revealed that the heated wall superheat is very influential on the bubble growth and detachment process. Finally, besides the novelty of the numerical approach, a last finding is the fact that the effect of the gravity level variation in the bubble detachment time and the volume diminishes with the increase of the ambient pressure.

Therefore, in order to further improve the existing predictive tools, a deeper physical understanding

198
The conservation of energy balance is given by the following equation:

201
where cp is the bulk heat capacity, T the temperature field, and λ is the bulk thermal conductivity.

202
The source term on the right hand side of the equation represents the latent heat of evaporation.

203
The volume fraction α is advected by the flow field by the following equation:

205
Interface sharpening is very important in simulating two-phase flows of two immiscible fluids.

221
It should be mentioned that the VOF method in OpenFOAM does not solve Equation 4 222 implicitly, but instead applying a multidimensional universal limiter with explicit solution algorithm

229
As it is known, the VOF method usually suffers from non-physical spurious currents in the 230 interface region. These spurious velocities are due to errors in the calculation of the normal vectors 231 and the curvature of the interface that are used for the calculation of the interfacial forces. These errors 232 emerge from the fact that in the VOF method the interface is implicitly represented by the volume 233 fraction values that encounter sharp changes over a thin region [74].

234
As previously mentioned in the introduction section of the present paper, the VOF-based solver 235 that is used in the present investigation has been modified accordingly in order to account for an

284
The evaporating/condensing mass flux is calculated from Eq. (9) and must be incorporated into 288 289

290
This initial sharp source term field (SSTF) is integrated over the whole computational domain to 291 calculate the "Net Mass Flow" through the entire liquid-vapour interface, using the following 292 equation: 293 294

295
This value is important for the global mass conservation, in order to ensure that the magnitudes 296 of the mass sources in the liquid and vapour parts are equal and correspond to the net evaporation 297 rate. The sharp source term field is then smeared over several cells, by solving the following diffusion 298 equation for the smooth distribution of source terms 299 300

301
Δτ is an artificial time step and Neumann boundary conditions are imposed for the smooth 302 source term field on all boundaries of the domain. Therefore, the integral values of the sharp and the 303 smooth source fields remain the same, in spite of the smearing. The width of the smeared source term 304 field is proportional to the square root of the product of the diffusion constant "D" and the artificial 305 time step "Δτ". It should be mentioned that the value of "D" must be adjusted to the mesh resolution 306 such that the source term field is smeared over several cells.

307
Then, the source terms in all cells that do not contain pure liquid or vapour (α < 1-αcut and α > 308 αcut , where αcut may be set to 0.05) are artificially set to zero. This cropping step ensures that source 309 terms are shifted into the pure vapour and liquid cells only in the vicinity of the interface. The

310
interface therefore is not subjected to any source terms and is only transported by the calculated 311 velocity field. Therefore, the transport algorithm for the volume fraction field as well as the associated remaining source term field is scaled individually on the liquid and the vapour side through the 314 application of appropriate scaling coefficients. This scaling step ensures that the mass is globally 315 conserved and that the evaporating or condensing mass flow, corresponds globally to the net mass 316 flow through the interface.

317
The newly proposed scaling coefficients Nl and Nv are calculated by integrating the smooth 318 source term field in each of the pure phases and comparing it to the net mass flow m z{| (Equation

319
12), utilizing the following equations: 320 321 Finally, the final source term distribution is calculated using the above scaling factors in the 322 following equation:  (2) is discretised using a "Gauss upwind" scheme. The • term of Eq.

342
(4) is discretised using the "Gauss vanLeer" scheme, while the • 1 − 8 term is discretised 343 using the "Gauss interfaceCompression" scheme that ensures the boundedness of the calculated 344 volume fraction field. Finally, all Laplacian terms are discretised using the "Gauss Linear Corrected"     this situation has been derived by Scriven [80]. According to this analytical solution the bubble radius 363 as a function of time is given by the following equation:

366
where β is a growth constant, details of which can be found in the work of Scriven [80], and D is the

489
As it can be observed the numerical model predictions are in very good agreement with the 490 corresponding experimental data. The numerically predicted spatial and temporal evolution of the 491 generated bubble matches very well with the corresponding experimental images (Fig. 9). Some small 492 deviations in the shape of the bubble especially after its detachment from the heated plate can be 493 attributed to the fact that the proposed experimental images were recorded after a few bubble cycles,

494
while the numerical simulation images represent the first bubble cycle. However, as it is indicated in

529
The temporal and spatial evolution of the bubble growth and detachment process for the base 530 case of Table 4, is depicted indicatively in Fig. 10 (Table 4). Details regarding the overall runs conducted are 555 summarised in Table 7. As it can be observed a total number of nine additional simulations were performed changing 559 in each case the initial temperature field. The reference/base case in Table 7 corresponds to the 560 validation run of Fig. 9. The prescribed ITBL in each case is illustrated diagrammatically in Fig. 11,   561 where the initial variation of temperature with respect to the vertical distance from the heated plate 562 is plotted for each run of Series A numerical simulations.

565
The spatial evolution of the generated bubbles for each of the above cases at the time of 566 detachment, is depicted in Fig. 12. As it can be observed, there is a substantial increase in the bubble

573
The bubble detachment time with respect to the ITBL thickness is plotted in Fig. 13a, while the 574 equivalent bubble detachment diameter with respect to the ITBL thickness is plotted in Fig. 13b. It     (Tables 4, 5 and 6). Details regarding the overall runs are summarised in Table 8.

631
Approximately the same trend can be observed also in the bubble detachment time.

688
As it can be seen, in the case of the hydrophilic surface the shape of the bubble before its 689 detachment is more close to case B4 of the present investigation (Fig. 14) while the case of the  Table 9.

825
Examining the diagrams of Figure 27, it can be concluded that increasing the ambient pressure 826 level of the system from 1 bar to 5 bar it seems that the previously identified effects of gravity level 827 ( Figure 26) are diminishing. Furthermore, it is evident that in general, increasing the pressure level 828 the bubble detachment characteristics decrease significantly.

829
Finally, in order to compare the relative importance of the overall examined controlling 830 parameters in the bubble detachment characteristics,

839
As it can be observed, according to the overall parametric numerical simulations, the heated 840 plate superheat seems to be the most influential parameter in the bubble detachment characteristics.

841
Quite important is also the influence of the ITBL and the surface wettability. Finally, the gravitational   presented, validated and applied in the present investigation, constitutes a quite promising and novel