Nucleate boiling is used in numerous engineering applications, such as the chemical, manufacturing, thermal, nuclear, and electronic industries. This research paper deals with the numerical analysis of bubble growth using a fluid flow model. This physical phenomenon of bubble growth has not been discussed mechanically and does not throw light on empirical models. We are discussing this phenomenon in another way to get the required results in this paper. Simulation of bubble growth is already published by measuring the volume of fluid flow; the method is known as VOFF tracking method. Lee’s model has already discussed the phase change that occurs due to evaporation and condensation of the fluid. We have used the method in which the equation terms involving energy and mass source caused by phase change are incorporated into the control equations by additional subroutines written in C language. We have mentioned in detail the results of simulating mass transfer caused by phase change and the effect of subcooling on bubble growth. The results thus obtained show that the subcooling effect prevents the growth of bubbles from growing due to the certain amount of bubble caps in the subcooled area. The effect of evaporation of the liquid increases the size of the bubbles in both the subcooled and the superheated zone.

1. Introduction

The multiphase flow that involves the transfer of heat and mass is the most important industrial phenomenon, and its fields of application include power generation, internal combustion engines, cooling and refrigeration systems, etc. The characteristics of boiling, evaporation, and condensation properties have a major impact on the design and performance of various technical systems. All these show the necessity for adequate numerical analysis of the processes and the need to develop global forecasting tools that are applicable to a wide variety of multistage flow configurations. Computational Fluid Dynamics (CFD) has been very successful in solving single-stage problems, and its industrial applications are now advanced. In contrast, the development of the multilevel CFD process is still ongoing and its applications are mainly limited to simple configurations [1].

We usually deal with phase change heat transfer processes. Heat transfer in which change of phase occurs is called latent heat. A large amount of heat in the form of latent energy is transferred to nucleated boiling, making it an efficient heat transfer mechanism. Nuclear boiling has applications in many areas of engineering, including nuclear, chemical, and microelectronic cooling [24]. Although there are many applications in many fields, the science of this phenomenon is not very clear and varies from application to application. It is also called “unsafe” science by some researchers [5]. The uncertain world is used to reveal the dependence of phenomena on a wide range of parameters. These parameters include the nuclear site density, the overheated liquid and the evaporation of the microlayer, the surface properties, the heat flux, and other variables of the system so that a comprehensive model cannot be obtained [6].

In recent decades, many empirical and semiempirical correlations are proposed to predict the relative dynamics of bubbles [7, 8]. However, the proposed correlation is limited to its own experimental results and cannot predict the conclusions of other researchers with reasonable accuracy. Therefore, an alternative method of numerical modeling is needed to describe the kinetics of the firing of nuclei with a wide range of parameters. Fluid mechanics and nanofluids are widely used in various heat transfer problems to improve their heat transfer properties and hence their performance. As one of the thermophysical properties, the specific thermal capacity of nanofluids plays a major role in the heat transfer of thermal area using nanofluids. In this regard, various studies have been conducted to investigate the factors that affect the specific heat of nanofluids [9].

Many engineering fields like nuclear, thermal, manufacturing, chemical, and electronics industries includes nucleate boiling. The phase change phenomenon is involved in nucleate boiling for which mass and energy transfer are modeled at the interface. For this purpose, interface tracking algorithms are needed to be adopted to check the boundary layer between two phases. Level Set, Volume of Fluid (VOF), and a combination of both Level Set and VOF which is Coupled Level Set Volume of Fluid (CLSVOF) are the common interface tracking methods in simulations [1012]. These interface tracking methods are used by researchers to simulate nucleate boiling with their own devised phase change models. The first numerical study on nucleate boiling was done by Son et al. [13] who used the level set method with mass and energy balance at the interface.

Son et al. [14] followed the model of Son et al. [13] to simulate vertical bubble merger from a single nucleation site. Mukherjee and Dhir [15] also followed the model of Son et al. [13] to simulate lateral bubble merger from different nucleating sites in pool boiling, respectively. A numerical model of Mukherjee and Dhir [15] was further used for studying the impact of dynamic contact angle and hydrophilic surfaces on vapor bubble growth [16, 17]. The CLSVOF interface tracking method was used by Ling et al. [18] for simulating nucleate boiling and the merger of vapor bubbles. Heat flux was measured and compared with experimental correlations which showed good agreement. The CLSVOF interface tracking method and phase change model were further used by Ling and Tao [19] to educate the effect of heat transmission coefficient in shallow liquid which showed an upsurge in heat transmission coefficient. Kunkelmann and Stephan [20] were the first who simulated the growth of spherical bubbles by coupling VOF with the phase change model of Hardt and Wondra [21].

More refinement of mesh and more bubble ebullition cycles were recommended to capture evaporation and validate the results with the experiment, respectively. The literature review shows that the level set is mostly used by researchers for interface tracking with their own devised phase change models. However, it has been reported that the level set method is not volume conservative and may lead to diffusion of mass in simulating phenomena involving phase change [22]. In recent years, significant progress has been made in developing technologies that use VOFs to capture interface characteristics, especially surface tension [23]. Haelssig et al. [24] simulated two-dimensional (2D) countercurrent and phase transition using approximate mixing formulas to process near-interface cells. Ali et al. [25] and Sun et al. [26] simulated 2D boiling, taking into account the saturation of one of the phases, and used the temperature gradient in the interface unit by face value, thus sacrificing process accuracy.

In the following paper, condensation with 2D axisymmetric bubbles was also proposed. Akhtar and Kleis [27, 28] used a mixing formula to simulate the growth of bubbles in a three-dimensional (3D) static liquid, but the calculated shape of the bubble was distorted and the growth rate was nonmonotonous. Raya and Kandlikar [29] accomplished 2D axisymmetric imitation on the same datum and observed only a small deformation of the temperature boundary layer. In addition to utilizing mixed thermal diffusivity between the interface cells, their method is also used to catch the user interface in a completely clear way. Raya and Kandlikar have successfully applied their method to the 2D axisymmetric simulation of nucleated boiling [30].

Recent study has revealed its importance of calculating the phase velocity required for the advection of the volume fraction in a nondivergent manner in the presence of a phase transition. For this reason, Malan proposed a method of extrapolation [31, 32]. It is implemented in an analyzer that combines the VOF engineering method and the phase change model of the acute interface: the formula of the mixture is used to obtain the enthalpy advection. Zhou et al. [33] used a combined active and passive approach to improve the thermal performance of smooth pipes. For this purpose, an experimental treatment with constant heat flow on smooth tubes was performed. In addition, they investigated both the arc wire turbulator as a passive method and the bubble-injection method as an active method. Variable parameters in the study included the distance between the arc, the flow of water fluid (as working fluid), and the ratio of water-air to air-fluid mixture. The study also investigates experimentally the thermal energy release behavior of unchanged two-phase air and water flows in a smooth pipe containing an arc wire turbulator.

Gong and Cheng [34] proposed a Boltzmann model of heat transfer by changing the liquid-gas phase. The model used two particle distribution functions, the density distribution function and the temperature distribution function. A new form of the term source was derived in the energy equation, and a proposed pseudo-potential model was used in the proposed model to increase its numerical consistency. The commonly used Peng–Robinson equation of state was incorporated into the proposed model. The problem of growing bubbles and leaving the horizontal surface was solved numerically based on the proposed model.

Zhog et al. [35] studied the effect of bubble injection on the operation of heat exchangers with helical casing, and the tube was investigated. As a change in the intake air volume fraction, the effect of changing the intake air flow was investigated. It was revealed that energy damage is more pronounced in parallel flow than in counter current. The results showed that bubble injection improved the performance of the heat exchanger. The efficiency factor increased sharply as the air volume friction increased. This performance improvement is more pronounced in counterflow exchangers. Pattern estimation of two irregular fluids using Sato and Nishino enthalpy spread [36] was originally developed for ITM (IT in Mechanical) smears.

Malan [31, 32] showed that simulation of bubble growth can be performed in 3 D static liquid medium that is processed to identify and predict the behaviour performance. The authors demonstrated that their simulations converge to an analytical solution with 1st-order accuracy and that the bubble shape and the temperature boundary layer are at least perpendicular to one of the coordinate axes. As far as we know, in the works of Malan et al., for the first time, the 3D geometric VOF method with a shrill phase interface transition model is verified; no more complicated problems have been tried in this context note. Scapin et al. [37] also proposed a technique of extrapolating speed; however, their work proved to be integrated only in the VOF algebraic solution.

Many researchers are actively improving PSI-BOIL (Parallel SImulator of BOILing); one is an internal CFD problem-solver with Direct Numerical Simulation (DNS) function. It has been certified and verified for many issues concerning polyphase flow and interfacial mass transmission [36, 38, 39]. This study focuses on simulating nucleate boiling by coupling a phase change model with a conservative mass interface tracking method. To accomplish this, the VOF method is used for the interface tracking and it is coupled with a phase change model to simulate nucleate boiling. A 2D axisymmetric geometry is considered for simulating nucleate boiling. The results of simulating the mass transfer caused by phase change and the effect of subcooling on bubble growth are studied. Bureš and Sato [40] and Lee [41] studied sharp-interface phase change model for Cartesian by using the VOF method.

Simulations of the Sucking-problem, Stefan-problem, and bubble growth in superheated liquid were also performed. Results obtained showed that there is a considerable compatibility with analytical solution results. In bubble separation experiment, velocity discontinuity is the major problem caused by the difference between the mass-weighted velocity and the volume-weighted velocity caused by the interfacial phase transition. Guo et al. [42] did some special treatments in this regard by using the volume of fluid level set (VOSET) method. It is a new method of capturing the phase transition combines the advantages VOF and level-set methods for tracking the interface.

2. Methodology

First of all, we selected the governing equations to apply for finding the solution to the problem. A phase change model was developed to apply the mass and energy transfer. The simulation was done by setting the input domain parameters that were applied to the control volume.

2.1. Governing Equations

We will be using partial differentials to apply the governing equations to the system [43], because the mass source due to phase change is dependent of both vapour velocity and volume fraction. The mass conservation equation solved during the simulation is given by

In the previous equation, Sm refers to the mass source due to phase change. The volume fraction equation is solved only for the second phase. In this example under discussion, the second phase is steam and the volume fraction equation is written as

In equation (2) and present the volume fraction and velocity of the vapor phase. The volume fraction is a sum of both phases, and it must satisfy the following criterion of

Both the phases share the same equation of momentum which is solved for both phases in the simulation and is given by

CST in equation (4) refers to the surface tension force at the surface, and it is modeled as a volume force by the model of the continuous surface force of Brackbill et al. [22]. The value of CST is given for two phases and is given bywhere represents the surface tension coefficient and k interface curvature which is expressed in the form of normal vector divergence which is further given by

The equation of energy is modeled aswhere refers to the source of energy due to phase change. The values of other variables such as viscosity, density, and thermal conductivity in each cell of a domain are based on their volume fraction. Equation (9) presents a general expression for them in terms of which can be replaced by density, viscosity, and thermal conductivity.

3. Simulation Procedure

3.1. Phase Change Model

In equations (1) and (8), Sm and Se refer to energy and mass sources due to phase change and are incorporated in governing equations to model the energy and mass transfer. The phase change model of Lee [41] is used to calculate these source terms of mass and energy. The model adds these terms based on the volume fraction of a cell and its temperature. If a cell contains both vapor and liquid phases with temperature exceeding saturation temperature of the liquid. Evaporation of liquid would take place in the cell for which the mass source term is given by

Similarly, if cell temperature is less than the saturation temperature of the liquid, condensation of vapor would take place. The source term for condensation is given bywhere is the temperature of a cell containing an interface while represent volume fraction and density of respective phases with subscripts l and standing for liquid and vapor. is a relaxation parameter for controlling the rate of mass transfer between phases. The value should be taken as such so that the temperature in the region containing interface remains near saturation temperature of the fluid. For this purpose, the value of taken here is 100/s. The energy source e due to phase change process is evaluated through

3.2. Computational Domain

A 2D axisymmetric domain is considered to perform the simulation. The length of the domain is 2 mm and 4 mm in height. The imposed boundary conditions include constant surface temperature at the bottom, pressure outlet at top of the domain, axisymmetric conditions on the LHS, and no flux condition on the RHS as shown in Figure 1.

To save computational time, nucleation from the surface and the development of thermal boundary layer processes are bypassed. To perform these processes, an initial bubble of a radius of 0.1 mm is positioned in the domain, while a current boundary layer is patched between a domain. The temperature of the thermal boundary layer varies linearly from 395 K at 0 mm to 353K at 1 mm. The region above 1 mm height is kept at a subcooled temperature of 353 K. The angle of contact of the bubble with the wall Ѱ considered here is 110.

3.3. Working Fluid

The working fluid used in this simulation is water in atmospheric conditions. The properties of the water in different phases are given in Table 1.

3.4. User-Defined Function Implementations and Phase Change

The mass and energy transfer rate due to evaporation and condensation between the two phases must be modeled at the interface. To track the interface between the two phases, additional user-defined functions (UDFs) are written in C language. For the computational range shown in Figure 2, the UDF first traverses all cells through a loop c.

During looping, the UDFs retrieve the volume fraction of both phases. A cell that has a volume fraction () value of vapor phase between 0 and 1 is considered as an interfacial cell as shown in the right part of Figure 2. Based on the temperature of the unit mass, the terms of the energy source due to evaporation and condensation are calculated by equations (10)–(12) and implemented in the governing equations.

3.5. Solver Settings

Fluent 16.0 is used to solve the equations that govern mass, momentum, and energy. Pressure and speed are combined with the PISO algorithm. The upstream second-order scheme is used for discrete momentum and energy equations, while the Geo-reconstruct scheme is used for the reconstruction of the interface. Transient simulation is performed with variable time steps by using a courant number. Courant number is a dimensionless quantity that refers to information wave traversed in a time step through a given cell element and is mentioned in equation (13). A constant value of the courant number 0.25 is considered for simulation. Based on this value, the smallest time taken for a cell in the overall domain is considered as a time step. The convergence criterion is set to 1e-8 for the energy equation, while for all other equations, it is set to be 1e-4.

4. Results and Discussion

4.1. Grid Independence Test

At the initial stage, the domain is discretized with a structured grid having a coarse element size of 14 μm. The grid is further refined with a continuous decrease in an element size up to 8 µm as shown in Table 2. The growth of the bubble is observed both in lateral and vertical directions. Figure 3 shows the contours of the bubble that are drawn at time t = 10 ms. It is observed that the results of grids c and d are very close to each other. Grid d takes more time for simulation due to a large number of cells than grid c. Therefore, grid c is selected for further simulation to save computational time and resources.

4.2. Mass Transfer

The contours of mass transfer for the bubble growth period are shown in Figure 4. The region of the domain considered has a height of 0.7 mm and a width of 1.0 mm. At time t = 0 ms, the bubble completely contains in superheated liquid due to which evaporation occurs at the interface of the bubble. The evaporation rate depends on the temperature distribution at the field interface. In the initial stage, the maximum evaporation rate is observed and the growth of vapor bubbles occurs suddenly. As the bubble grows in the vertical direction, its upper portion proceeds towards the region near the saturation temperature. This leads to a decrease in the rate of evaporation which is observed at time t = 2 ms and continues up to time t = 5 ms. After time t = 5 ms, the bubble crosses the saturated temperature region of 373 K and enters the subcooled region. Condensation starts to take place at the interface of the bubble when, time t = 6 ms. The rate of condensation is the minimum at the initial stage due to the presence of a small portion of the bubble interface in the subcooled region. As the growth of bubble size continues, the rate of condensation also increases in the subcooled region. The condensation halts the growth of the bubble but still the size of the bubble increases. It is due to the presence of the major bubble interface portion in the superheated region as compared to the subcooled region. As a result of this, more evaporation occurs at the interface as compared to condensation.

5. Conclusion

We have discussed the numerical simulations of bubble growth and the effect of subcooling on bubble growth and found that the growth rate and the dimension of the bubbles are affected by the local temperature as the mass of the bubbles is increased or reduced due to the condensation and evaporation, respectively. The subcooling does not allow the bubbles to change the directions. The bubbles face the obstruction in moving, and it depends on the volumes of the bubble that exists in the subcooling area. The overheated and subcooled areas affect the size of the bubble due to the evaporation or condensation. It would be interesting to transfer the results of the numerical simulation process to microscopic simulations of boiling fluids in future numerical studies of multiscale boiling experiments and subsequent projects of boiling experiments with subcooled fluids.

Data Availability

All the data used to conduct this study are available in the manuscript in tabulated and graphical form.

Conflicts of Interest

The authors declare no conflicts of interest.


The authors greatly acknowledge the support of all colleagues and friends for the completion of this research.