Numerical study of bubble growth and merger characteristics during nucleate boiling

Progress in Nuclear Energy 112 (2019) 7–19

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

Numerical study of bubble growth and merger characteristics during nucleate boiling

T

Jingliang Bia,∗∗, David M. Christopherb, Dawei Zhaoa, Jianjun Xua,∗, Yanping Huanga a b

CNNC Key Laboratory on Nuclear Reactor Thermal Hydraulics Technology, Nuclear Power Institute of China, Chengdu, 610041, China Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Thermal Engineering, Tsinghua University, Beijing, 100084, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Bubble departure Bubble merger Heat flux variations Nucleate boiling

This paper numerically investigates bubble nucleation, growth, merger, and departure characteristics as well as the heat transfer distributions on the surface, and compared the simulation results with experimental results. Micro-heaters are utilized in experiments to provide constant temperature as well as measure the time-resolved heat flux variations. Bubble images and corresponding heat flux distributions under bubbles are obtained in experiments. A 3-dimensional computational domain with two symmetry planes is built and 5 heating elements are arranged at the bottom surface of the domain to simulate the real heaters used in experiments. The temperatures of these elements are set constant and the local heat fluxes were calculated. Liquid-vapor interface is captured with a volume of fluid (VOF) method in numerical simulations. Bubble growth and merger characteristics and heat flux distributions under bubbles calculated in numerical simulations are analyzed and compared to experimental observations. The heat transfer contributions of different heat transfer modes are calculated. The heat transfer mechanisms during nucleate boiling are better interpreted by the coupling of bubble dynamics and time- and space-resolved heat flux variations on the surface.

1. Introduction Bubble generation and merger is a common phenomenon in many industrial applications including chemical engineering, biochemistry, petrochemical, nuclear power, food industry et al. (Gu and Guo, 2016; Xu et al., 2014; Kuang et al., 2015). In nuclear reactors, the steam bubble dynamics are very important for the reactor safety and control analysis. Due to the various complex phenomena such as interface changes and contact line movement involved during boiling, the mechanisms of nucleate boiling heat transfer are still not fully comprehended. The heat and mass transfer characteristics at liquid-vapor interface and contact line area are important for the investigation of heat transfer mechanisms during boiling. Researchers summarize that microlayer model, transient conduction model and micro-convection model are three most important heat transfer models during nucleate boiling (Rosenhow, 1952; Han and Griffith, 1965; Judd and Hwang, 1976; Kim, 2009). There are numerous studies on nucleate boiling heat transfer and most studies focus on single bubble growth dynamics. Zhang and Shoji (2001) experimentally investigated bubble behaviors during bubble formation and departing. The results reveal that bubble departing

∗

periods are influenced by interactions between successive bubbles and incorporation of the wake effect of the previous bubble. Zhao et al. (2017) carried out a visualization study on micro-bubbles and found that bubbles experienced dramatic deformation before collapse. Due to the complexities of nucleate boiling, numerical investigations of bubble shapes and movements and the related heat transfer are fewer compared to experimental studies. Micro-region model was adopted by some researchers to simulate the heat transfer and vapor evaporation in the microlayer during single bubble growth. Sato and Niceno (2015) developed a new microlayer model and applied it to CFD (Computational Fluid Dynamics) simulation. Single bubble growth characteristics were simulated and the temperature field beneath the growing bubble was measured. The numerical results were compared with previous pool boiling experimental results with water and the comparison results were good. Wu and Dhir (2010) simulated single bubble growth on two-dimensional scale with level-set function and moving mesh method. The computational domain was divided into micro region and macro region. They observed two saddle points in the isotherms at high subcooling condition and explained the phenomena as the condensed liquid along the interface flew downward towards the wall and thins down the thermal layer near the bubble at high

Corresponding author. Corresponding author. E-mail addresses: [email protected] (J. Bi), [email protected] (J. Xu).

∗∗

https://doi.org/10.1016/j.pnucene.2018.12.001 Received 24 December 2017; Received in revised form 5 November 2018; Accepted 2 December 2018 0149-1970/ © 2018 Published by Elsevier Ltd.

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Nomenclature

A Cp D F F¯σ hev hfg k keff L ¯ M m ˙ m mv P q" Qlatent Qinput Q˙ Rb Rh r S Sm t tb T Tcell → u

V x, y, z

Heater area specific heat at constant pressure Bubble diameter force surface tension force evaporation coefficient latent heat of vaporization thermal conductivity effective thermal conductivity characteristic length of each element gas molecular mass mass The unit mass transfer rate in each cell at the interface total mass of the bubble before departure pressure heat flux latent heat required for liquid vaporization total heat supplied by the working heaters heat source term in unit volume Radius of the bubble heater resistance radial coordinate Source term the mass transport among phases time bubble life time temperature fluid temperature in each element velocity

voltage Cartesion coordinates

Greek symbols

Δt α δ κ μ ρ σ σˆ

the current time step size volume fraction microlayer thickness interface curvature viscosity density surface tension accommodation coefficient

Subscripts 0 b d c h int l m mc me q v ref sat tc

initial bubble departure Contact line heater interface liquid mass microconvection microlayer evaporation qth fluid, liquid or vapor vapor reference saturation transient conduction

In recent years, some researchers started to correlate the detailed heat transfer characteristics with bubble dynamics with the aid of high technology in MEMS (Micro-electro-mechanical Systems) and had a deeper understanding of nucleate boiling heat transfer mechanisms. Moghaddam and Kiger (2009) used a micro-sensor array to measure the temperature distributions under single bubbles during one ebullition cycle and calculated the contributions of different heat transfer mechanisms. Gerardi et al. (2010) investigated bubble growth through the synchronization of high-speed video and infrared thermometry from the bottom surface. Bubble merger is an important phenomenon in nucleate boiling and needs to be investigated in detail. Bi et al. has conducted experiments with a microheater array to study the mechanisms involved in the bubble coalescence processes by measuring the heat flux distributions under bubbles and recording bubble dynamics from bottom view (Bi et al., 2013, 2014). The influence of contact line on nucleate boiling has been investigated by some researchers through numerical simulations. Jia et al. (2015) conducted 2-dimensional numerical investigation based on VOF method on single bubble growth and departure at a constant temperature surface at saturated condition. The local heat flux profile in the vicinity of the contact line was plotted. They proposed that the peak heat flux moved with contact line where microlayer evaporation occurred. Mukherjee and Kandlikar (2007) numerically investigated the effect of dynamic contact angle on single bubble dynamics and vapor volume growth rate using level-set method. Kunkelmann et al. (2012) investigated the effect of three-phase contact line speed on evaporative heat transfer using VOF solver of the OpenFOAM CFD package. Literature review shows that numerical research on nucleate boiling mostly focuses on single bubble growth process. However bubble merger and departure processes and the coupling heat transfer characteristics are not intensively studied. The detailed heat flux distributions under bubbles and the correlations with bubble dynamics need to

subcooling condition. Yuan et al. (2016) investigated the bubble growth and sliding process of subcooled flow boiling in narrow rectangular channel using VOF model. The results reveal that a large amount of heat is transferred through microlayer at the initial period during bubble growth. Eddies around the growing bubble accelerates contraction of the vapor-liquid interface at the bottom of the bubble. Deen and Kuipers (2013) presented a direct numerical simulation in dispersed gas-liquid two-phase flow using VOF approach. They investigated the heat transfer rate between the liquid and a hot wall kept at a fixed temperature. The results show that when a bubble departs from the surface, the heat exchange rate between the liquid and the wall is considerably enhanced, which is due to bubble induced agitation in the liquid. Bubble merger processes in both vertical direction and horizontal direction were studied numerically by some researchers. Owoeye and Schubring (2016) conducted a numerical analysis of bubble coalescence in an upward subcooled flow boiling and a microlayer model was adopted to perform near-wall treatment. The bubble shape and rise velocity were compared with experimental results. The bubbles coalesce slower when bubble contact angle increases, while the bubbles coalesce faster when the trailing bubble is larger than the leading bubble. Son et al. (2002) numerically investigated the bubble merger process in vertical direction at one nucleation site with level set method. The vapor removal rate increases almost linearly with the wall superheat regardless of bubble merger process. The total wall heat flux and bubble departure diameter are not much influenced by the vertical merger process. Mukherjee and Dhir (2004) investigated lateral bubble merger phenomena with level set functions and compared the simulated bubble dynamics with experimental observations. The temperature fields around the bubble were obtained. The researchers found that bubble merger can increase the vapor removal rate and the merger of multiple bubbles increases the total wall heat transfer significantly. 8

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

u) ∂ (ρ→ T + ∇⋅(ρ→→ u u ) = −∇ p+ ∇⋅(μ(∇→ u + ∇→ u )) + ρ→ g + F¯σ ∂t

be investigated in depth. Three-dimensional numerical simulations on bubble growth, merger and departure processes compared with experimental results will be conducted in this paper. The detailed heat flux distributions under bubbles and the correlating bubble dynamics will be obtained and compared with our own experimental results obtained with a micro-heater array.

The Boussinesq approximation is used for the buoyancy term. A single set of momentum equations is used for both phases with appropriate changes in the thermophysical properties. The source term is surface tension force acting on vapor-liquid interface and contact line on the wall. The surface tension force term is modeled with CSF (Continuum Surface Force) model. The surface tension force is converted into the following volume force,

2. Numerical model 2.1. Computational domain

F¯σ = σ

A 3 mm × 2 mm × 2 mm computational domain is set up in Fluent 6.3 to model the bubble growth and merger process as seen in Fig. 1. Since the nucleation sites are arranged symmetrically, X = 0 and Y = 0 planes are specified as symmetry planes, and only 1/4 of the total domain, which is 1.5 mm × 1 mm × 2 mm, is specified as the real computational domain. Thus, only one bubble grows in the computational domain, but it combines with its mirror image when it reaches the X = 0 plane. The liquid is set as FC-72, and the contact angle is 53°. The computational domain includes five 0.1 mm × 0.2 mm heaters, and the temperature of each heater can be controlled as in experiments and the local heat flux can be calculated. Heater 1 is most close to Y-Z symmetry plane, and heaters 2, 3, 4, 5 are farther. By controlling the surface temperature, bubble nucleation is controlled.

αl ρl κ v ∇α v + α v ρv κl ∇αl 0.5(ρl + ρv )

(3)

where κ is the curvature of the interface, and σ is the surface tension coefficient. Energy Equation is:

∂ (ρE) + ∇⋅(→ u (ρE + p)) = ∇⋅(k eff ∇T) + Q˙ ∂t

(4)

Where the energy is a mass averaged energy based on the vapor and liquid fractions in each cell, ρ and keff (effective thermal conductivity) are shared among the phases, Q˙ is heat source term in unit volume transporting through interface due to vaporization and condensation. The VOF model considers energy and temperature, as mass averaged variables, and in this case, there are two phases, and T is expressed as:

2.2. Meshing method There are 90 × 40 × 60 computational cells in X, Y, and Z directions of the computational domain. The whole computational mesh includes 216,000 hexahedral cells. The cell size near the heater is very small and the size of the cell adjacent to the heater is only 0.015 mm. Five 3-μm-high boundary layers are set on the surface for better calculation near the surface. 120 × 60 × 90 computational grids are also tested in the calculation. The results are similar to the results using this mesh, so 90 × 40 × 60 grids are used in the research to save computing time and space. A time step size of 1 × 10−6 s is utilized in the simulations for better convergence.

2

T=

∑q = 1 αq ρq Tq 2

∑q = 1 αq ρq

(5)

The above governing equations are solved using finite volume method and the advection terms are solved with QUICK method in Fluent 6.3. The Geo-reconstruct scheme is selected for VOF model. The PISO method is used for the pressure-velocity coupling equation and the pressure terms are solved with Presto! method.

2.3. Governing equations The analysis is based on transient, two-phase, three-dimensional laminar flow with the Boussinesq approximation for the buoyancy term. When the liquid superheat on heaters reaches a certain value, the superheated liquid starts to vaporize. Basic governing equations describing momentum, mass and energy transport are solved for an equivalent fluid whose physical properties are calculated as functions of the physical properties of its constituent phases and their volume fractions. The bubble dynamics were modeled based on VOF with separate liquid and vapor phases. The VOF model (Hirt and Nichols, 1981), assumes that all immiscible fluid phases presenting a control volume share the same velocity, pressure and temperature fields. In VOF model, the volume fraction of liquid αq is introduced. When the cell is filled with fluid, i.e. when αq = 0, there is no fluid in the cell. When αq is between 0 and 1, there is a mixture of liquid and vapor in the cell, i.e. the interface. In each control volume, the volume fractions of two phases sum to unity. Continuity Equation is:

∂ (αq ρq ) + ∇⋅(αq ρq → u ) = Sm ∂t

(2)

(1)

q represents liquid or vapor, source term Sm represents the mass transport among phases. The mass transport is due to liquid vaporization and vapor condensation. Momentum equation is:

Fig. 1. The computational domain. 9

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

¯ is Where hfg is liquid latent heat, pv is pressure inside the bubble, and M gas molecular mass. The ideal vaporization heat transfer coefficient is 2.48 × 105 W/m2, however some studies show that the evaporation and condensation accommodation coefficient σˆ may be much less than 1 especially for water (Carey, 1992). σˆ is assumed to be 0.1 for FC-72, indicating the non-ideal effect of non-condensable gas at the interface. The vapor temperature is equal to the saturation temperature for the pressure inside a bubble of the given size based on the ClausiusClapeyron equation

2.4. Initial conditions Bubble nucleation requires a certain superheat of adjacent liquid as can be observed from experiments. So the heaters were heated for 10 ms before two-phase model was applied. The time step size at natural convection stage is 1 × 10−4 s, larger than the time step in twophase calculations since there is no complex phase change in the preheating period. 2.5. Boundary conditions

Tv = Ts + Working heaters are set at a constant temperature of 373 K, and the other part of the bottom surface is set as adiabatic wall. Pressure outlet boundary is applied at the top and outer surfaces of the computational domain and the temperatures are set at 304 K. The system pressure is kept at 38 kPa and the corresponding FC-72 saturation temperature is 304 K. The heat flux variations with time of all heaters were exported to a file.

˙ ≤ Vaporization: m

Se = (6)

⎜

(10)

(Tsat − T ) hc (ρl cP, l − ρv cP, v )(Tsat − Tref ) hfg L ρl

(11)

where the second fraction is the enthalpy change in both phases relative to the liquid mass transfer given in Eq. (10). The momentum source terms were then:

→ → ˙ Sm = mU

(12)

Those computed equations including the liquid source term, energy source term, and momentum source term were implemented in Fluent as UDF subroutines. 3. Heat transfer models during nucleate boiling In nucleate boiling, heat transfer mainly includes three parts, microlayer evaporation, transient conduction, and microconvection. Those three mechanisms were illustrated in Fig. 2. The heat transfer due to these three mechanisms was analyzed below. 3.1. Microlayer evaporation

1

2 ¯ 2 2σˆ hfg ρv ⎛ M ⎞ ⎡1 − pv ⎤ ¯ v⎠ ⎢ 2 − σˆ Tv ⎝ 2πRT 2hfg ρv ⎥ ⎣ ⎦

ρv αl Δt

(9)

where ρv is vapor density, α v is vapor volume fraction, αl is liquid volume fraction and Δt is the current time step size. The corresponding energy source term was then calculated as:

Where L is the characteristic length of each element which is assumed to be the cubic root of the element volume, and Tcell is two-phase fluid temperature in each element. Bubble nucleation requires a certain superheat as observed in previous experiments, and the superheat is found to be 40 K according to previous experimental results. In UDF, before the nucleation occurs, when the liquid superheat is larger than 40 K, liquid vaporizes, and minus sign is applied in Eq. (6). When there is already vapor on the surface, liquid vaporizes as long as the liquid superheat is large than the saturation temperature. When the vapor temperature is lower than saturation temperature, vapor condenses to liquid, and plus sign is applied in Eq. (6). The source term in Eq. (6) was then used as a liquid source term in Fluent with positive values representing condensation and negative values representing evaporation. The boundary condition at the bubble interface used to model the evaporation and condensation at the interface is based on a convective heat transfer coefficient according to Carey:

hev =

ρv α v Δt

˙ ≤ Condensation: m

The source terms including the heat, momentum and mass transfer rates at the vapor-liquid interface are defined in User Defined Functions (UDF). The unit mass transfer rate in each cell at the interface is defined as,

(Tsat − Tcell ) hev hfg L

(8)

To avoid unreasonable results, the unit condensation in each cell during each time step is limited not to exceed the total vapor amount in the cell, and the maximum vaporization not to exceed the total liquid amount in the cell.

2.6. Modeling of source terms

m˙ = ±

2σTs hfg ρv Rb

⎟

(7)

Microlayer exists beneath a growing bubble, and this thin layer

Fig. 2. Bubble growth and departure illustrations. 10

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

flux during the rewetting process can be expressed as:

evaporates during bubble growth. The evaporation of microlayer generates the vapor for bubble growth. However, since this thin layer is too small, it is difficult to directly observe and measure the microlayer. The microlayer is considered within a cylindrical coordinate system. The continuity equation for the microlayer can be expressed as:

∂vf 1 ∂ (ruf ) + =0 r ∂r ∂y

q"tc =

(13)

q (tR) =

(14)

In the dryout region, Tint = Tw, κ = 0 and q" = 0 . Substituting for the heat flux in the above equation, the thickness of the dryout area can be expressed as:

δ03

=

∫R

R c

q"tc

2πrdr = tR − tr

∫R

R c

kf ΔT 2πrdr πα tR − tr

(19)

Assuming the bubble contact line receding velocity, uc , is constant, and then the time when the contact line approaches R will be tR = (R c − R)/ uc , and the total rewetting time will be tre = R c / uc . Substituting tR and tre into Eq. (19), and the heat transfer rate at the moment tR is:

2

q″ A + δ3 ρv hfg2

(18)

When the liquid front approaches to a position where the radius is R from Rc (0 ≤ R ≤ R c ), the heat transfer rate at time tR in this rewetting process can be obtained by integrating the heat flux from Rc to R where the liquid front is encountered, resulting in

Pressure drop across the interface is calculated through the modified Young-Laplace equation which is obtained by a force balance which includes the surface tension, disjoining pressure, and the vapor recoil force:

Pf = Pv − σκ −

kl (Tw − Tl ) παt

q (tR) =

4πkf ΔTuc πα

2 ⎛R c tR1/2 − uc tR3/2 ⎞ 3 ⎝ ⎠

(20)

A

(

Tint Tv

)

− 1 hfg ρv

3.3. Microconvection and natural convection

(15)

The slope of the microlayer at Rmax is tan(θ) and the radius of the outer edge of the microlayer is R sin θ hence the boundary conditions at the microlayer can be presented as:

δ (R 0) = δ0,

During nucleate boiling, the rapid bubble growth perturbs the liquid and enhances the natural convection around the bubble. In addition to the bubble growth perturbing the liquid, bubbles can also oscillate slightly on the heating surface during bubble growth. The natural convection is enhanced due to bubble growth, oscillation and departure. This part of enhanced convection around the bubble is referred to as micro-convection. The local Nusselt number for forced convection over a flat plate is used to calculate micro-convection. According to Forster and Zuber, the bubble Reynolds number was defined as:

δ′ (R 0 ) = 0

δ′ (Rmax ) = tan(θ), Rmax =

D sin(θ) 2

(16)

The heat flux through the microlayer is based on Fourier's law,

q″me = kf

(Tw − Tint ) δ

(17)

Reb = 3.2. Transient conduction during bubble departure

ρv Ub D μf

(21)

Where the equivalent diameter of bubble is defined as the diameter of a sphere with the same volume as the bubble, and Ub = dRb / dt . The Nusselt number due to microconvection is expressed based on laminar flow over a flat plate by

Transient conduction occurs when the bubble is ready for departure and the contact line starts to shrink and the surface under the bubble is rewetted. The heat flux increases sharply during bubble departure, which has been observed in previous experiments. The bubble radius is R c when it starts to recede, and R c is the radius of the contact line. The rewetting process lasts until the bubble base radius decreases to zero as the bubble departs. According to Mikic and Rohsenow, the surface heat

1

Numc = C1 Rem b (Pr)

3

(22)

where C1 is 0.332 as is the case for forced convection over a flat plate, and m is 0.67 as in Rohsenow's model. The heat flux due to

Fig. 3. Experimental setup schematic. 11

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

2000 Hz per channel. The bubble dynamics were recorded by a highspeed CCD camera from the side visualization windows at the same time. The total heat flux for each heater was then calculated from the voltage across the heater and calibrated resistances:

microconvection is calculated as: 1

q"mc =

0.332 Reb0.67 (Pr) 3⋅kf 4f ⋅Dd

(23)

4. Experimental setup and results

q″ = (V 2/ Rh)/ A

Experiments were also conducted to investigate the coupling effects of bubble dynamics and heat transfer characteristics. Bubbles were generated on a micro-heater array which can record the heat fluxes at the same time. The synchronization study of bubble shapes and the heat flux distributions under bubbles was conducted.

where V is the voltage across the heater and A represents the heater area. The heat losses due to substrate conduction and heater conduction were deducted from the total heat flux. The heater resistances were calibrated by measuring the resistances at different temperatures. The detailed heat flux data processing and heater calibration process can be seen in previous published paper by Bi et al. (2013). The bubble size measured with the image analysis software had an uncertainty of about 2%. The uncertainty in the voltage measurements was less than 2% due to uncertainties in the measurement equipment. The heater resistances were measured with an accuracy of better than 0.5%. The overall boiling heat flux calculations then had an uncertainty of less than 4.8%.

4.1. Experimental apparatus The schematic of the experimental system in this study is shown in Fig. 3. The boiling rig is an 80 mm high transparent quartz glass vessel so bubble dynamics can be observed from side view. FC-72 is used as the boiling fluid because it is dielectric and will not damage the heaters. The saturation temperature of FC-72 is 304 K for the system pressure of 38 kPa. The system pressure was controlled by a vacuum pump and was measured by a pressure transducer through an access port on the top of the boiling rig. The bulk temperature was measured by a K-type thermocouple placed in the bulk liquid FC-72 at the boiling rig, and was recorded by a high-speed Agilent data acquisition system. The heater temperature can be kept constant through a Wheatstone bridge circuit. The heater temperature was controlled at different values by setting the resistance of variable resistor through a C++ program. The detailed heater control process can be seen in previous publication by Bi et al. (2013). Voltages across the heaters were measured by a high speed data acquisition system with a sampling rate of

(24)

4.2. Experimental results The bubble merger dynamics from bottom view are shown in Fig. 4. Two bubbles grow large enough to contact each other. When two interfaces contact, the interfaces merge and the two bubbles start to merge into one. The merging bubble shrinks at axial direction, while elongates at the perpendicular direction as shown from 1.5 ms to 7.0 ms in Fig. 4. The bubble surface was still irregular during this process due to the rapid motion. Then the bubble stretches at axial direction and contracts at perpendicular direction as shown from 7.0 ms to 11.0 ms. The surface tension forces and the pressure in the bubble have caused

Fig. 4. Bubble merger dynamics from bottom view in experiments. 12

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Fig. 5. Experimental heat fluxes for (a) Heater 5 (b) Heater 4 during bubble coalescence events.

Fig. 6. Comparison of numerical bubble growth dynamics with experimental results (a) numerical simulations (b) experimental observations.

process, the heaters on the edge of the bubble, for example heater 5, will be wetted and small bubbles generate on the surface. After 11 ms, the bubble necking process begins and the contact line is shrinking. More heaters are exposed to cold liquid and transient conduction heat

the interface motion in the axial direction to reverse direction and begin moving back inwards just as the interface motion in the perpendicular direction is reversing direction and moving outwards. The bubble surface is more regular in this process. During the bubble deforming 13

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Fig. 7. Velocity vectors at y = 0 plane during bubble growth and necking stages (a) t = 1.1 ms, (b) t = 16.6 ms.

at 373 K and the other heaters are set as adiabatic wall. Due to the small time step size in numerical simulations, the bubble interface variations can be observed during a short time interval of 1 × 10−6 s. Hence the simulated bubble dynamics are more detailed than the images obtained from experiments. Fig. 6 compares the numerical bubble growth characteristics in one ebullition cycle with experimental observations. The bubble grows rapidly after nucleation, especially for the first 1 ms. The bubble shape is long in horizontal direction for the first 8 ms after nucleation, and then the bubble stretches in the vertical direction and shrinks in the horizontal direction. The bubble starts necking from 13.2 ms and the contact line shrinks sharply during the necking process, and then the bubble soon departs at 17 ms. Fig. 7 shows the velocity field around and inside the bubble at 1.1 ms and 16.6 ms, representing growth and necking stages. At 1.1 ms, the bubble starts to grow and the growth rate is relatively large, and the velocity inside the bubble can reach 0.54 m/s. The vapor-liquid interface moves fast and pushes the liquid in the vicinity of the fast growing bubble. There are eddies around the growing bubble as shown in Fig. 7(a), and these eddies and the movement of the adjacent liquid enhances the micro-convection. At 16.6 ms, the vapor velocity in the bubble neck is rather large, up to 3.59 m/s, much larger than the vapor velocity during growth, indicating a very large tendency to depart from the surface. During the necking period, cold liquid flows to the heating surface and pushes hot liquid away. Transient conduction takes place and the heat transfer is strongly enhanced.

Fig. 8. Comparison of numerical and experimental heat fluxes during single bubble ebullition cycle.

transfer is very high in the wetted region. Heater 4 is wetted during the necking process. The bubble departs from the surface at 15 ms. The experimental heat fluxes for different heaters were shown in Fig. 5. During bubble growth period, the heat flux was about 10 W/cm2. After two bubbles coalesce, two bubbles merge into one and the merged bubble shrinks and deforms on the surface for a long time as shown in Fig. 4. During this deformation process, the heating surface is rewetted repeatedly and the heat transfer outside of the bubble contact line is greatly enhanced, so heat flux spikes occur for the wetted heaters after bubble coalescence. The heat flux spikes for heater 5 last longer and were higher than the heat flux spikes for heater 4. This is because heater 5 is rewetted more than heater 4 during the bubble deforming process after bubble merger.

5.1.2. Heat flux characteristics during single bubble growth and departure The heat fluxes calculated in numerical simulations are compared with the heat fluxes observed in experiments. The comparison results are shown in Fig. 8, and the numerical and experimental plots show similar trends. The heat flux for heater 2 increases during the whole growth period due to the shrinking of the bubble base line in the horizontal direction. According to numerical results, since 13.2 ms, the bubble starts to neck, and the heat fluxes of heater 2 start to increase sharply. After 0.1 ms, the heat flux for heater 1 also increases. Then two small bubbles start to nucleate on the two sides of the large necking bubble. At 15.2 ms, the heat flux for heater 2 reaches 89 W/cm2. Then the heat flux for heater 1 reaches 40 W/cm2 at 16.1 ms. The heat fluxes have some oscillations during the bubble necking period. At 17 ms, the bubble departs from the surface and vapor generates on the surface. The heat fluxes for heater 1 and 2 drop quickly, as seen from the numerical plots. The comparison results show that the heat fluxes are relatively low for most of the ebullition time in both numerical simulations and experimental results, as seen in Fig. 8 from 0 ms to 13 ms. The explanation for the low heat fluxes is that during the whole single bubble ebullition cycle, heater 1 is in the center area of the bubble and the heat

5. Results and discussions 5.1. Single bubble dynamics and heat transfer characteristics Single bubble growth and departure characteristics are numerically investigated in this paper. Heater 1 and 2 are set at 373 K and other heaters are set as adiabatic conditions. The boundary conditions in the numerical model are the same with experimental boundary conditions. 5.1.1. Single bubble growth and departure In numerical simulations, the liquid saturation temperature for the system pressure is 304 K, and the bulk liquid temperature is set at 304 K as in experiments, so the bulk liquid is saturated. Heater 1 and 2 are set 14

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Fig. 9. Comparison of bubble merger simulations and experimental results (a) numerical simulations (b) experimental observations.

conduction from the heater to vapor is small. The heat flux increases sharply for 4 ms before the bubble departs from the surface. Heater 2 is on the edge of the bubble during bubble growth and departure period. The heat fluxes for heater 2 during bubble growth period increase with time in both numerical calculations and experimental results, which is due to the shrinking of the bubble contact area in the horizontal direction during bubble growth. The increasing rate is small before 13.2 ms, and then the heat flux starts to increase sharply. There are some differences between the numerical and experimental heat flux plots because the contact line movement is not exactly the same in numerical simulations and in experiments.

5.2. Bubble merger characteristics 5.2.1. Bubble dynamics during merger events In numerical simulations, the bulk liquid is saturated and the temperature is 304 K. VOF model was applied in the simulations to track the bubble interface during merger process. When two bubble interfaces reach to one same computational cell, or the spacing between adjacent interfaces falls below the grid spacing, the two interfaces merge into one. Heater 1 is set as adiabatic wall, and the other heaters are set at 373 K. Fig. 9 displays the comparisons of bubble merger and departure

15

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Fig. 10. Velocity vectors at y = 0 section planes during bubble merger events (a) 0.2 ms, (b) 1.2 ms, (c) 4.9 ms, (d) 11.8 ms.

Two bubbles have been completely merged at 1.2 ms. Fig. 10(b) shows the velocity vectors when two bubbles start shrinking to the middle along X axis. The velocity of the vapor inside the bubble is upwards along the positive direction of Z axis and towards the center axis of the bubble, which explains the bubble neck shrinking along X axis and elongating along Z axis. At 4.9 ms, two small bubbles on heaters 5 at two sides are being absorbed by the large deforming bubble. As seen from Fig. 10, the vapor velocity inside the small bubbles is very high and reaches 1 m/s. The vapor moves towards the large bubble, indicating the directions of bubble movement. At 11.8 ms, the merged bubble soon starts necking. The vapor velocity in the bubble neck is rather large and is about 1–2 m/s. The vapor moves towards the axis of the large bubble, i.e. the Z axis, indicating the rapid shrinking of the large bubble neck. The vapor in the bubble also moves in positive direction of the Z axis, since the bubble is about to depart from the surface.

Fig. 11. Numerical simulations of heater heat fluxes during bubble merger events.

process between the numerical simulations and experimental results. Fig. 9 (a) shows the simulated bubble contours at y = 0 section plane and the picture size is 3 mm × 2 mm. Fig. 9 (b) shows the side view of bubble images obtained from experiments, and the maximum width of the picture is 1.66 mm. As is shown in Fig. 9 (a), two bubbles coalesce at 0.1 ms, and then the two bubble interfaces start merging. The merged bubble starts shrinking to the middle in X axis direction, and expands in the Y and Z direction. Since 4.7 ms, the bubble stops shrinking and starts expanding slowly from the bottom, and finally forms an ellipsoid which is very close to a spherical. The bubble starts necking since 10.4 ms, and detaches at 15.4 ms. The numerical results are compared with experimental pictures, and the merged bubble shapes are very close to what are observed from experiments. Fig. 10 displays the velocity vectors at y = 0 plane at various moments. The bold green lines represent the two-phase interface. Fig. 10(a) shows the velocity contours of two phases when the two bubbles start to merge at 0.2 ms. As the interface is merging, the liquid above the merged interface is pushed up towards the positive direction of Z axis by the deforming bubble. The fluid in the trapped liquid layer flows at about 0.4 m/s along the heating surface. The flow enhances the heat transfer of the wall under the trapped liquid layer and results in a small heat flux spike.

5.2.2. Heat transfer characteristics during bubble merger events The heat fluxes of different heaters calculated with numerical methods are also compared with experimental measurements, as are shown in Fig. 11. The trends of the heat flux variations are similar; however the heat flux values are not exactly the same. This is probably because the contact line movement is not exactly the same in numerical simulations and experiments, especially during merger events. In experiments, the bubble is more easily disturbed by the liquid flow around the heating area. The merged bubble moves violently, and the contact line shapes are irregular and the contact line variations are even larger than in the necking phase. Fig. 12 displays the contact line movement on heaters. The corresponding bubble dynamics and contact line variations are also analyzed to investigate their effects on heat flux distributions under bubbles. The controlling heat transfer mechanisms during bubble growth and merger process are hence revealed. The bubble interfaces are merging into one from 0.1 ms to 1.2 ms. From 1.2 ms, the merging bubble starts to shrink to the middle when heater 5 is rewetted as seen in Fig. 12(a), and the heat flux increases sharply at this moment. The merged bubble keeps shrinking to the middle along X axis direction and heater 5 is completely rewetted at 3.2 ms. The heat flux increases from 10.4 W/cm2 at 1.2 ms to 50–70 W/ 16

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

Fig. 12. Contact line variations during bubble merger (a) 1.2 ms, (b) 2.3 ms, (c) 2.9 ms, (d) 4.9 ms, (e) 11.8 ms, (f) 15.4 ms.

more areas of heaters 2 and 3 are wetted, as seen in Fig. 12(e). The heat fluxes increase sharply during bubble necking. The heat fluxes for all heaters drop after bubble detaches from the surface since vapor generates on the surface and the heat conduction from surface to vapor is small. Heat flux variations for heater 4 and heater 5 measured from experiments are compared with numerical results in Fig. 11. The heat flux is relatively low before bubble coalescence. However, the heat flux starts to increase rapidly from the bubble coalescence moment. The edge heaters have higher heat flux peaks in both numerical simulations and experiments. The heat flux peaks are not exactly the same because the contact line shape and movement is very complex and the contact line can hardly be accurately predicted. The heat flux data indicates that during bubble merger events, the heat flux is low before bubble coalescence, while the oscillation and deformation of the merged bubble results in the rewetting of the heaters. The cooler liquid and hot surface has large heat transfer rate in the form of transient conduction. The moments of heat flux peaks are the same with the heater rewetting moments, which proves that the rewetting is the main reason of the heat flux peaks. Like in single bubble case, the edge heaters have higher and longer heat flux peaks since the edge heaters are more easily exposed to the cooler liquid. Transient conduction due to heater rewetting is dominant heat transfer mode during bubble merger.

Fig. 13. Comparison of single bubble growth plots during one ebullition cycle.

cm2 and remains high. Heater 4 is partially rewetted from 2.9 ms as seen in Fig. 12(c) and heat flux increases slightly. As seen from Fig. 12, a small bubble generates on heater 5 at 2.9 ms, when the heat flux of heater 5 increases sharply from 68.7 W/cm2 to 85.9 W/cm2, indicating that bubble nucleation requires much energy and leads to high heat flux. The new bubble is absorbed by the deforming bubble at 3.4 ms and the heat flux drops to 68 W/cm2. At 3.7 ms, a small bubble generates on heater 5 and the heat flux increases to 84.2 W/cm2, and the small bubble is soon absorbed by the adjacent large bubble. This process repeats for many times and results in many heat flux spikes of heater 5. As small bubbles generate and are absorbed by the large bubble, the merging bubble shrinks to the middle and heater 4 is partially rewetted, resulting in an increase of the heat flux for heater 4 from 10 W/cm2 to 40 W/cm2. The heat flux for heater 4 remains at about 40 W/cm2 since 4.6 ms, and varies slightly with the oscillation of bubble contact line, as seen in Fig. 12(d). Bubble starts necking at 10.4 ms, and the necking becomes fast after 11.9 ms and

5.3. Heat transfer portion during bubble growth and merger process The bubble growth plots from both numerical simulations and experimental observations are shown in Fig. 13. The bubble equivalent diameters are plotted against the growth time. The plot shows that for the first 1 ms, the bubble grows very fast and the bubble growth rate is very large, and then after the bubble diameter reaches 0.36 mm, the bubble growth rate keeps constant. The observed bubble diameters and numerically calculated bubble diameters agree well. The heat transfer mechanism during nucleate boiling is analyzed by solving for the energy required for bubble growth and the total energy provided by the heaters during one bubble life cycle. The liquid absorbs 17

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al.

events, the total heat supply from the powering heaters is 1.42 × 10−3, while the heat required for the net liquid vaporization is 1.3 × 10−4. The ratio is 9.1%, also less than 10%. The calculated results show that for nucleate boiling, the latent heat required for the bubble growing to the observed departure diameter is only a small part of the heat provided by the heater, and the majority of the heat transfer is due to transient conduction and enhanced microconvection during bubble necking, merger, and departure processes. 6. Conclusions Bubble shapes and the corresponding surface heat fluxes are numerically investigated during bubble nucleation, growth, merger, and departure events. The heat, momentum and mass transfer rates at the vapor-liquid interface are defined in the simulation. The synchronized experimental study of bubble dynamics and heat flux distributions is also carried out with micro-heaters and high speed camera. The detailed heat flux distributions with bubble movement can be analyzed. The results show that numerically predicted bubble shapes and heat flux distributions were in good agreement with time-resolved experimental data. During single bubble events, the heat flux in the dryout area is relatively low, while the heat fluxes around the contact area are much higher. Heat flux peaks occur during bubble necking and departure stages in both simulations and experiments. In bubble merger cases, the oscillation and deformation of merged bubble results in the rewetting of edge heaters and constant nucleation of small bubbles which further leads to high heat fluxes of edge heaters. The heat fluxes of central heaters increase only during bubble necking and departing stage. The calculation results show that latent heat required for bubble growth only accounts for less than 10% of the total energy provided by the powering heaters. Most of the heat transfer is due to surface rewetting taking place during bubble necking, deformation and departure process, indicating that transient conduction during bubble merger and departure is the main mechanism of high heat transfer rate during nucleate boiling.

Fig. 14. Comparison of the total heat supplied by heaters and the required energy for bubble growth in one (a) single bubble growth and (b) bubble merger.

Acknowledgments

latent heat and vaporizes during the whole bubble life cycle. The total latent heat the bubble needs to grow to departure diameter is based on the total vaporization mass:

This work was supported by the National Natural Science Foundation of China, Grant No. 11475161, No. 11505177 and No. 51476091.

Qlatent = m v hfg

Appendix A. Supplementary data

(25)

Where Qlatent is the total latent heat the bubble needs to grow for departure, and mv is total mass of the bubble before departure. The heat provided by all the heaters during this bubble life time is:

Qinput =

∑ q"Ah tb

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.pnucene.2018.12.001. References

(26)

where q" is the heat flux of the heater, Ah is the heater area and tb is the bubble lifetime. The ratio of the latent heat to the total heat input is then:

ratio = Qlatent / Qinput

Bi, J.L., Lin, X.P., Christopher, D.M., Li, X.F., 2013. Analysis of coalescence phenomena on microheaters at two surface superheats. Int. J. Heat Mass Tran. 67, 798–809. Bi, J.L., Christopher, D.M., Lin, X.P., Li, X.F., 2014. Effects of nucleation site arrangement and spacing on bubble coalescence characteristics. Exp. Therm. Fluid Sci. 52, 116–127. Bi, J.L., Vafai, K., Christopher, D.M., 2015. Heat transfer characteristics and CHF prediction in nanofluid boiling. Int. J. Heat Mass Tran. 80, 256–265. Carey, V.P., 1992. Liquid Vapor Phase-transition Phenomena. Hemisphere, New York, pp. 112–120. Deen, Niels G., Kuipers, J.A.M., 2013. Direct numerical simulation of wall-to liquid heat transfer in dispersed gas–liquid two-phase flow using a volume of fluid approach. Chem. Eng. Sci. 102, 268–282. Gerardi, C., Buongiorno, J., Hu, L., McKrell, T., 2010. Study of bubble growth in water pool boiling through synchronized, infrared thermometry and high-speed video. Int. J. Heat Mass Tran. 53, 4185–4192. Gu, H.Y., Guo, L.J., 2016. Modeling of bubble shape in horizontal and inclined tubes. Prog. Nucl. Energy 89, 88–101. Han, C.Y., Griffith, P., 1965. The mechanism of heat transfer in nucleate pool boiling. Int. J. Heat Mass Tran. 8, 887–914.

(27)

Fig. 14 shows the required energy for phase change and total heat supplied by the working heaters during bubble growth and merger events. The total heat supplied by the working heaters is 9.1 × 10−4 J, while the latent energy required for bubble nucleation and growth is 7.18 × 10−5 J. The calculated heat input ratio shows that the ratio is only 7.9%, which means that the microlayer evaporation heat accounts for less than 10% of the total energy while transient heat due to bubble sliding accounts for more than 90% of the total heat transfer over a bubble ebullition cycle. The numerical results are consistent with previous theoretical analysis (Bi et al., 2015). During bubble merger 18

Progress in Nuclear Energy 112 (2019) 7–19

J. Bi et al. Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201–225. Jia, H.W., Zhang, P., Fu, X., Jiang, S.C., 2015. A numerical investigation of nucleate boiling at a constant surface temperature. Appl. Therm. Eng. 88, 248–257. Judd, R.L., Hwang, K.S., 1976. A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation. ASME J. Heat Transfer 98, 623–629. Kim, J., 2009. Review of nucleate pool boiling bubble heat transfer mechanisms. Int. J. Multiphas. Flow 35, 1067–1076. Kuang, Y.W., Wang, W., Zhuan, R., Yi, C.C., 2015. Simulation of boiling flow in evaporator of separate type heat pipe with low heat flux. Ann. Nucl. Energy 75, 158–167. Kunkelmann, C., Ibrahem, K., Schweizer, N., Herbert, S., Stephan, P., GambaryanRoisman, T., 2012. The effect of three-phase contact line speed on local evaporative heat transfer: experimental and numerical investigations. Int. J. Heat Mass Tran. 55, 1896–1904. Moghaddam, S., Kiger, K., 2009. Physical mechanisms of heat transfer during single bubble nucleate boiling of FC-72 under saturation conditions1-experimental investigation. Int. J. Heat Mass Tran. 52, 1284–1294. Mukherjee, V.K., Dhir, 2004. Study of lateral merger of vapor bubbles during nucleate pool boiling. J. Heat Tran. 126, 1023–1039. Mukherjee, A., Kandlikar, S.G., 2007. Numerical study of single bubbles with dynamic contact angle during nucleate pool boiling. Int. J. Heat Mass Tran. 50, 127–138.

Owoeye, E.J., Schubring, D.W., 2016. CFD analysis of bubble microlayer and growth in subcooled flow boiling. Nucl. Eng. Des. 304, 151–165. Rosenhow, W.M., 1952. A method of correlating heat transfer data for surface boiling of liquids. Trans. ASME J. Heat Transfer 4, 969–975. Sato, Y., Niceno, B., 2015. A depletable micro-layer model for nucleate pool boiling. J. Comput. Phys. 300, 20–52. Son, G., Ramanujapu, N., Dhir, V.K., 2002. Numerical simulation of bubble merger process on a single nucleation site during pool nucleate boiling. J. Heat Tran. 124, 51–62. Wu, J., Dhir, V.K., 2010. Numerical simulations of the dynamics and heat transfer associated with a single bubble in subcooled pool boiling. J. Heat Tran. 132 111501-1-15. Xu, J.J., Chen, B.D., Xie, T.Z., 2014. Experimental and theoretical analysis of bubble departure behavior in narrow rectangular channel. Prog. Nucl. Energy 77, 1–10. Yuan, D.W., Xiao, Z.J., Chen, D.Q., Huang, Y.P., 2016. Numerical investigation on bubble growth and sliding process of subcooled flow boiling in narrow rectangular channel. In: Science and Technology of Nuclear Installations, Article ID 7253907, 12 pages. Zhang, L., Shoji, M., 2001. Aperiodic bubble formation from a submerged orifice. Chem. Eng. Sci. 56, 5371–5381. Zhao, L., Mo, Z., Sun, L.C., Xie, G., 2017. A visualized study of the motion of individual bubbles in a venturitype bubble generator. Prog. Nucl. Energy 97, 74–89.

19