Numerical Investigation of Nucleate Boiling Flow in Water Based Bubble Bumps
1Unit of Thermic and Thermodynamics of the Industrial Processes, National Engineering School of Monastir, University of Monastir, Monastir, Tunisia
2Department of Mechanical Engineering, Universitat Rovira i Virgili, Tarragona, Spain
To cite this article:
R. Garma, M. Bourouis, A. Bellagi. Numerical Investigation of Nucleate Boiling Flow in Water Based Bubble Bumps. International Journal of Fluid Mechanics & Thermal Sciences. Vol. 1, No. 2, 2015, pp. 36-41. doi: 10.11648/j.ijfmts.20150102.14
Abstract: In this paper, numerical simulations of nucleate boiling flow bubble pumps are conducted with the commercial CFD (Computational Fluid Dynamics) package Ansys-Fluent. The Eulerian multiphase flow framework model was used to model the phase’s interaction. User-Defined Functions (UDFs) are provided to compute the wall heat transfer and to calculate inter-phase heat and mass transfer. The heat flux from the wall is divided into three parts according to a wall heat partitioning model based on three mechanisms including convective heat for heating the bulk liquid, evaporative heat for generating vapor and quench heat for heating the liquid in the nucleation sites. The rate of vapor formation is obtained by adding the mass exchange at the bubble surface and the bubble formation due to heat flux at the wall. Constant heat fluxes are applied to the stainless-steel made tube wall. In the simulation results we discuss the radial temperature distribution and the radial and axial profiles of the vapor void fraction in the pipe to localize the onset of vapor generation in the pump tube.
Keywords: Boiling Flow, Bubble Pump, CFD
In the diffusion absorption refrigerators (DARs) invented by the Swedish engineers von Platen and Munters  in the 1920s, the bubble pump is considered as one of the most important components. Therefore, significant attention has been devoted to this device to improve its performance, which will contribute to the performance of the whole system. The thermally driven bubble pump is a simple vertical tube which can be powered by waste heat or solar thermal energy. When the liquid solution is heated up, vapor bubbles are generated causing pressure difference between the bottom and top of the tube. As result, natural circulation of the liquid solution is driven up through the tube and so boiling flow takes place inside the bubble pump.
Several researches have foxed on numerical investigation of nucleate boiling flow for various applications (thermosyphon, nuclear engineering, PWRs, BWRs, electronics cooling…). CFD modeling of flow and heat transfer in a thermosyphons was carried out by Alizadehdakhel et al. . In this research, the effects of the fill ratio (amount of working fluid) and the heat pipe input were considered. This study shows that the CFD code is a useful tool to model and explain the flow coupled to heat transfer in a thermosyphon. Using the same software CFD code Fluent, numerical simulation of nucleate boiling for power electronics cooling applications was conducted by Narumanchi et al.  to enhance heat removal from electronics packages. Li et al.  modified the basically by incorporating new closure correlations in order to investigate numerically the boiling flow of nitrogen in a vertical tube. In the same way, multidimensional modeling of vertical upward sub-cooled boiling flow was undertaken by Koncar et al. .
In this work a general-purpose CFD code Ansys-Fluent was used for simulations of upward nucleate boiling flow. For modeling the boiling phenomena user defined subroutines (UDFs) are employed. Radial temperature profiles and void fraction distributions for different heat fluxes applied to the stainless-steel made tube wall are discussed.
2. Model Description
Mathematical model utilized in this paper was developed and then applied in CFD codes to be finally implemented in Ansys-Fluent via user-defined functions (UDFs) in conjunction with the Eulerian multiphase model in which the conservation equations are written for each phase, liquid and vapor. The following is a summary of the main model equations for a certain phase q.
2.1. Conservation Equations
The mass conservation equation for q phase is given by:
where α is the phase volume fraction of phase q, ρ, the density, v the velocity vector, mpq characterizes the mass transfer rate from phase q to phase p and n, the number of phases.
The momentum conservation equation for q phase is given by:
where τ s is the stress-strain tensor, Fpq an interaction force between phases, Fq is an external body force, Flift,q is a lift-force, Fvm,q is a virtual mass force, p is the pressure shared by all phases, and g is the gravitational acceleration vector.
The energy conservation equation for q phase is given
Where h is the specific enthalpy, is the heat flux vector, is the source term, Q is the intensity of heat exchange between the different phases, and hpq is the difference in the formation enthalpies of phases p and q.
The heat exchange between phases must comply with the local balance conditions, as well as the interfacial mass, momentum interfacial exchange:
2.2. Mass Equation
The rate of vapor formation per unit volume in Eq. 1 can be written as:
The first term is the mass exchange at the bubble surface and the second term represents the bubble formation due to heat flux at the wall. Here:
is the interfacial heat transfer coefficient calculated using the Ranz-Marshall correlation, where is the diameter of the secondary phase (vapor bubble), is the liquid heat conductivity. Re, Pr and Nu are Reynolds number, Prandtl number and Nusselt number respectively.
is the interfacial area density, with is the vapor volume fraction and ;
· = evaporating heat flux calculated from the RPI model [6, 7];
· = latent heat per unit mass
· = interfacial area density of the wall surface.
Subscripts l, v and s mean liquid phase, vapor phase and saturation state, respectively.
2.3. Momentum Equation
The interfacial drag force between liquid and vapor phases per unit volume is calculated as:
where Cd is the drag coefficient determined by choosing the minimum of the viscous regime and the distorted regime:
The lift coefficient is calculated as (Moraga et al. ):
This lift coefficient combines the opposing action of two lift forces:
• Classical aerodynamic lift force resulting from the interaction between bubble and liquid shear,
• Lateral force resulting from the interaction between bubble and vortices shed by the bubble wake.
is the bubble Reynolds number and,
is the bubble shear Reynolds number.
2.4. Turbulence Model
The mixture turbulence model, default multiphase turbulence model, was used. It represents the first extension of the single-phase k-ε model. In the present case, using mixture properties and mixture velocities is sufficient to capture important features of the turbulent flow. The equations describing this model are respectively:
where and are the mixture density and velocity, is the turbulent viscosity, is the production rate of turbulence kinetic energy, is the turbulent kinetic energy, is the dissipation rate. and are constants.
This model contains two additional terms describing additional bubble stirring and dissipation. is the bubble-induced turbulence in the turbulent kinetic energy equation and the bubble-induced dissipation in the dissipation rate equation:
The turbulent diffusion force is calculated as :
with is the turbulent dispersion coefficient . This force simulates liquid turbulence induced diffusion of bubbles from the wall into the liquid bulk.
2.5. Wall Boiling Model
According to the basic Rensselaer Polytechnic Institute (RPI) model , the total heat flux from the wall to nucleate boiling consists of three different components, namely the convective heat flux, the quenching heat flux, and the evaporative heat flux:
The heated wall surface is subdivided into a portion Ω (0≤ Ω≤1) covered by nucleating bubbles and the rest of the surface area (1-Ω), is covered by liquid.
The convective heat flux is expressed as:
where is the single phase heat transfer coefficient, and and are the wall and liquid temperatures, respectively.
The quenching heat flux models the cyclic averaged transient energy transfer related to liquid filling the wall vicinity after bubble detachment, and is expressed as:
where is the bubble departure frequency, , the thermal conductivity, , the speciﬁc heat, and , the density.
The evaporative flux is given by:
where is the bubble departure diameter, , the vapor density, , the active nucleate site density.
These equations need closure for wall boiling parameters:
2.5.1. Bubble Departure Diameter
The default bubble departure diameter (mm) for the RPI model is based on empirical correlations and estimated as:
with is the sub-cooling temperature.
2.5.2. Nucleate Site Density
The nucleate site density is represented by a correlation  based on the wall superheat as follow:
2.5.3. Frequency of Bubble Departure
The bubble departure frequency is calculated as:
2.5.4. Area Density
The effective wall area occupied by boiling sites definition is based on the departure diameter and the nucleate site density:
where , and is the Jacob number .
2.5.5. Bubble Diameter
The bubble diameter in the free stream is correlated with the local sub-cooling temperature
3. CFD Modeling
The commercial CFD code Ansys-Fluent  was used as a frame to perform the simulation. The interfacial forces models and the wall boiling model described previously were implemented in the code through User-Defined Functions (UDFs). The flow domain is shown schematically in Fig. 1. The pipe is 1m in length and 10 mm in diameter. The sub-cooled water enters the system at the bottom, and then boils due to the constant heat flux supplied from the pipe walls.
3.1. Mesh Geometry
Only half of the pipe is considered due to the symmetry at the pipe axis. GAMBIT software was used to build and mesh a two-dimensional computational domain (2-D geometry) with 5x500 uniform rectangular cells.
3.2. Initial Boundary Conditions
In the present work, the saturation temperature is fixed to 425.15K. Fully-developed profile of velocity and sub-cooling temperature (Tsub = 5K) is applied at the inlet. The gravity acceleration is 9.8 m/s2. Wall thickness is fixed to 2mm.
3.3. Solution Techniques
Unsteady state calculations with a time step of 0.1 s were fixed for all cases. SIMPLE algorithm was carried out for the calculations of the pressure velocity-coupling and the first order upwind calculation scheme was performed for the discretization of momentum, energy and volume fraction equations.
4. Results and Discussion
The Prediction of axial and radial void fraction distributions and radial temperature profiles for boiling flow in a vertical tube is conducted using the previous described model. Specially seven different wall heat flux are considered from 25 to 150 kWm-2.
Fig. 2 illustrates the void fraction distribution along the tube for 90kW/m2. It is seen that no boiling occurs at the inlet zone which corresponds to a single phase heating zone. In fact, the sub-cooled liquid enters the tube, is heated at the wall but vapor is not yet formed at this point. Then, when the sub-cooled water reaches the saturation temperature boiling is initiated.
Fig. 3 depicts the effect of the heat inputs on the void fraction profiles along the pipe axis. One can see that the higher the heat flux is, the lower is the start-up of boiling and the higher is the void fraction. In fact, when the heat input is increased from 25 to 150kW/m2, the onset boiling length is reduced from 43 to 1.6 cm and the void fraction is increased from 0.01 to 0.7 at the tube outlet.
Fig. 4 illustrates radial void fractions at different axial locations for two different heat inputs. One can observe, as expected, that the void fractions are higher adjacent the wall and progressively fall down towards the center for the reason that bubbles form near to the heated wall and then move towards the middle part of the tube.
Radial temperature profiles at different tube locations for 75kW/m2 are depicted in Fig. 5. One can notice that the liquid temperature profiles show the same behavior for the various locations. Actuality, the liquid at the center of the tube is sub-cooled even at the outlet (z = 1 m) and it is superheated near to the wall. We can see also that the sub-cooling level is decreasing throughout the tube’s length.
Liquid and vapor velocities arrangement for 90kW/m2 are illustrated in Fig. 6. One can remark that the liquid velocity increases slightly, it goes from 1.2 to 1.3 m/s from the entrance to 0.4 m, then increases sharply to the rest of the tube it reached 2.9 m/s at the exit. The steam follows almost the same profile as the liquid .
Nucleate boiling flow of water in vertical tube was carried out in this paper using commercial CFD (Computational Fluid Dynamics) package Ansys-Fluent. User defined Functions (UDFs) are employed to model the boiling phenomena. Radial temperature profiles at different tube locations and void fraction distributions for different heat fluxes are calculated and discussed. It was found that the onset boiling point is reduced from 43 to 1.6 cm and the void fraction at tube’s outlet is increased from 0.01 to 0.7 when increasing the wall heat input from 25 to 150kW/m2. The liquid temperature in the middle part of the tube is still sub-cooled along the pipe and it is superheated near to the wall so bubbles migrate from the heated wall to the center of the tube.