Lift force on single bubbles in subcooled pool boiling systems: Experimental measurements and CFD simulations

Lift coefficient has been evaluated for a subcooled pool boiling system with water as continuous phase and vapor (steam) bubble as dispersed phase generated on vertical walls with constant heat flux for a specific Rayleigh number (Ra) range (8x1012 < Ra < 2x1013) for the first time. Experimental investigation has been done with High Speed Camera while two phase CFD simulations have been performed using Volume of Fluid (VOF) method. The Morton number (M) for this system is log10M = -10.6 and Eötvös number (Eo) range is 0.2 < Eo < 1.3 while the shear rates have been maintained in the range 2 ≤ Sr ≤ 2.8. Bubble size of the bubbles arising from the sites on vertical surface have been determined by both experiments and simulations. The CFD predictions for bubble size and lift coefficients agree well with experimental measurements with a deviation of less than 10%. Further, empirical expressions from experimental data for the estimation of bubble size and lift coefficients have also been presented. The data for lift coefficients obtained is different from the ones of bubbles generated by single sparger in air-water systems.


Introduction
Boiling phenomena has been studied predominantly on horizontal plates or tubes in saturated stream both using experimental and numerical techniques.Studies of boiling on vertical surfaces have been limited in the reported literature.Modeling boiling using CFD is carried out using two approaches nmely Euler-Euler and Euler-Lagrange approaches.In the Euler-Euler approache the interactions between vapor and liquid phase is accounted by a source term representing an overall volumetric force FGL in the momentum equation.This force is the summation of all interfacial forces namely drag, lift, turbulent dispersion force and wall lubrication force.The vapor bubble formation has the following steps: (i) it forms a nucleation site on the vertical surface (ii) grows to a certain size (iii) departs from the wall after attaining a certain size.The overall volumetric force is given by the equation below,: .
where drag ( ), lift ( ), turbulent dispersion force and wall lubrication force ( ) Researchers have observed that as the bubble departs from the nucleation site (in case of orifices placed at the center at bottom of a liquid column), the bubble may take a central path parallel to the wall at the center or deviating away or nearer to the wall depending on the shear induced on the bulk liquid.The side way movement during its rise is due to the fact that during its motion the bubble experiences a lift force.The lift force investigations in the presence of thermal gradients however differ from the shear induced flows.Both systems however have velocity gradients in the continuous phaseand lift force plays a crucial role in determining the lateral migration of bubbles.

Previous work
Estimation of lift force in terms of experimental measurements have been done in the gas-liquid flows while works with vapor-liquid flows have been limited.The lift force exerted on a single bubble rising in a liquid is generated by three different mechanisms: magnus lift due to bubble rotation [1][2][3], Saffman lift force on non-rotating bodies due to shear in the liquid [4][5][6] and wake induced lift [7] caused by bubble wakes [8].The most important parameters in the lift force evaluation are lift coefficient, bubble diameter and the velocities of both phases.The details of gas/vapor liquid (GL) systems on which experimental measurements and analytical modelling have been carried out, the experimental operating conditions and geometric details, typical range of bubble sizes in the GL systems and the range of lift coefficients estimated has been tabulated in Table 1 while the prominent correlations for lift coefficients have been listed in Table 2.The shear induced lift force is defined as [9][10][11]: ( where is the lift force, CL is the lift coefficient, ρL is the liquid density, is the velocity and dB is the bubble diameter.Auton [10] carried out theoretical investigations to arrive at an expression for CL.He assumed a clean spherical bubble in an inviscid rotational system.However, predicted value of lift coefficient (CL=0.5)by Auton [10], did not match with the experimental observations of Kariasaki [7].The variation in the value of CL indicated that the lift coefficient was strongly dependent on the flow condition.Later, Tomiyama et al. [6] performed experiments to estimate the lift force acting on a single spherical/ellipsoidal air bubble of different sizes (2.8 ≤ dB ≤ 5.68 mm) rising in a simple shear flow in a viscous liquid solution ( = 0.018 -0.089 kg/m.s).The authors proposed different correlations to calculate CL, which were valid only for viscous systems.The applicability for less viscous liquids like water was however not verified.Bothe et al. [12] performed VOF simulations of air bubble rising in viscous liquid imposed with linear shear corresponding to the experiments reported by Tomiyama et al. [6].The predicted CL for viscous system was in good agreement with the measured CL of Tomiyama et al. [6].The authors also investigated the effect of bubble diameters (range of 2.8 -10 mm), liquid properties ( = 0.018 -0.1 kg/m.s;= 0.061 -0.8 N/m) on CL and proposed a correlation for CL as a linear function of EoH based on the horizontal bubble diameter.Unfortunately, this correlation was found to be valid for only higher liquid viscosities (4 ≤ EoH ≤ 8, -3.6 ≤ log Mo ≤ -5.7), and could not be applied to less viscous fluids such as water.Rabha and Buwa [13] investigated the effect of bubble size/shape and the effect of neighbouring bubbles on the magnitude and direction of lift force.Further, they investigated the lift force acting on single bubbles rising in low viscosity fluids (corresponding to the wobbling regime, 1.1 ≤ Eo ≤ 8.72, -10.6 ≤ log Mo ≤ -14.6).They observed that unlike the steady lateral migration of single bubbles with single characteristics value of CL (+ve or -ve depending on dB) for viscous systems, the bubbles were found to oscillate around the center line and the instantaneous CL was found to fluctuate in both +ve and -ve directions.Van Helden et.al. [14] carried out experimental measurementson the bubble detachment near a plane in a vertical rectangular channel.The authors used two different fluids for bubble generation namely (1) vapour bubble (2) nitrogen bubble.The authors observed that the steam bubbles traveled towards the bulk liquid (away from the wall) while the nitrogen bubbles travel parallel to the wall in the bulk liquid.The authors deduced two notable conclusions from their study namely (i) the bubble detachment depended significantly on the bulk density of continuous phase ith radius of bubble detachment decreasing with increase in bulk density (ii) theoretically derived expressions for lift coefficients (Auton [10]) were unable to predict lift coefficients for bubbles moving near the wall Felton and Loth [15] carried out experimental investigations on dynamics of spherical air bubbles in water on the basis of bubble velocities in terms of Stokes number (ratio of characteristic time scale of bubble to that of flow which depends on fluid velocity and length scale (example boundary layer thickness)).The authors concluded that for high Stokes number, the bubbles got accumulated along the wall and a Gaussian distribution of bubbles across horizontal plane for low Stokes number.Further, three different types of bubble trajectories namely sliding bubbles, bouncing bubbles and freedispersion bubbles were also observed.Tomiyama et.al. [6] experimentally investigated the estimation of transverse lift force acting on single bubbles in simple shear flows on glycerol-water system.The authors observed that the net transverse lift coefficient for small bubbles is a function of the bubble Reynolds number.However, for larger bubbles the lift coefficient was found to be a function of modified Eötvös number.The authors estimated expressions for critical bubble diameter as a function of lift coefficient.Rabha and Buwa [13] carried out numerical investigations on the effects of the size of bubbles, the number of neighbouring bubbles and the shape of bubbles on the magnitude and direction of the lift force.The authors found that the direction of lift coefficient depended on the viscosity of the system with oscillatory behaviour in low viscosity systems and steady behaviour in high viscosity systems.Further, as per the authors, the magnitudes of lift coefficient depended on bubble-bubble interactions.An increase in the magnitude was observed when bubble-bubble interactions were higher irrespective of the properties of the continuous/stagnant liquid phase.Aoyama et.al. [16] experimentally investigated the effects of the bubble shape, the liquid velocity gradient and the fluid properties on the magnitude of lift coefficientfor two different bubble shapes namely spherical and ellipsoidal.The authors observed that the lift coefficient of spherical bubbles at low Reynolds number increases with an increase in Shear Rate while it decreases with an increase in Reynolds number.Hessenkemper et.al. [17] experimentally investigated the effects of inorganic surfactant (NaCl), on the lift force for single bubbles.The authors concluded that the lift coefficients for small bubbles got affected by the addition of inorganic surfactantwhile those for larger bubbles got affected only upto a certain critical concentration (0.1 mg•l-1), after which it weakly affected the CL values.Hidman et al. [18] numerically investigated two different parameters: (1) the lift force acting on a freely deformable bubble rising in a linear shear flow (2) the lift force scaling with the undisturbed shear rate for different lift force mechanisms.The authors investigated the lift coefficient as a function of various non-dimensionelss numbers Eo, Ga, Sr and Re.The authors observed that for highly viscous flows and for significant bubble deformations, the Shear Rate becomes an important study parameter for modeling of the lift force coefficient.Subburaj et al. [21] carried out numerical studies for lift coefficient evaluation considering for the variation in the bubble spreading mechanism for different aspect ratios of the column with homogenous and heterogenous bubbly flows.The authors observed that the local gas fraction and the Shear Rate were simultaneously correlated with the lift coefficient.The authors also observed that decreasing decrease in the gas fraction reduced the magnitude of the lift force, while a bubble in the proximity of another bubble experienced higher lift.This is a contradictory observation as compared to [22].Table 2 Correlations for Lift coefficients

Objective of present work
The review of published literature shows that most of the work carried for evaluation of lift coefficients is for air-water systems and those for vapor-liquid systems are not available.Further, the lift forces are evaluated for higher viscous systems for different ranges of Eötvös and Morton numbers.The generation of bubbles have been mainly though single spargers and the shear has been induced using mechanical systems like moving belt etc.In the present work, a system different than the conventional GL system has been chosen.Here, the shear induced was due to the thermal boundary layer on a vertical wall maintained at constant heat flux.Further, the bubbles are generated from the nucleation sites created on the wall.The three main objectives of the evaluation of lift coefficient are 1.Determination of bubble size as a function of Rayleigh number 2. Determination of empirical correlation for determination of lift coefficient as a function of dimensionless numbers 3. Building a robust CFD model capable of predicting lift coefficients for vapor-liquid systems.

Experimental and CFD details
For experimental measurements of the bubble size, shape and trajectory a High speed Camera (Photron Fastcam 10KC, Resolution 512×480 pixels) was used.The images were analyzed by means of MATLAB software.All the instruments were calibrated before testing.For CFD simulations, a bubble on the vertical wall (represented as line in a 2D axisymmetric geometry) has been patched and the VOF multiphase model has been implemented with turbulence given by SST k-w model.The governing equations of momentum, energy, turbulence and VOF are same as used by [25] in Ansys Fluent 14.The spatial derivatives for the governing equations of momentum were discretized using second order upwind scheme, energy was discretized using first order upwind scheme while turbulence were discretized using the QUICK scheme [26].A first-order implicit method was used for the discretization of the temporal derivatives.The Pressure Implicit with Splitting of Operator (PISO) algorithm [27] was used for the pressure-velocity coupling in the momentum equation.The procedure for transient simulation is as follows: 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 = 3.76E+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 used.A no slip boundary condition was used at the stationary wall.The effect of column depth (D) on the bubble rise velocity (VR) and lift coefficient (CL) has been studied.The time step was kept as 1x10 -6 s for all simulations.

Effect of Ra number on bubble size
In the present investigation of vapor-liquid flows it has been observed that bubble size during inception and movement is a function of Rayleigh number.It was also observed that during the 55% of the total time the vapor bubble nearly retains its size with negligible condensation since it remains in the region of thermal boundary layer for all Ra numbers except the highest Ra number considered.For the highest Ra number, the time duration is 38% of the overall time after it starts condensing.Hence, the lift force calculation has been done considering the above aspects during bubble dynamics and the time of condensation and change in bubble size as a function of both spatial and temporal variation has not been considered.Figure 1 shows the bubble size as a function of Ra number for the Ra number range 8 x 10 12 < Ra < 3.76 x 10 13 .As the Ra increases, the bubble size linearly increases till Ra = 1.8E+13.The bubble size variation tapers down with further increase in Ra number.The experimental and simulation predictions show a good agreement with a deviation of 4-7% with CFD simulations underpredicting the experimental data.This is attributed to the two dimensional nature of simulations.An empirical correlation for the bubble size as a function of Ra is given below:  !/ !,#$% = 0.203 &.(() 8 x 10 12 < Ra < 3.76 x 10 13 ; dB, max = 3 mm (3)

Evaluation of lift coefficient
In the present section we present the results of the lift force exerted on the bubble when a specific amount of heat input is provided are reported.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.9e+12 < Ra < 3.76e+13.The lift force has been calculated using the expression below The lift coefficient calculations (both experimental and numerical cases) have been carried out according to Eq. 4. The major difference between air bubbles in water in wall regime is the existence of thermal boundary layer and the intrusion layer that forms at the top surface.Further, since the bubble motion starts from the wall and the inception of the bubble, its formation are dependent on nucleation site, it is different than single bubble motion from the center.For the present case the bubble motion has been explained in our previous work (Ganguli et al. [25]).The lateral migration, length or distance from the wall depended on the relative velocity gradient.A positive lift coefficient has been observed for the time period considered and for all the heat fluxes under consideration.In the present case of subcooled pool boiling the bubbles should be having positive values in the wall regime.Although positive values are observed both in simulation and experiments and the bubbles move in rectilinear motion at low heat fluxes, the important observation is that after a certain vertical distance the bubbles depart away from the wall.This suggests that there is a good balance between the buoyancy, surface tension and viscous force upto a certain vertical distance after which the buoyancy w forces overcome the viscous resistance and move away from the wall.Thus, similar to Tomiyama et al. [6] we can predict that the bubble is in the wall regime depending on the shear due to the hydrodynamic and thermal boundary layers.Figure 1A shows the lift coefficient variation with Reynolds number for Ra number range 8 x 10 12 < Ra < 3.76 x 10 13 .The simulation results agree reasonably well with the experimental measurements with 6% under-predictions which is attributed to the 2D simulations carried out.An empirical correlation for lift coefficient as a function of Re has been given below:  * = −3 − 05 ) + 1 − 03 + − 0.00026 − 0.0034 log10 M = -10.6;0.2 < Eo < 1.3; 8 x 10 12 < Ra < 3.76 x 10 13  (5) Interestingly, the highest values lift coefficient (CL) are around 0.1 for a bubble size of 3 mm.Further, the Reynolds numbers (based on the resultant velocity and bubble diameters obtained) are in the range of 5 < Re < 30 which is similar to the Reynolds numbers considered by Tomiyama et al. [6] but due to the shear mechanism and different Morton number the lift coefficients magnitudes are much lower than the ones observed by Tomiyama et al. [6].By the lower values of lift force it may be attributed to several factors like the mass transfer taking place between the liquid and vapor bubble, the vapor shear rate and the unsteady movement of the bubble (since the lift coefficient is averaged over time and the vapor bubble after detaching has lower lift coefficients than when it has traveled a certain distance) hence giving us values which are one third of the values if CL in published literature.Figure 2B shows the variation of CL with respect of Eö.The tapering of CL values after Eö value of 0.6 can be observed.This is a significant result considering the fact that the CL values for this range have not been reported for such lower range of Eö numbers at the specified Morton number and the range of Ra numbers considered.The lower magnitude of lift force and positive sign is attributed to the explanation provided previously in Figure 2.However, this aspect needs to be looked in more detail and has not been included presently.

Figure 1
Figure 1 Variation of bubble diameter ratio with Rayleigh number (Ra).Symbols denote experimental data while line denotes CFD predicted data

6 )Figure 2
Figure 2 Variation of lift coefficient with (A) Reynolds number.(B) Eötvös number (Eo).Symbols denote experimental data while line denotes CFD predicted data

Table 1
Details of various geometric and operating parameters for lift coefficient measurement