Bubble Dynamics of a Single Condensing Vapor Bubble from Vertically HeatedWall in Subcooled Pool Boiling System : Experimental Measurements and CFD Simulations

Bubble dynamics of a single condensing vapor bubble in a subcooled pool boiling system with a centrally heated cylindrical tank has been studied in the Rayleigh number range 7.9 × 10 < Ra < 1.88 × 10. The heat source in the system is steam condensing inside a vertical tube. The tube was placed in the center of the tank (300 mm i.d., 300 mm height) which is well filled with water. Experimental investigation has been carried out with High Speed Camera while Computational Fluid Dynamics (CFD) investigation has been performed using Volume of Fluid (VOF) method. The heat source has been modeled using simple heat balance. The rise behavior of condensing bubbles (change in size during rise and path tracking) was studied and the CFD model was validated both quantitatively and qualitatively.


Introduction
Pool boiling is a complex nonlinear dynamic process.The microscopic aspects of boiling process, for example, embryo nucleation and surface roughness determine its nature largely.The microscopic effects eventually exhibit in a macroscopic behavior.Further, boiling is a highly stochastic process and the precise prediction of the location and time of the generation, collapse, coalescence is difficult with the present status of knowledge.The condensation of bubbles in subcooled boiling systems is extremely important for studying the hydrodynamics during subcooled pool boiling.The simple reason for this is that condensation changes the shape and size of the bubbles.When a bubble is formed at the interstitial site in the case of nucleate boiling the bubble grows to a certain diameter and then departs from the interstitial site.As the bubble departs from the nucleation site, it tends to slide along the hot surface and detaches from the surface.The bubble then does not necessarily follow a straight path but moves near or away from the wall.In addition to experimental investigations, Computational Fluid Dynamics (CFD) is useful to obtain better understanding of bubble dynamics.
Experimental investigations on dynamics of single vapor bubbles condensing in a subcooled boiling liquid have been carried out by Voloshko and Vurgaft [1] and an analytical model was developed.The authors found the relation between the dynamics of vapor bubble during condensation and intensity of heat exchange at vapor bubble interface.Bubble dynamics of vapor bubbles in subcooled liquids have been investigated (both experimentally and numerically) by researchers [2][3][4][5][6][7][8][9].Numerical investigations using volume of fluid (VOF) method have been recently carried out by Ganguli [2].Further, Jeon et al. [3] developed a CFD model for condensation of the bubble and sensitivity tests for different bulk temperatures and liquid velocities have been performed.Various parameters of the bubble, like bubble velocity, interfacial area, bubble diameter, and bubble volume obtained from CFD simulations were compared.
The unavailability of data on the bubble dynamics in pool boiling systems is the main motivation of the current work.Further, the hydrodynamics of bubbles in subcooled International Journal of Chemical Engineering boiling systems have been studied by bubbles being formed on a horizontal plate.The dynamics of bubbles on vertical tube are; however, different from the one on a horizontal plate.The flow generated due to heating of a vertical tube is different than the one on horizontal plate.Another important aspect is the Rayleigh numbers considered in this work.It may be noted that this range of Rayleigh numbers has not been investigated for subcooled pool boiling.The current work is an effort to develop a CFD model and validate the results with experiments.

Experimental Setup and Measurements
2.1.Experimental Set-Up.The experimental setup is similar to the one used by Ganguli et al. [10].For producing a single bubble care has been taken for the bubble formation.The entire pipe has been buffed and an interstitial site has been made at the center of the pipe (Area = 3 mm 2 ) for bubble formation as shown in Figure 1.The bubble size was dependent on the heat flux and for understanding the bubble sizes trial runs were performed and reproducibility of bubble sizes were noted.

Photographic Measurement.
A digital video camera (Photron Fastcam 10KC, Resolution 512 × 480 pixels) together with halogen lamp (located at the opposite side of the camera) was used for taking pictures of the condensate films.All signals except video images were processed using the data acquisition system, consisting of A/D converter and a PC.The video images were analyzed by means of MATLAB software.All the instruments were calibrated before testing.

Computational Model
In this study, the VOF method [11] was used to simulate the bubble rise behavior in liquids having nonuniform velocity gradients.As the local pressure fluctuations were negligible compared to atmospheric pressure, the gas phase compressibility was neglected.Considering extremely low values of compressibility of water it was also considered as an incompressible fluid.

Assumptions.
The following assumptions were made for the study.
(2) No internal mass source are assumed.Hence the source terms for liquid and vapor have been considered to be equal.
(3) Bubble generation and behavior is a stochastic phenomenon.Yet the bubble generation at the same interstitial site has been tracked for the single bubble analysis, and the effect of other bubbles on the single bubble has been neglected.
(4) The bubble on the tube surface slides on the surface vertically upward for some distance and lifts off.
(5) The shear due to the thermal boundary layer as soon as it is developed is assumed to be constant for the lifecycle of the bubble since very little change in velocity gradients occurs in 500 ms which is the maximum time period for any bubble to travel the vertical distance under consideration.

Volume of Fluid (VOF) Method.
A two-phase model has been formulated in the Euler-Euler framework using VOF approach.The VOF model is a surface-tracking technique applied to a fixed Eulerian mesh.It is designed for two or more immiscible fluids where the position of the interface between the fluids is of interest.In the VOF model, the fluids share a single set of momentum equations and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain.The field for all variables and properties is shared by the phases and represent volumeaveraged values, as long as the volume fraction of each of the phases is known at each location.In each control volume the volume fractions of all the phases sum up to unity.Thus the variables and properties in any given cell are either purely representative of one of the phase, or representative of a mixture of the phases depending upon the volume fraction values.In other words, if the qth fluid volume fraction in the cell is denoted as ∈ q , then the following three conditions are possible: ∈ q = 0 the cell is empty (of the qth fluid), ∈ q = 1 the cell is full (of the qth fluid), 0 < ∈ q < 1 the cell contains the interface between the qth fluid and one or more other fluids.
Based on the local value of ∈ q the appropriate properties and variables are assigned to each control volume within the domain.Further, interphase mass transfer has been tracked by using the phase change model incorporated through a User Defined Function (UDF).For turbulence closure, realizable K-ε turbulence model has been incorporated.The governing equations for the entire region of the liquid and vapor are given later ( 5), ( 8), ( 9), ( 21) and (23).The explanation is given as follows.
The term on the right-hand side of (8) describes the mass transfer source term.This term can be calculated by using the mass transfer model (see the details in Section 3.3 below) and decides the extent of phase change (condensation).In the present system, we consider condensed water as the secondary phase and saturated steam as the primary phase.In the VOF modeling approach, the volume fraction equation is solved for the secondary phase.The primary phase volume fraction is computed based on the following constraint: ( The terms on the right-hand side of ( 5) describe the forces acting on the fluid element in the control volume.These are the overall pressure gradient, the stresses, the gravitational force, and the body force due to surface tension effect.The momentum equation ( 5) is dependent on the volume fraction of all the phases through the properties density ρ and viscosity μ: The effective thermal conductivity k e f f in the energy equation ( 9) and effective viscosity μ e f f in the momentum equation ( 8) is shared by both phases.S E is the source term in the energy equation that includes all the volumetric heat sources [W/m 3 ].In the present application, this term represents a sink that accounts for the latent heat of condensation during the phase change.Energy E is mass averaged as shown below: Further, the properties appearing in the transport equations are determined by the presence of the component phase in each control volume.In the present vapor-liquid system the volume fraction of the liquid is being tracked, the density and viscosity in each cell are given by ( 2) and (3).
A single set of Navier-Strokes equation for an incompressible Newtonian flow was solved.Momentum equation is given below: ( The CSF model of Brackbill et al. [12] was used to compute the surface tension force as a source term in the momentum equation for the cells containing the interface.The expression for F q is given by: where σ is the coefficient of surface tension, n is the surface normal which is proportional to the gradient of volume fraction, κ is the local surface curvature calculated as follows [12]: International Journal of Chemical Engineering The geometric reconstruction scheme [5] based on piecewise linear approach [4] was used for the reconstruction of the interface.In the VOF approach, the reconstructed interface is advected by solving an advection equation for the volume fraction of the dispersed phase.The equation for volume fraction is given as: Energy equation is given by: The liquid phase volume fraction is computed from the following equation:

Computation of Interface Mass Transfer and Turbulence
Modeling.A model for phase change in the present study has been used to simulate the condensation phenomena, which describes the interphase heat transfer of intrinsic flows and considers the heat transfer processes on each side of the phase interface.Therefore, the mass exchange source terms in mass equations and latent heat transfer term in energy equation have been incorporated in the modeled equation through UDF.At the liquid-vapor interface, the temperature should be equal to the saturation temperature.Heat content on vapor side across the interface is given by: where q v vapor heat flux across the interface, ∈ v is volume fraction of vapor, k v is thermal conductivity of vapor, ∇T is the temperature gradient at the liquid-vapor interface, and where "j" represents the volume fraction of vapor at each cell of the vapor-liquid interface.
Similarly, the Heat Flux on the Liquid Side across the Interface is Given by: where q l liquid heat flux across the interface, ∈ l is volume fraction of liquid, k l is thermal conductivity of liquid.The heat liberated during condensation (q cond ) in terms of latent heat is given by: where ṁ represents the total mass transfer from vapor to liquid.Equation ( 13) can also be written as The total mass transfer rate is the sum of each cell's mass transfer rate in the interface region.Let ṁ j represent each cells mass transfer rate proportional to vapor volume fraction in the jth cell in the interface region and V j is the volume of each cell.consisting of the vapor-liquid interface will have a spatial variation of the volume fraction ∈ at the region consisting of the vapor liquid interface, the mass transfer rate for two phase mixture is represented as follows [3]: where ∈ q, j is the volume fraction of the jth cell.The heat liberated due to condensation can also be expressed as: From ( 11)-( 16) we have Thus, the mass source in the vapor phase is represented by: Since there is no internal mass source, mass source for the liquid phase becomes ( Latent heat source for the energy equation becomes Turbulence equations are given by the SST k-ω model as follows.
Turbulent kinetic energy is Energy dissipation rate equation is given as follows Generation of turbulence due to buoyancy is given as: Dissipation of K and ω is given below: where, Production of ω is given by: Equations, (1) to (27) were solved using the commercial flow solver Fluent 6.2.

Solution Domain, Boundary Conditions, and Numerical
Schemes.The spatial derivatives in ( 5), ( 8), ( 9), (21), and (23) were discretized using the QUICK scheme [13] while a first-order implicit method was used for the discretization of the temporal derivatives.The pressure implicit with splitting of operator (PISO) algorithm [14] was used for the pressure-velocity coupling in the momentum equation.
A 2D rectangular solution domain (H = 150 mm, R = 150 mm) as shown in Figure 2 with nonuniform structured grid was used.A single phase simulation has been carried out for the first 20 seconds when the first bubble departs as observed in the experiments of Ra = 1.88 × 10 13 .A bubble of suitable diameter as calculated from the measured bubble size by high-speed camera has been patched as an initial condition and two-phase simulation using VOF method has been carried out.A no slip boundary condition was used at the stationary wall.A constant heat flux was given as a boundary condition, so that the spatially averaged heat flux can be monitored for longer times and the behavior of heat transfer coefficient can be monitored in time.The constant heat flux q cond has been calculated by using (13).
The effect of time step on bubble size was investigated by performing the simulations with four different time steps (1 × 10 −4 s, 1 × 10 −5 s, 1 × 10 −6 s, and 1 × 10 −7 s) for Ra = 1.8 × 10 13 .The difference in predicted bubble sizes for 1 × 10 −6 s and 1 × 10 −7 s was about 2% and therefore all further simulations were carried out using a time step of 1 × 10 −6 s.

Results and Discussion
In the present section we present the results of bubble movement and the lift force exerted on the bubble when a specific amount of heat input is provided.The present simulations have been provided after 60 seconds of run when the first bubble appears as per experimental observations.Rayleigh number (Ra) has been varied in the range 7.9 × 10 12 < Ra < 1.88 × 10 13 .For formation of single bubbles only small part of the whole pipe (Area = 3 mm 2 ) has been considered with interstices and the rest of the circumferential area has been buffed so that no bubbles can form at the surface.Simulations have been performed as per experimental conditions and the vertical bubble motion has been compared.The effect of heat flux which influences the bubble diameter and distance from wall has been studied.4.1.Grid Independence.The grid used for the simulations is shown in Figure 2.For single phase run the grid independence has been checked.It has been found that there is very little variation in the axial velocities at different axial locations as shown in Figure 3. Further, the grid independency has been tested for various grid size to bubble diameter ratio ((i) 1/10 and (ii) 1/20) for resolving the surface deformations.It has been observed that the deformation of the bubble surface is not much influenced by increasing the grid size.However, decreasing the number of cells caused divergence in the 2D simulations.

Variation of Axial
Velocity with Radial Distance.Figure 4 shows the variation of axial velocity with dimensionless radial distance.As the heat flux increases the axial velocity increases for a distance of 1 mm, decreases upto a distance of 4 mm, and then remains constant.A reasonably good agreement with experimental data has been observed.The experimental data was obtained with a hot film anemometer probe at various vertical positions near the wall.The constant temperature anemometer mode has been used which gives velocity series with time.

Effect of Heat Flux on Bubble Diameter and Flow Patterns.
In the case of subcooled boiling the effect of heat flux on the bubble diameter for a single bubble has been carried out in the present work.Figure 5 shows the variation in bubble diameter with increase in heat flux.It can be observed that, as the heat flux increases, the bubble diameter and the distance of the bubble from the wall after departure also increase.The increase in the bubble diameter is steep during the initial heat fluxes while the slope decreases as the flux increases.Figure 6 shows the velocity vectors due to the movement of the bubble at a time of 56 ms.As can be clearly seen significant amount of increase in flow is observed due to the bubble rise.The velocity vectors also match the photographs as in Figure 7(a).The rushing of the bulk fluid from the bulk can also be observed.

Bubble Tracking.
Movement of the single bubble has been tracked both with the help of high-speed camera (125 frames/s) and 2D CFD simulations.For different heat inputs (7.9 × 10 12 < Ra < 1.88 × 10 13 ) the bubble motion has been tracked.The bubble diameter was measured and the similar conditions (bubble diameter, heat input) were provided for the numerical runs.A single phase run was performed before the first bubble appeared and then a bubble of diameter as observed in the experimental run has been provided.Figure 7(a) shows the sequence of images of bubble movement (d B = 1.6 mm, Ra = 7.9 × 10 12 ) by both experimental and computational techniques.An area of liquid (30 mm × 40 mm) of the total area has been chosen International Journal of Chemical Engineering for investigation.Initially bubble adheres to the tube wall after which it moves in the upward direction.The bubble departs from the surface and migrates away from the wall after which it moves steadily without any movement upto a certain distance and then moves away from the wall.The shear for this case was measured as 1.3 s −1 .It can be observed that in corresponding VOF simulations the bubble gets deformed into ellipsoidal-shaped bubbles after departing from the surface.This denotes the need for a surface-tracking mechanism better than VOF.However, in experimental photographs (Figure 7(a)) the bubble shape is spherical till 72 ms while at t = 112 ms and t = 128 ms the experimental photographs also show ellipsoidal shape.Similarly, Figure 7(b) shows the sequence of images for higher bubble diameters (d B = 2.2 mm).It can be seen that the movement away from the wall after bubble departure is much more pronounced for 2.2 mm bubble.The shear rate considered here is 2.5 s −1 .However, several other reasons like the use of better turbulence models (Large Eddy simulations (LES)) and three dimensional (3D) simulations may also be needed.The movement of the bubbles is; however, well matched with the experimental observations.The bubble trajectory as shown in Figure 8 depicts that a smaller diameter bubble for a low-heat input moves steadily along the wall while a higher diameter bubble moves steadily along the wall and departs away from the wall after a certain vertical distance has been traversed.Tomiyama et al. [15] found that, in bubble columns, larger bubbles depart away from the wall while smaller bubbles moved towards the wall.The shear in the case of Tomiyama et al. [15] was 3.8 s −1 while here the shear and bubble diameter both are function of heat flux.Experimental visualization showed that, for the same shear rate, when runs were carried out at prolonged time periods (t = 120 s) the bubble tended to move towards the wall adhere to the surface and again depart.This process continued till it reached the top free surface of the pool.The analysis presented in this subsection suggests that a more detailed study of the dynamics occurring in the thermal boundary layer and the turbulent structures near the wall has a major role to play in the bubble movement.The study is; however, out of scope of the present work.

Conclusions
Bubble dynamics in subcooled pool boiling systems have been evaluated for (7.9 × 10 12 < Ra < 1.88 × 10 13 ).The bubble movement has been found to be lateral for a certain vertical distance and bubbles departed when they approached the top surface.Bubble trajectories were plotted for the vertical distance in which the bubbles were found to be in the wall regime with a shift to the neutral regime as defined by Tomiyama et al. [15].Bubble diameter increases steeply with increase in Ra number (Ra < 1.5 × 10 13 ) with no significant increase with an increase in Ra number.The single bubble considered moves within the boundary layer until it is near the top surface where it departs.The developed CFD model can be extended to study various phenomena like (a) study the effect of multiple bubbles on hydrodynamics and heat transfer in subcooled boiling systems, (b) calculation of non-drag forces like lift and virtual mass forces in such systems, and (c) studying coalescence and breakup phenomena in such systems.

Notation d B : Equivalent bubble diameter, mm E:
Total energy supplied to the system, kJ E q : Energy of fluid, kJ F q : Force because of the pressure jump at the interface, N g: Gravitational constant, m s −2 G b : Generation of turbulence due to buoyancy, m 4 s −1 G K : Generation of turbulent kinetic energy due to mean velocity gradients (kg•m −1 s 3 ) H: Column height, mm H: Specific enthalpy as in energy equation in (9), kJ/kg I: Identity matrix, (-) K: Turbulent kinetic energy as in (21), m 2 /s 2 k: Thermal conductivity, W m Thermal conductivity of liquid, W/mK k v : Thermal conductivity of vapor, W/m K ṁ: Total mass flow rate of steam, kg/s ṁ j : Mass flow rate at the jth cell, kg/s n: Surfacenormalvector,(-) n: Unit normal vector, (-) p: Totalpressure,Pa p : Time averaged of total pressure, Pa q cond : Condensation heat flux in (18), W/m 2 q l : Liquid heat flux across the interface, W/m 2 q v : Vapor heat flux across the interface, W/m 2 R: Radius of the tank, m r: Spatial co-ordinate in radial direction S ∈l : Mass source in liquid phase, kW/m 3 S ∈q : Mass source in qth phase, kW/m 3 Volumefraction ∈ q : Volume fraction of fluid, (-) ∈ q,J : Volume fraction of fluid at jth cell, (-) ∈ l : Volume fraction of the liquid phase, (-) κ: Surface curvature λ cond : Latent heat of condensation kJ/kg

Subscripts and Superscripts
eff : Effective i, j, k: Axis indexes of spatial co-ordinates l: Liquid q: Phase v: Vapor r: Cylindrical space coordinate for radial direction t: T urbulent z: Cylindrical space coordinate for axial direction.

Figure 1 :Figure 2 :
Figure 1: (a) Schematic of a dimensionless rectangular enclosure with leftside wall heated.(b) Schematic of a vertical rectangular enclosure with leftside wall heated and a bubble patched at the interstice.

Figure 5 :
Figure 5: Variation of bubble diameter and maximum radial distance of bubble from the wall with Rayleigh number.radial distance from the wall bubble diameter.

Figure 6 :
Figure 6: Velocity vectors showing the movement of the bulk fluid due to bubble movement.