Abstract

In order to investigate single bubble evolution, a boiling phase change model in subcooled flow boiling is proposed in this paper, and VOF model combined with phase change model is adopted to simulate the single bubble growth and movement. The effects of flow velocity, liquid subcooling, wall superheat, and vapor-liquid contact angle are considered in this model. The predicted bubble growth curve agrees well with the experimental result. Based on the analysis of bubble shape evolution and temperature field, it is found that the average bubble growth rate, flow velocity, and dynamic contact angle have significant effect on the bubble shape evolution during the bubble growth and movement while the temperature gradient in superheated liquid does not change with bubble growing. The character of dynamic contact angle during bubble growth and movement is also obtained in different working condition.

1. Introduction

In order to understand and optimize the process of boiling heat transfer, it is vital to understand the mechanism of bubble nucleation, growth, and departure. As a complex process, the investigation of model for bubble formation and disappearance has gotten only limited success. In recent years, there are some investigations about single bubble growth model for nucleate boiling.

Many researches about the bubble growth have been done using numerical methods, such as Mei et al. [1], Welch [2], and Heo et al. [3]. Mei et al. [1] have performed a numerical analysis to study the bubble growth under the condition of saturated heterogeneous boiling. The simultaneous energy transfer between vapor bubble and liquid micro-layer or heating wall was considered having effects on the bubble growth. Welch [2] carried out a study on the growth of axis-symmetric vapor bubble. The moving unstructured grid was considered by using the interface tracking method with the combination of the finite volume method. The control volume continuity, momentum, and energy equations were modified to include surface tension and discontinuous pressure and velocity. The flow-directional local grid (MPS-MAFL) method has been used by Heo et al. [3] to do the numerical study about the growth of bubble in transient pool boiling through moving particle semi-implicit with meshless advection. The growth process of a bubble with different initial radii was calculated under the condition of high heat flux and high subcooling condition.

Besides, some new models were used to study the bubble growth. A new understanding was provided by Li et al. [4] to study the bubble which was about the interfacial transport characteristics of inviscid spherical bubble with different geometric parameters, rising in a stagnant hot or bisolution liquid. The flow and temperature fields around bubbles and similarly sized rigid spheroids were computed numerically while the development of the physical model for vapor bubble growth in the condition of saturated boiling was provided by Liao et al. [5], including heat transfer through the micro-layer and the surrounding foam body superheated liquid and bubble growth. Both asymptotic and numerical methods were carried out to study the liquid temperature field surrounding a hemispherical bubble and it indicated that there was a thin unsteady thermal boundary layer existing adjacent to the bubble dome. During the early stages of bubble growth, heat transfer to the bubble dome through the unsteady thermal boundary layer constituted a substantial contribution to vapor bubble growth.

The region around a single growing bubble was considered to be subdivided into three subregions by Genske and Stephan [6], which were micro-region, bubble area, and surrounding liquid. Their results showed that the flow pattern in the liquid around a growing vapor bubble was determined by not only the movement of the bubble surface, but also the vapor flow that flied inside the bubble. Heat conduction would be the dominant factor on the bubble growth in the regions away from the bubble. Except that, the contact angle has effects on bubble growing. A static contact angle model and a dynamic contact angle model were proposed by Mukherjee and Kandlikar [7], which were based on the contact line velocity and the sign of the contact line velocity which was used to sign the dynamic contact angle model on the heating wall. The effects of dynamic contact angle on bubble dynamics and vapor volume growth rate were compared with results obtained with the static contact angle model.

The heterogeneous boiling on a horizontal plate in stagnant and slowly flowing fluid was simulated by Hazi and Markus [8] by using the method of lattice-Boltzmann approach. The study has found that, in a stagnant fluid, the bubble departure diameter was proportional to while the release frequency scales with where was the gravitational acceleration. According to the simulation results, it is found that there was no correlation between the bubble departure diameter and the static contact angle, but with the increase of the static contact angle, the release frequency increases exponentially.

Numerical researches about bubble growth are mostly based on the pool boiling. However, the character of bubble growth in subcooled flow boiling is quite different from that of the pool boiling. The bubble growth process in subcooled flow boiling is affected by the thermal characteristics of heating wall and mainstream liquid and so forth. Due to the existence of the mainstream liquid velocity, the motion characteristic of bubble interface and the flow field around the bubble are more complex than those of pool boiling. Thus, distortion and deformation of the bubble may occur which will lead to characterized features during the growth.

In the present paper, based on the characteristic of bubble growth process in subcooled flow boiling, a bubble growth model is proposed and a numerical simulation is performed to understand the process of bubble growth. The simulation results reflect bubble growth process as well as a dynamic result of evaporation and condensation which agree the experimental results very well and the relationship between the vapor-liquid interface motion and velocity field inside the bubble presented here. The character of dynamic contact angle during bubble growth and movement process is also obtained in different working condition.

2. Numerical Simulation Methods

The VOF formulation relies on the fact that two or more fluids (or phases) are not interpenetrating. For each additional phase that was added to model, a variable is introduced: the volume fraction of the phase in the computational cell. In each control volume, the volume fractions of all phases sum to unity. The fields for all variables and properties are shared by the phases and represent volume-averaged values, as long as the volume fraction of each of the phases is known at each location. This paper investigates the bubble growth process by VOF model [9, 10] and UDF interface of the CFD software Fluent.

2.1. VOF Model

Generally, a VOF algorithm [10] solves the problem of updating the phase volume fraction field given the fixed grid, the velocity field, and the phase volume fraction are determined in previous time step. In two-dimensional problem, the interface is considered to be a continuous and piecewise smooth line. The problem is reduced to the reconstruction of an approximation of the interface in each cell, knowing only the volume fraction of each phase in the cell itself and in the neighboring cells.

During all simulation cases in present work, a piecewise linear interface calculation (PLIC) [12] interface reconstruction method has been used for interpolation in a cell. In the existing CFD code, this scheme is the most accurate one and it is applicable for general unstructured meshes as used here; this interpolation scheme assumes that the interface between two fluids has a linear slope within each cell and this linear shape is used for the calculation of the advection of the fluid through the cell interfaces.

In VOF model, the th fluid’s volume fraction in the cell is denoted as and then the following three conditions are possible:: the cell does not contain the th fluid;: the cell is full of the th fluid;: the cell contains the interface between the th fluid and one or more other fluids.In present model, there are two phases: vapor and liquid. When , it represents vapor bubble region; when , it is liquid region. When , , there is a bubble interface exiting in this region. Based on the local value of , the appropriate properties and variables will be assigned to each control volume within the domain.

2.1.1. Continuity Equation

In VOF model, the volume fraction of each fluid is calculated by tracking the interface between different phases throughout the solution domain. Tracking of the interfaces between different phases present in the system is accomplished by solving continuity equations of the phase volume fraction for phases. A physical interpretation can be given by various terms in this continuity equation. For the vapor phase, this equation has the following form:where is the mass source term of vapor in the bubble growth process, which can be estimated by (8) and (10) through the UDF interface.

2.1.2. Momentum Equations

where is the surface tension force, which can be estimated through CSF (continuum surface force) model [13] as follows:

2.1.3. Energy Equation

The physical properties of and (effective thermal conductivity) are shared by both of the two phases. Source term contains contributions from all other volumetric heat sources.

2.1.4. Vapor-Liquid Interfacial Density

The properties appearing in the transport equations are determined by the presence of the component phases in each control volume [9]. In a two-phase system, for example, if the phases are represented by the vapor and liquid and if the volume fraction of the second phase has been tracked, the density in each cell would be given by

2.2. UDF Procedures

The User Defined Functions (UDF) interface of the commercial CFD code of Fluent [9] is employed to simulate the heat and mass transfer of the phase change. UDF subroutine allows us to customize Fluent to fit particular modeling needs. UDF is a function that user program can be dynamically loaded with the Fluent solver to enhance the standard features of the code. They are defined by DEFINE macros which are supplied by Fluent Inc. They can access data from the Fluent solver using predefined macros and functions. Present UDF programs are based on following described bubble generation theories of micro-layer evaporation to simulate the bubble growth.

3. Bubble Growth Physical Models in Subcooled Flow Boiling

3.1. Bubbles Growth Process in Subcooled Flow Boiling

Figure 1 shows the schematic diagram of bubble growth process on heating wall in subcooled flow boiling. The channel is vertically arranged, and the fluid flows vertically from bottom in -direction. The temperature of the mainstream region is lower than that of the liquid saturation temperature corresponding to local pressure. is the temperature of single-side heated wall which is higher than . There is a temperature boundary layer close to heating wall. In this region, the liquid temperature is higher than saturation temperature but lower than the heating wall temperature, that is, . There is a liquid micro-layer region near the bottom of bubble whose temperature is higher than . is the contact angle of the bubble with the heating wall. The growth times sequence of bubble is , , and .

At the initial time, , bubble radius is relatively small. The whole bubble is submerged in superheating liquid layer. The heat and mass supplied for bubble growth are derived from the liquid micro-layer evaporation at the bottom of bubble and superheating liquid around bubble surface. As the bubble diameter increases, at the time of , the top of bubble passes through the superheating boundary layer and enters the subcooled mainstream region. Consequently, the condensation will occur at the top of the bubble. At this period, as the mass into the bubble caused by evaporation at the bottom and the interface exceeds the condensation mass at the bubble top, the bubble continues to grow up. At the time of , the bubble diameter is greatly larger. Because a large portion of bubble top region is submerged in subcooled mainstream liquid, condensation effect becomes significant and bubble growth rate will be gradually slow. When the condensation mass of vapor is approximately equivalent to the mass added to the bubble by evaporation at the bottom and the interface, the bubble growth rate will slow down to zero and bubble diameter remains unchanged.

3.2. Bubble Growth Model

Based on above analysis, there are four assumptions in current model: (i) there is a very small initial bubble existing on the heating wall in which the radii will be estimated through classic bubble nucleation theory; (ii) the vapor in the growing bubble is ideal gas; (iii) the temperature at the initial bubble growth time is equal to saturation temperature corresponding to the working pressure; (iv) the pressure in the initial bubble is based on Yong-Laplace equation:where is surface tension coefficient, N/m.

3.2.1. Micro-Layer Liquid Evaporation Model on the Bottom of Growing Bubble

Experimental results [14] showed that there is a liquid micro-layer existing at the bottom of bubble during bubble growth, and it is very important to the process. The convectional heat transfer through this layer can be ignored due to its small thickness. Thus, the heat flux through the micro-layer liquid can be estimated aswhere is the thermal conductivity of the liquid, W/(m·K); is the temperature of the heating wall, K; is the vapor temperature of growing bubble, K; is the micro-layer thickness under the bubble bottom, m.

The micro-layer thickness for growing bubble on nucleation site was calculated by [15]where is vapor density, kg/m3; is latent heat of vaporization, kJ/kg; is the specific heat at constant pressure, kJ/(kg·K); is liquid density, kg/m3; is wall superheated, K; is bubble diameter, m; Pr is the liquid Prandtl number.

The micro-layer thickness for sliding bubble on nucleation site was calculated by [16]where is liquid dynamic viscosity, m2/s; is the height of sliding bubble, m; is the velocity of bubble slide, m/s.

The mass transferring from the micro-liquid layer to vapor bubble can be estimated bywhere is the latent heat, kJ/kg; is the heat flux through the micro-layer, kW/m2. is the net mass flow rate per unit area, kg/(m2·s).

Because Fluent is the CFD software which is based on the FVM (finite volume method), the mass source which is added to the continuity equation is the mass flow rate per volume, kg/(m3·s),where is the mass transfer rate proportional to vapor bubble volume fraction in the th cell of vapor-liquid interface region. is the area of the th grid. is the volume of the th grid.

3.2.2. Interfacial Bubble Heat and Mass Transfer Model

Vapor-liquid interface is quite important for bubble growth. When the radius of bubble is relatively small, the whole bubble is submerged in superheated liquid layer. Heat will transfer to bubble from the superheated liquids through the bubble interface, and the diameter will be increased consequently. When the diameter of bubble increases to a certain scale, the bubble interfacial area will be divided into two parts; one is submerged in the superheated liquid layer where the heat would be transferred from the superheated liquid to the bubble; the other one is submerged in subcooled liquid where the vapor in the top of bubble would be condensed and the heat would be transferred from bubble to the subcooled liquid through the interfacial areas.

According to classical analysis of the evaporation and condensation [15], the net mass flux at the phases interface could be estimated by Hertz-Knudsen-Schrage Equation [16]:where is the net mass flow rate per unit area, kg/(m2·s); is the universal gas constant, J/(mol·K); is the mass per mole, (kg/mol); is the liquid temperature of the vapor-liquid interface, K; is the saturation liquid pressure corresponding to , Pa; is the vapor temperature within bubble, K; is the vapor pressure in the growing bubble, is the evaporation coefficient of the vapor-liquid interface, and is the condensation coefficient of the vapor-liquid interface, respectively.

Under the conditions of thermodynamic equilibrium, the mass transfer mechanism representing evaporation coefficient will be not distinguished from the ones of the condensation coefficient; and ones often have the simplification of [17]. For a curved interface, an analogous equation can be deduced [1820]:When , the net mass from liquid will enter into the bubble through the vapor-liquid interface; when , the net mass from bubble will enter into the liquid through the interface.

The heat flux due to the phase change in the th cell of vapor-liquid interface region can be defined byIt should be known that is the mass source term of the continuity equation and the heat flux is the energy source term of the energy equation according to VOF model, W/m3.

3.2.3. The Dynamic Contact Angle Model

The three-phase moving contact line is a fundamental problem in multiphase flow. A material contact angle (given by Young’s equation) can be defined at a position microscopically close to the contact line. Due to the hysteresis, there are many stable contact angles for a given system, of which the largest is called the advancing angle and the smallest the receding angle. When the contact line is moving, the dynamic contact angle differs from its equilibrium counterpart and can be beyond the range limited by the static advancing and receding angles . The apparent contact angle is generally used as an auxiliary concept in studying the contact line dynamics, which is the angle formed between the wall and a line tangent to the interface at a certain distance from the apparent contact line.

Contact angle is important for flow boiling, including static contact angle and dynamic contact angle. The static contact angle reflects wetting property of fluid. When one fluid on the wall was displaced by another, the relationship between the contact angle and wet property becomes very complex due to the contact angle hysteresis phenomenon, this time the contact angle named dynamic angle, which was influenced by fluid property, wall property, and the velocity of three-phase moving contact line. The contact angle was continued to be studied, while there are still no certain correlations for analysis [2123].

Dynamic contact angle at bubble base during nucleate pool boiling was discussed in many present papers, while dynamic contact angle in the flow boiling was less studied. Because of the quasi-steady drag in the flow direction and buoyancy force, the character dynamic contact angle in flow boiling is very different from that in pool boiling. Figure 2 shows the dynamic contact angle in flow boiling. and were the dynamic contact angles of bubble behavior on the flow boiling [24]; was the incline angle. In bubble detachment heating wall, the contact area at bottom of growing bubble reduces and the three-phase contact line moves towards vapor phase, the contact angle named advancing contact angle, , while it can be seen that the receding contact angle and advancing contact angle at different time are in the flow boiling system.

3.3. Geometry, Boundary Conditions, and Working Conditions

In present work, a single bubble growth process has been simulated in subcooled flow boiling in a narrow channel. Due to the symmetrical characteristic in the wide side of the bubbles, it can be simplified to a three-dimensional geometry structure Figure 3(a). In the simulation geometry structure model, the vertical rectangular narrow channel has an inlet section of 2 mm × 2 mm and a length of  mm.

Based on the assumption of the bubble growth model, there is a small bubble on heating wall at the initial stage of the simulation. The diameter of the small bubble is calculated by the Han and Griffith correlation equation [11]In this paper, is ranged from 20 μm to 50 μm in the subcooled flow boiling.

Mesh around the bubble has been refined to obtain the details of bubble growth; Figure 3(b) shows the mesh through the center of bubble in plane. According to mesh sensitivity testing, the number of meshes (120000) is suitable for simulation. The heating wall temperature keeps constant while the other wall is adiabatic. The inlet of the channel is set as velocity-inlet condition, and the outlet is set as pressure-outlet condition. The viscous model is the laminar model. Table 1 is the working conditions as flowing.

4. Results and Discussion

4.1. Bubble Growth Process
4.1.1. Bubble Shape Evolution during Bubble Growth and Detachment Process

The evolution of bubble shape in subcooled flow boiling would be impacted by a couple of factors. Figure 4 shows the evolution process on case . Figure 5(a) is the curve of bubble diameter in the time of bubble growth and detachment. Figure 5(b) is the comparison of simulation results with experimental results obtained by Chen et al. [25], and it is found that the bubble growth curves based on CFD simulation agree well with the experimental results with uncertainty less than ±25%. Figure 6 is typical images of bubble shape in bubble growth and detachment process in forced convective subcooled flow boiling [26].

As shown in Figure 4, there is a growing bubble at the nucleation site on the heating wall. The nucleation site coordinates are  m,  m. At the initial time of bubble growth (0.1–0.4 ms), due to high pressure in the bubble, bubble growth is in the stage of inertial control. A large amount of heat and mass is transferred from the micro-layer to the bubble. The bubble grows fast as a hemisphere. The value of bubble diameter increases from 0.2 mm to 0.54 mm. After 0.4 ms, the pressure inside bubble decreases quickly; bubble growth in the stage is the thermal diffusion control. As the condensation intensity gradually enhances at the top of bubble, bubble growth rate becomes much slower than that of the former. From 0.6 ms to 0.8 ms, bubble expands slowly in hemisphere shape. However, with increased bubble diameter, the force balance acting on the growing bubble has been broken. As a result, the shape of bubble evolves from a hemisphere to an ellipse at the time of  ms. At  ms, bubble diameter increases to the maximum value 0.71 mm. After that bubble begins to detach from nucleation point. From 1.2 ms to 2.6 ms, the shape of bubble evolves from ellipse to sphere and bubble diameter decreased gradually. Bubble remains as a rough sphere until it lifts off from the heating wall.

4.1.2. Pressure Distribution Evolution during Bubble Growth

Figure 7 is the pressure distribution of bubble growth at 0.4 ms and 1.0 ms on case . As shown in Figure 6, in the bubble growth process, the pressure is identical in different regions inside vapor bubble. The pressure inside the bubble is higher than that of the flow liquid around the bubble. Farther away from the bubble, the pressure of the flow liquid becomes much higher. At  ms, the value of gauge pressure inside bubble is 1930 Pa. Away from the top of bubble, there is a flow liquid region with the minimum pressure, whose gauge pressure value is 1880 Pa. At  ms, the value of gauge pressure inside bubble is 1539 Pa. The minimum pressure of flow liquid region is at the two sides of the growing bubble.

Figure 8 presents the pressure evolution inside a bubble during bubble growth and detachment on case . As a whole trend, the pressure inside the bubble is decreased continuously. From 0.1 ms to 0.6 ms, the bubble growth is in the inertial control stage. Bubble expands rapidly leading to the dramatic pressure drop inside bubble. From 0.6 ms to 1.2 ms, the bubble growth is in the thermal diffusion control stage. The bubble growth rate is relatively low and stable. The pressure inside the bubble fluctuates from 1500 Pa to 2000 Pa. From 1.2 ms to 2.6 ms, the bubble is in detaching. Because of the intensive condensation of a departing bubble caused by subcooled liquid, the bubble diameter dramatically decreases and the pressure inside bubble drops rapidly also.

4.1.3. Velocity Field of Growing Bubble

Figure 9 shows the velocity vector field at  ms,  ms on case . According to the analysis of the pressure filed around a bubble, eddies are formed by the pressure discrepancy around the bubble. As shown at  ms, eddies can be observed on both sides of the growing bubble. The eddy accelerates contraction of the vapor-liquid interface at bottom of bubble, which leads to bubble shape changing from ellipsoid to approximate sphere. As the eddies transfer subcooled liquid to the heating wall on both sides of bubble, the micro-convection in the vicinity of bubble is enhanced, which would increase heat transfer significantly. At  ms, the difference of the fluid velocity on both sides of the departing bubble is quite evident, which accelerates bubble departing from the heating wall.

Figure 10 shows the velocity curve along the centerline of bubble in different bubble growth times. It shows that, along the centerline of bubble, from inside bubble to the mainstream liquid region around bubble, the velocity increases gradually and at the vapor-liquid interface the velocity fluctuation appears. From 0.2 ms to 1.2 ms, with increase of bubble growth time, the velocity inside the bubble gradually reduced.

4.2. The Dynamic Contact Angle Bubble Shape Evolution

Figure 11 is the evolution of bubble shape from bubble growth to movement on different conditions of case to case . It can be seen that the bulk velocity, liquid subcooling, wall superheat, and vapor-liquid contact angle will lead to different bubble dynamic. Figure 12 shows the dynamic contact angle and the incline angle changes with different working conditions. It can be seen that the dynamic contact angles and are keeping constant in bubble growth; meanwhile is equal to ; . For example in case , the dynamic contact angles and are 45° in bubble growth; in case , the dynamic contact angles and are 60° in bubble growth; in case the dynamic contact angles and are 60° in bubble growth; in case the dynamic contact angles and are 63° in bubble growth. According to contact angle hysteresis theory, the dynamic contact angle values varied between the receding contact angle and the advancing angle on the heating wall in the heterogeneous boiling system. The numerical results of dynamic contact angle in bubble growth results agreed with this theory. Based on the numerical results, the incline angles are all 0° in case to case . It means that the incline angle is little influenced by working conditions on bubble growth process.

When bubble was on the stage of detachment or sliding on the wall, bubble shape and the dynamic contact angle were strongly influenced by working condition. As shown in Figure 12(a), bubble changes slowly from the ellipsoidal shape on the heating wall to spherical shape; the contact area at bottom of bubble on the heating wall slowly decreases during bubble departure progress. As shown in Figure 12(a), during the bubble departure from heating wall on case , the dynamic contact angles and gradually reduced compared with bubble growth. In a certain moment and keep constant and was equal to . The characteristic of dynamic contact angles and on bubble departure and sliding progress still followed the theory of contact angle hysteresis phenomenon; the numerical results of dynamic contact angle on the progress of bubble growth or detachment or sliding from the heating wall varied between receding contact angle and the advance contact angle which was defined as boundary conditions. The incline angle kept constant during bubble determent: °. As shown in Figure 12(b), the characteristic of dynamic contact angle and incline angle in case was the same as in case . As shown in Figure 12(c), in case bubble shape was ellipsoidal shape at the beginning of bubble sliding on the heating wall and persistently changed during bubble sliding on the heating wall. The dynamic contact angles and changed rapidly at the beginning of bubble sliding and kept constant when bubble sliding velocity became stable, while was not equal to . The incline angle varied from 0° to −40° at the beginning of bubble sliding on the wall. As shown in Figure 12(d) the characteristic of dynamic contact angle and incline angle in case was the same as in case while the incline angle varied from 0° to −40° at the beginning of bubble sliding on the wall.

5. Conclusions

In this paper, a numerical investigation on the bubble growth in subcooled flow boiling has been proposed, and several conclusions are obtained as follows.(1)A mechanism model for bubble growth in subcooled flow boiling has been proposed in this study, and the bubble growth curve is obtained based on CFD simulation; the simulation results agree with the experimental results with uncertainty within ±25%.(2)The evolution of bubble shape is complicated. At the initial time of bubble growth, a large amount of heat and mass transfers through the micro-layer into the bubble and the bubble grows fast and expands rapidly as a hemisphere. As the condensation gradually enhances at the top of bubble, bubble growth rate becomes much slower than that of before. The forces acting on a growing bubble have significant effect on the bubble evolution.(3)There are eddies around the growing bubble. The eddy accelerates contraction of the vapor-liquid interface at the bottom of the bubble, which leads to bubble shape changing and bubble detachment, and this will enhance the of micro-convection heat transfer in the vicinity of a bubble.

Nomenclature

:Energy
:Surface tension force
:Latent heat, kJ/kg
:The net mass flow rate per unit area, kg/(m2·s)
:Mass per mole, kg/mol
:Pressure, Pa
:Heat flux, kW/m2
:Bubble radius, m
:Universal gas constant, J/(mol·K)
:The mass source term
:Time, s
:Temperature, K
:Velocity, m/s
:Coordinate axis
:Micro-layer thickness under the bubble bottom, m
:Surface tension coefficient, N/m
: Coefficient of the vapor-liquid interface
:Volume fraction
:Density, kg/m3
:Thermal conductivity of liquid, W/(m·K)
:The temperature difference, K.
Subscripts
:Liquid
:Vapor
:Wall
:Saturated
:Evaporation
:Condensation
:th fluid component.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors are grateful for the support of the National Science Foundation of China (no. 51106142 and no. 51206199) and the project funded by China Postdoctoral Science Foundation (no. 2014M562337).