Progress in Nuclear Energy 89 (2016) 120e131

Contents lists available at ScienceDirect

Progress in Nuclear Energy

journal homepage: www.elsevier.com/locate/pnucene

Numerical simulation of bubble condensation using CF-VOF

N. Samkhaniani, M.R. Ansari*
Faculty of Mechanical Engineering, Tarbiat Modares University, P.O. Box 14117 13116, Tehran, Islamic Republic of Iran

a r t i c l e i n f o a b s t r a c t

Article history: In present study, subcooled boiling is simulated using color function volume of fluid (CF-VOF) method.
Received 3 November 2015 For this purpose, energy equation and Tanasawa mass transfer model accompanied with some suitable
Received in revised form source terms are implemented in OpenFOAM solver (interFoam). The surface tension between vapor
27 January 2016
eliquid phases is considered using continuous surface force (CSF) method. In order to reduce spurious
Accepted 3 February 2016
current near interface, a smoothing filter is applied to improve curvature calculation. The variation of
Available online 23 February 2016
saturation temperature in vapor bubble with local pressure is considered with ClausiuseClapeyron
relation. The numerical model is validated with one-dimensional Stefan problem.
Keyword:
Condensation
The shape and life time history of single vapor bubble condensation are verified with existing
Subcooled boiling experimental data. Computational study shows bubble life time is nearly proportional with bubble size
VOF and it is prolonged at bubble swarm motion. The present study reveals some fundamental characteristics
OpenFOAM of single and multiple vapor bubble condensation and is expected to be instructive for further
Bubble applications.
© 2016 Elsevier Ltd. All rights reserved.

1. Introduction analyzed the bubble deformation and life time (Kamei and Hirata,
1990a,b). Chen and Mayinger (Chen and Mayinger, 1992) studied
The bubble condensation is one of the fundamental issues in the the heat transfer at the interface of condensing vapor bubble in a
subcooled flow boiling to describe the heat and mass transfer. It is subcooled liquid of the same substances with ethanol, propane,
encountered in many industrial applications such as nuclear re- R113 and water. Harada (Harada et al., 2010) carried out the visu-
actors. For nuclear reactors, the bubble dynamics can greatly in- alization experiments to investigate the dynamics of vapor bubbles
fluence the reactivity feedback characteristics of coolant which generated in water pool boiling. In the experimental studies, bubble
brings more challenges for reactor safety analysis. The bubble size, behavior with bubble size history, shape, velocity, collapse time and
bubble shape and void fraction change continuously in bubble interfacial heat transfer coefficient was investigated in quiescent
condensing process, and it affects flow structure around bubble. In (Brucker and Sparrow, 1977) and upward flow (Lucas and Prasser,
order to understand the subcooled flow boiling, it is a challenge to 2007). It was found that the collapsing points of time and height
obtain an extensive knowledge on the condensing bubbles are proportional with the variation of pressure and temperature
behavior. difference in quiescent liquid (Brucker and Sparrow, 1977). Bubble
There have been many experimental analyses on the condensing collapse during condensation in immiscible liquids can be
bubble behavior (Kamei and Hirata, 1990a,b; Sideman and Hirsch, controlled either by inertia or by heat transfer mechanisms. At a
1965; Sudhoff et al., 1982). In experimental investigations, visuali- high subcooled temperature, bubble rapidly collapses. This process
zation is a common method to analyze the bubble condensing controlled by inertial. On the other hand, if subcooled temperature
process. The surface area, volume and vapor content of a rising is relatively low, the bubble life time is longer and the process will
bubble are determined through visualization. Sideman and Hirasch be controlled by the heat transfer at the interface. The bubble
(Sideman and Hirsch, 1965) studied free rising of the isopentane collapse rate is generally assumed to be controlled by the internal,
condensing bubbles at subcooled water. Kamei and Hitara have external thermal resistances and the temperature driving force.
However, it can be affected by many parameters such as: working
fluid, miscibility, bubble shape, bubble size, fraction of non-
condensable gases, surface mobility and etc. For short duration of
* Corresponding author.
E-mail addresses: [email protected] (N. Samkhaniani), mra_1330@ experiment and complexity of phenomenon, it is impossible to
modares.ac.ir (M.R. Ansari). obtain detailed information about condensing process through

http://dx.doi.org/10.1016/j.pnucene.2016.02.004
0149-1970/© 2016 Elsevier Ltd. All rights reserved.
 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 121

Nomenclature !
Ub Bubble velocity, m/s
!
Uc Compressive velocity, m/s
Ca Compression factor, e !
Cp Specific heat, J/kgK Ur Relative velocity, m/s
dL Liquid thermal diffusivity: k/rCp, m2/s
D Equivalent diameter, m Greek symbols
D0 Bubble Initial diameter, m a Volume fraction factor
!
g gravity acceleration, m/s2 k Interface curvature, m1
HLG Latent heat, J/kg m Dynamic viscosity, Pa s
k Thermal conductivity, W/mK r Density, Kg/m3
M Molar mass, Kg/Kmol
000 Subscripts
m_ Condensate mass flow rate per unit volume, Kg/m3s
G Gas (vapor) phase
P Pressure, Pa
L Liquid phase
Prgh Dynamic Pressure, Pa
sat Saturation condition
R Universal gas constant, J/mol K
sub Subcooled condition
T Temperature, K
!
U velocity, m/s

experimental data because the shape and the area of bubble

8 !
interface are in exposure of rapid changes. Therefore, numerical
! VLiquid < 1 x 2Liquid
simulation of condensing bubble is vital as complement to exper- aL ð x ; tÞ ¼ ¼ 0 < aL < 1 !x 2interface (1)
iments. Tian et al. (Tian et al., 2010) simulated single steam bubble V : !
0 x 2gas
condensation behaviors in subcooled water using the Moving Par-
ticle Semi-implicit (MPS) method in axisymmetric domain. Chen The thermo physical properties of two immiscible fluids such as
et al. (2010); Chen et al. (2011) applied MPS method in 2- viscosity (m), density (r) and thermal conductivity (k) are calculated
dimensional domain for simulation of pair bubble rising in stag- using a weighted average:
nant liquid and simulation of nucleating boiling in subcooled liquid.
Pan et al. (Pan et al., 2012) numerically investigated the behavior of y ¼ aL yL þ ð1:0  aL ÞyG ; y2½r; m; k (2)
condensing single vapor bubble in subcooled boiling flow within
The global continuity equation of two phase flow is given by:
two different vertical rectangular channels using the Volume of
Fluid (VOF) multiphase flow model. Zeng et al. (2015) investigated
the bubble condensation process for different initial bubble sizes v  !
ðrÞ þ V$ r U ¼ 0 (3)
and subcooled temperatures using the couple of level set method vt
(LS) and volume of fluid (VOF) in 3D domain. In two phase flow, the interface moves with flow. Thus, the
In present study, the vapor bubble condensing process in transport equation should be solved for VOF function to keep
quiescent water is simulated using volume of fluid method in interface. This transport equation is derived from global continuity
OpenFOAM CFD package (H. G. Weller et al., 1998). The OpenFOAM equation by substitution of density from equation (2) and defined
library allows implementations of fields, equations and operator as:
discretization using high level Cþþ(Weller et al., 1998). In the
previous study VOF method in interFoam solver has been used to !
vaL ! ! r V$ U
simulate single bubble rising (Samkhaniani et al., 2012). Deshpande þ U $VaL þ aL V$ U ¼  G (4)
et al. (Deshpande et al., 2012) evaluate the performance of this vt ðrL  rG Þ
solver using a variety of verification and validation test cases and !
found out interFoam is generally comparable with the recent where V$ U is calculated from local continuity equations for each
algebraic VOF algorithms. In order to simulate phase change pro- phase:
cess in bubble condensation phenomena, interFoam solver in  
OpenFOAM220 is extended to solve energy equation and calculate ! 000 1 1
V$ U ¼ m_  (5)
mass flux. rG rL
000
where m_ (kg/m3s) is volumetric transferred mass rate. Momentum
equations are given by:
2. Mathematical model
 !   !T
v rU  !! !
2.1. Governing equation þ V$ r U U  V$ m V U þ V U
vt
!
¼ VP þ r g þ skVaL (6)
In present study, the two phase flow is treated as incompressible
and immiscible Newtonian fluid. Interface between two phases are Last term in right hand of equation (6) indicates surface tension
resolved using color function volume of fluid method (CF-VOF). force between two phases. s is surface tension and k is interface
Following discontinuous scalar function is applied for interface curvature. The surface tension is accounted by Continuum Surface
tracking in fixed Eulerian grids. This scalar function is the ratio of Force model (CSF) without the density averaging proposed by
one fluid volume to the volume of cell and defined as: (Brackbill et al., 1992). Curvature is defined as:
122 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131

0 1 equation. It is assumed that the vapor phase in saturation tem-

perature has an ideal gas behavior.
@ VaeL A
k ¼ V$ (7)  

VaeL Psat;1 MHLG 1 1
ln ¼  (12)
Psat;0 R Tsat;1 Tsat;0
where aeL is calculated from the VOF function a by smoothing it over
a finite region around the interface. In the VOF method, the fluid
interface sharply changes over a thin region. This abrupt change of
2.2. Numerical details
the VOF function creates errors in calculating the normal vectors
and the curvature of the interface, which are used to evaluate the
The present solver is implemented using OpenFOAM finite
interfacial forces. These errors induce non-physical spurious cur-
volume library (H. G. Weller et al., 1998). It is based on interFoam
rents in the interfacial region (Klostermann et al., 2013; Scardovelli
solver. This solver supports simulation of incompressible two phase
and Zaleski, 1999). The spurious currents create extra heat con-
adiabatic flows without phase change using algebraic VOF. In order
vection around interface which increases local mass transfer. An
to simulate the subcooled flow boiling, the base solver is extended
easy way to suppress these artifacts is to compute the interface
by adding a thermal energy transport equation and modifying the
curvature from a smoothed VOF function. In this study the
base solver (the volume fraction transport equation and mo-
smoother proposed by Lafaurie et al.(Lafaurie et al., 1994) is
mentum equations) with phase-change source terms. The overall
applied:
solver algorithm can be described in following steps:
Pn
f ¼1 aLf Sf
fp ¼ Pn
a (8) 1. Define vector and scalar fields for the multiphase flow problem
!
f ¼1 Sf including: U , P, T and aL.Note that pressure used in the VOF
solver is the dynamic pressure Prgh where Prgh ¼ Prgh, where h
where Sf is the magnitude of face area, the subscript P denotes the is the liquid height.Prgh is used to avoid any sudden changes in
cell index and f denotes the face index. The interpolated value (aLf) the pressure at the boundaries for hydrostatic problems
at the face centre is calculated using linear interpolation. The (Rusche, 2003)
application of this filter can be repeated m times to get a smoothed 2. Start the time loop and solve the volume fraction advection.
field. It should be stressed that smoothing tends to level out high
curvature regions and should therefore be applied only up to the The base solver uses algebraic VOF which means no interface
level that is strictly necessary to sufficiently suppress parasitic reconstruction method is employed to locate the exact position of
currents. Hoang and et al. (Hoang et al., 2013) found that the interface in each cell. Therefore, the interface is diffuse among two
magnitude of parasitic current decreases up to one order from or three cells. In order to limits the smearing of the interface
m ¼ 0 to m ¼ 2 and only a slight further decrease was observed for because of the compensation of the diffusive fluxes. The extra
m > 2. Therefore, in present study m ¼ 2 is employed in all !
divergence term V$ðaL ð1  aL Þ Uc Þ is added to equation (4), which
simulations. contributes only in the region of the interface (0 < aL < 1.0) and
Energy equation is given by: diminishes in liquid or vapor phases.
    
v    000     vaL ! ! 1 1 1
rCp T þ V$ rCp UT  V$ðkVTÞ ¼ m_ HLG þ CpL  CpG TSat
000

vt þ U $VaL þ V$ aL ð1  aL Þ Uc ¼ m_  aL 
vt rL rL rG
(9) (13)
The term in the right hand of equation is due to phase change. !
HLG is the latent heat coefficient. The transferred mass flux should Uc is compressive velocity. It is calculated in the normal direction
be computed using an appropriate mass transfer model. In this to the interface to avoid any dispersion. Moreover, a compressive
study phase change model proposed by Tanasawa (Tanasawa, 1991) factor (Ca) is used to increase compression:
is employed. Mass flux (kg/m2s) at the liquid vapor interface is
! VaL
calculated as: Uc ¼ minfCa jUj; maxðjUjÞg (14)
jVaL j
rffiffiffiffiffiffiffiffiffi
00 2g M rG HLG ðT  Tsat Þ In present study the compression factor Ca ¼ 1.0 is considered.
m_ ¼ (10)
g1 2pR T
3=2 The coefficient Ca controls the weight of the compression flux and
sat
should be usually in the range of unity (1.0 < Ca < 4.0) (Berberovic
where M is molar mass of fluid, R ¼ 8.314 J/mol K is the universal et al., 2009; Samkhaniani et al., 2013).
gas constant, TSat(P) is local saturation temperature, and g is the The volume fraction advection equation is solved using the
fraction of molecules transferred from one phase to the other multidimensional universal limiter with explicit solution (MULES)
during phase change. Marek and Straub (Marek and Straub, 2001) method which is based on the method of flux corrected transport
determined g based on published data. They recommended (FCT) (Zalesak, 1979) where an additional limiter is used to cutoff
g ¼ 0.1  1 for dynamically renewing water surfaces such as jets or the face-fluxes at the critical values. This solver is included in
moving films, and g < 0.1 for stagnant surfaces. More details on the OpenFOAM library, and performs conservative transport equation
different phase change model can be found in (H. Lee et al., 2015). of hyperbolic transport equations with defined bounds (0 and 1
For equations (5) and (9), the volumetric mass source term in kg/ foraL).
m3s is needed which is determined from:
3. Update the fluid physical properties and the fluxes using the
000 00 volume fraction function aL using equation (2).
m_ ¼ m_ VaL (11)
4. Solve the Navier Stokes system of equations for velocity and
The variation of saturation temperature based on local pressure pressure using the Pressure Implicit with Splitting of Operators
P is calculated using a simplified version of ClausiuseClapeyron (PISO) (Issa et al., 1986). In the present study, three pressure
 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 123

Table 1
Discretisation schemes.

Term Discretisation scheme Method

v ! v ! Euler the first order bounded implicit scheme

vt ðr U Þ; vt ðr U TÞ
!! vanLeerV Similar to VanLeer scheme (Van Leer, 1974) modified for vector field
V$ðr U U Þ
! ! vanLeer See (Van Leer, 1974)
V$ð U aL Þ; V$ðr U TÞ
! InterfaceCompression See (H. Weller, 2008)
V$ð Uc aL ð1  aL ÞÞ
Vc* Linear central difference schemes
** corrected surface normal gradient with correction on non-orthogonal meshes (Jasak, 1996)
V1f c
V$(c1Vc2) Linear corrected face values (c1) approximated by central difference scheme, and the resulting
surface normal gradient is calculated using central difference schemes with
non-orthogonal correction

Term Interpolation scheme Method

cf Linear Default interpolation schemes for getting face values from cell values

correction steps were used and ensured that the pressure re-
v 000    
sidual remained always below 1010. ðTÞ þ V$ðUTÞ  V$ðDk VTÞ ¼ Dc m_ HLG þ CpL  CpG Tsat
vt
First, the matrix equation for the momentum equation is (16)
formed. Then the inner pressureevelocity correction process is
initiated. In PISO, an intermediate velocity field is first obtained, where Dk and Dc are defined as:
and the cell-face volume fluxes (f) are evaluated and corrected for
1
gravitational forces, the continuum surface-tension force, and Dc ¼ (17)
rL CL aL þ rG CG ð1  aL Þ
boundary conditions. The pressureePoisson equation is then
formed and solved. Following the approach of (Weller et al., 1998),
the coefficients of the pressure equation are obtained from the kL aL þ kG ð1  aL Þ
Dk ¼ (18)
diagonal entries of the momentum matrix equation (AD). For pre- rL CL aL þ rG CG ð1  aL Þ
sent condensing flows, the pressure equation would be:

6. Move to the next time step (starting from 2).

For the stability of the solution procedure, the calculations are

performed using a self-adapting time step based on the user
    defined maximum Courant number (Comax) and maximum time
1 000 1 1
V$ VPrgh ¼ V$f  m_  (15) step (Dtmax) (Berberovi c et al., 2009). Here, Comax ¼ 0.1 is employed
AD rG rL
for all simulations.
where the last term in the right hand of equation (15) is added to The governing equations are discretized based on a finite vol-
consider phase change process. The pressure correction equation is ume method (FVM) formulation. The discretization is performed on
solved in present study using a conjugate gradient iterative (PCG) a structured grid. All the variables are stored at the cell centers,
solver preconditioned with a geometric algebraic multi grid where a collocated grid arrangement is used. In order to avoid a
method (GAMG). checker boarding effect in the momentum equation, the Rhie-Chow
momentum interpolation (Rhie and Chow, 1983) is employed.
5. Solve energy Equation for temperature The applied discretisation schemes and the parameters of the
numerical model are summarized in Table 1. For convenience, the
The numerical interfacial flow computations become more corresponding terminology of OpenFOAM is given.
challenging when the imbalance of material properties between
the phases increases (Scardovelli and Zaleski, 1999). In order to 2.3. Validation of phase change model
reduce material properties imbalance through interface and in-
crease numerical robustness, equation (9) is redefined as: The one dimensional Stefan problem was introduced for

Fig. 1. Schematic of Stefan problem and boundary conditions, TSat ¼ 380.26[K], DTsub ¼ 10[K], CpG ¼ CpL(psat).
124 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131

Table 2
Water thermo physical properties at saturation temperature (Psat ¼ 0.130 MPa,
Tsat ¼ 380.26).

Dimension Liquid Vapor

Density Kg/m3 953.13 0.75453

Viscosity Pa.s 2.6E-04 1.25E-05
Thermal conductivity W/mK 0.68106 0.025905
Specific Heat kJ/kgK 4.224 2.110
Latent Heat kJ/kg 2237.41
Surface Tension N/m 0.05753

Fig. 3. Schematic of bubble condensation in subcooled boiling, initial and boundary

conditions.

Fig. 2. Stefan problem, comparison between analytical solution and CFD model,
Psat ¼ 0.130 MPa and DT ¼ 10 K.

Table 3
Stefan problem convergence.

64 128 256 512 Exact

d(t ¼ 100)(mm) 7.97001 7.9691 7.9689 7.9688 7.95

E(mms) 1.17 0.51 0.30 0.2 e

solidification at first (Alexiades, 1992) but it became a well-known

bench mark for boiling simulation (Guo et al., 2011; Hardt and
Wondra, 2008; Welch and Wilson, 2000). However, it can be con-
ducted as a validation case for any one dimensional phase change
phenomenon. In present study, Stefan problem is solved to validate Fig. 4. Bubble shape at t ¼ 1 ms for different grids, Psat ¼ 0.130[Mpa] and DTsub ¼ 25[K].
condensation. The schematic of Stefan problem for condensation
process is illustrated in Fig. 1. Heat is transferred by conduction
from saturated vapor phase to liquid phase and it is rejected
through subcooled wall. The vapor phase condensation leads to a
motion of the interface to the right. It is assumed that during the
condensation process the interface stays flat. Hence, the problem
can be regarded as one-dimensional. The analytical solution of this
problem is given by (Welch and Wilson, 2000):

pffiffiffiffiffiffiffi
xðtÞ ¼ 2h dL t (19)

where x is the interface position from cold wall, dL is the liquid

thermal diffusivity and h is determined from:

  cpL ðTsat  Twall Þ

h exp h2 erf ðhÞ ¼ pffiffiffi (20)
pHLG

where HLG and cpL is latent heat and liquid specific heat, respec- Fig. 5. Bubble diameter history, compare between present numerical simulations on
tively and erf is the error function. different grids and experimental data (Kamei and Hirata, 1990a,b), D0 ¼ 1.008 mm,
A quasi 1D computational domain with only one grid cell in the Psat ¼ 0.130[Mpa] and DTsub ¼ 25[K].
direction of translational invariance is considered. The
 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 125

Fig. 6. Bubble shape, comparison between present numerical result and experimental data (Kamei and Hirata, 1990a,b), D0 ¼ 1.008 mm, Psat ¼ 0.130[MPa] and DTsub ¼ 25[K].

liquidevapor thermo physical properties for water at saturation 3. Result and discussion
pressure 0.130 MPa is chosen and displayed at Table 2. In order to
ensure that the coefficient of mass flux in energy equation is con- 3.1. Problem description and validation
stant in CFD model during phase change process, liquid and vapor
phases specific heat are assumed equal (CpG ¼ CpL(Psat)). No slip In present study, the rising of a single vapor bubble in quiescent
boundary condition is employed for velocity boundary condition at subcooled water is simulated. The geometry of the considered
the walls. Temperature of cold wall is 10 less than saturation problem is illustrated in Fig. 3. The 2-dimensional space domain is
temperature. set as 2D0  4D0 where D0 is the initial diameter of vapor bubble.
Since there is no sharp interface at color function VOF, iso- The bubble is located in the position of (D0,D0) at the beginning of
contour aL ¼ 0.5 is applied as interface. Result of the comparison simulation. It should be noted that the solver is capable of 3D
between exact solution and CFD model for Stefan problem is simulation. In order to simulate 2D simulation using OpenFOAM,
depicted in Fig. 2. special boundary condition called empty is applied in front and
The integrated simulation error can be estimated as the film back faces.
thickness error (jdsimdanj) summed over time stepsi, weighted by
Dt ¼ 1E05(see Table 2).

3.2. Mesh resolution dependence

A mesh resolution analysis is performed on 4 different grids to

figure out the appropriate grid size for simulation of bubble
X
E¼ jdsim  dan jDt (21) condensation in Fig. 4. These grids are 50  100, 75  150,
i 100  200 and 150  300 which respectively represent 25, 37, 50
and 75 cells across the initial diameter. It is obviously seen that
There is an excellent agreement between present numerical bubble shape becomes converged with mesh 100  200 and
result and exact solution based on illustrated result in Fig. 2 and 150  300. The variation among two lateral grids is negligible.
Table 3. Therefore the bubble is resolved using 100  200 cells in this study.

Fig. 7. Comparison of bubble diameter history, fixed saturation temperature vs vari- Fig. 8. The influence of mass transfer coefficient (g) in Tanasawa model on CFD result.
able saturation temperature based on local pressure, experiment data is provided by Experimental data of bubble condensation at D0 ¼ 1.008 mm, Psat ¼ 0.130[MPa] and
(Kamei and Hirata, 1990a,b). DTsub ¼ 25[K] is from (Kamei and Hirata, 1990a,b).
126 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131

water is a complex phenomenon. The mass transfer flux is pro-

portional with the difference of vapor local saturation temperature
with bulk liquid temperature. Local saturation temperature is a
function of thermodynamic pressure which is non-uniform in
quiescent water due to surface tension and gravity. Surface tension
causes higher pressure in vapor bubble than surrounding liquid and
gravity increases liquid pressure proportional to water column
depth. Furthermore, various forces act on bubble such as capillary
and buoyancy resulting different regimes of bubble.
From numerical point, real time simulation of two phase flow
with phase change is a challenging task. For instance, the difference
in the values of thermo physical properties across the narrow
interface is high; the ratio of density rrL in present simulation is
G
about 1263. Moreover, for precise computation of mass flux in
phase change process, an accurate mass transfer model is required.
Some mass transfer models such as Tanasawa (Tanasawa, 1991) and
Fig. 9. The influence of mass transfer coefficient (g) in Tanasawa model on CFD result. Lee (Lee, 1980) are semi-experimental correlation and must be
.Experimental data of bubble condensation at D0 ¼ 0.95 mm, Psat ¼ 0.101[MPa] and tuned up for numerical simulation using experimental data.
DTsub ¼ 12.8[K] is from (Kamei and Hirata, 1990a,b).
Furthermore, accurate interface tracking method is necessary to
compute surface tension force on vapor bubble interface correctly.
3.3. Validation of bubble condensation The poor estimation of interface curvature causes unphysical
spurious current near interface which may impair the mass flux
In order to validate present numerical simulation, the shape and calculation. Unfortunately, recent published articles (Bahreini et al.,
the bubble diameter history is compared with published experi- 2015) (Pan et al., 2012; Zeng et al., 2015) in this field have not
mental results (Kamei and Hirata, 1990a,b) under saturation pres- addressed these issues. In this section the importance of different
sure 0.130 MPa. Thermo physical properties are displayed in aspects of numerical parameters on final CFD result is discussed.
Table 2. The comparison of bubble diameter history is shown in
Fig. 5 and the bubble shape at Fig. 6. Bubble diameter in 2D 3.4.1. Local saturation temperature
simulation is the diameter of equivalent circle. It can be seen that To study the influence of local saturation temperature on the
numerical result is in reasonable good agreement with the exper- numerical result, the test case similar to previous section is carried
imental results. out using constant saturation temperature for vapor; the result is
compared with experimental data in Fig. 7. It shows when
3.4. Discussion on the single bubble condensation simulation Tsat ¼ const is assumed, bubble diameter reduces almost linearly;
because the temperature difference between saturated vapor
Even the rising of a single vapor bubble in quiescent subcooled bubble and subcooled water is constant during simulation. The area

Fig. 10. Bubble shape sequences in condensation process started from different initial shape, DT ¼ 25k, Psat ¼ 0.130 MPa, D0 ¼ 8 mm.
 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 127

Fig. 11. Bubble diameter history, DT ¼ 25k, Psat ¼ 0.130 MPa, D0 ¼ 8 mm. Fig. 13. The effect of spurious current magnitude on bubble life history, s ¼ 0.05753 N/
m.

to volume ratio is becoming bigger in small bubble diameter;

therefore, bubble is condensing faster in last stage. However, while hitting the interface at the lower saturation pressure and subcooled
a bubble is getting smaller in condensation process, the bubble temperature. Hence, the value of mass transfer coefficient should
inside pressure is increasing. Hence, if saturation temperature is be re-evaluated in each simulation by comparing it with experi-
computed as a function of thermodynamic pressure, it is slightly mental data. However, in this study the value of g ¼ 0.9 is employed
increased and accelerates bubble condensation. In the present for the rest of simulations, it is assumed the dependence of g to
study, the small variation of local saturation temperature due to subcooled temperature is negligible.
pressure is calculated using ClausiuseClapeyron relation (equation
(12)). 3.4.3. Effect of initial shape
It is assumed in CFD model that bubble initial shape is circle. It is
3.4.2. Mass transfer coefficient often true for small bubbles, where surface tension is predominant
The importance of variable g in Tanasawa phase change model is and deviation of bubble shape from circularity is small. But the
investigated by comparison of bubble diameter history corre- deviation can be big in large bubbles and it may affect the nu-
sponding different value of g in Fig. 8. g ¼ 0.9 is suitable for present merical result and lead to the wrong interpretation while
simulation. When a vapor molecule collides on interface, it may comparing with experimental data. In this section, 3 different
reflect from interface or convert to liquid molecule. g ¼ 0.9 means initial shapes of bubble (circle, horizontal ellipse and vertical el-
only 90% of vapor molecules are converting to the liquid and the lipse) are compared with each other in Fig. 10 and Fig. 11. These
remain is reflecting. The result shows bubble diameter history bubbles have the same area. The elliptic shape is defined as:
derived from numerical data is highly sensitive to mass transfer
coefficient. Another drawback of Tanasawa phase change model is ðx  x0 Þ2 ðy  y0 Þ2 D2
2
þ 2
¼ 0 (22)
that mass transfer coefficient is a function of temperature and m n 4nm
pressure. As seen in Fig. 9, for the vapor bubble condensation at
where for horizontal ellipse (m ¼ 1, n ¼ 1.5) and for vertical
saturation pressure 0.101 MPa and subcooled temperature 12.8 K,
(m ¼ 1.5, n ¼ 1) is considered.
the suitable value of g is around 0.4. The reason is that smaller
The CFD result shows bubble shape sequences obviously vary
fraction of vapor molecules convert to liquid molecules while
due to different initial shape, but the bubble life time history is not
sensitive to initial shape.

3.4.4. Spurious current

The interfacial tension force has been applied to Eulerian grids
using different approaches a) continuous surface stress (CSS)
method (Gueyffier et al., 1999), b) continuous surface force (CSF)
(Brackbill et al., 1992) method and finally (c) sharp surface force
(SSF) (Renardy and Renardy, 2002). In present study, CFS model is
employed. However, an often reported problem of all these
implementations (Brackbill et al., 1992; Gueyffier et al., 1999;
Renardy and Renardy, 2002) is the existence of spurious currents
in the flow field of the numerical simulations in the vicinity of the
interface. In order to investigate the influences of spurious current
on CFD results, the condensation of a single vapor bubble in zero
gravity in quiescent subcooled water is simulated. In the absence of
external force, flow must be still and motionless. Heat transfers
from saturation vapor bubble to subcooled water by conduction
heat mode and moves the interface. Spurious currents generate
Fig. 12. Spurious currents in the vicinity of interface vapor bubble in zero gravity small vortex near interface as seen in Fig. 12. These small currents
condition, DT ¼ 25k, Psat ¼ 0.130 MPa, Tsat ¼ const. increase the heat transfer coefficient and lead to extra mass flux
128 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131

Fig. 14. Contour of mass flux, temperature and relative velocity of bubble condensation at Psat ¼ 0.130 MPa, DT ¼ 20 K, D ¼ 4 mm.

across the interface. In this series of simulation, saturation tem- spurious current on the numerical simulations of the single vapor
perature is assumed constant to omit the effect of pressure inside bubble condensation.
bubble on life time history. As seen in the Fig. 13, bubble life time Here, a smoothing filter method (equation (8)) is applied which
reduces by the increase in the surface tension value. Because the improves interface curvature calculation and reduces the spurious
magnitude of largest spurious currents around a bubble is linear current magnitude in one order. The filter makes remedy the
proportional with the magnitude of surface tension (Lafaurie et al., problem but it cannot totally fade it out.
1994). This test shows the significant influences of unphysical

Fig. 15. Vapor bubble condensation at Psat ¼ 0.130 MPa, DT ¼ 20 K, D ¼ 4 mm.

 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 129

produces higher mass flux. In the back face of bubble, relative ve-
locity is upward. It is higher in the middle part of back face. A thin
thermal boundary layer is generated in the behind of condensing
vapor bubble, thus, subcooled temperature reduces leading to
lower mass flux at back face.

3.5.2. Multiple bubble condensation

The Fig. 15 shows the performance of multiple vapor bubbles
condensation on each single bubble. Three simulations are con-
ducted in this section. In case A, two similar vapor bubbles are
placed in computational domain at the beginning of simulation.
The distance between bubbles centre is 1.1D0. In cases B and C the
condensation of each lateral bubble are considered separately.
Since thermodynamic pressure increases by water column
depth, local saturation temperature is higher in deeper side. Hence,
lower bubble experiences higher mass flux and condensates faster
Fig. 16. Bubble life time for different sizes of bubble at PSat ¼ 0.130 MPa, DT is sub- (compare Fig. 15-Case B and 15-Case C). The comparison of Fig. 15-
cooled temperature. Case A and 15-Case B indicates when a bubble moves afterward
another bubble; the lower bubble life time becomes longer. Despite
of single vapor bubble condensation, in this situation, the lower
3.5. Bubble condensation bubble is placed in the thermal boundary layer of above bubble.
Thus, it is in exposure of less subcooled water which leading to
3.5.1. Single bubble condensation smaller mass flux and longer life time.
In the condensation process of vapor bubble, mass flux is higher
in frontal face of bubble as illustrated in Fig. 14, because it is in
exposure of lower subcooled temperature and higher relative ve- 3.5.3. Bubble life time
locity. Relative velocity is defined as: In a single vapor bubble condensation, bubble life time is
dependent to bubble size and subcooled temperature. Fig. 16 in-
dicates bubble size increases life time almost linearly. Fig. 17 shows
! ! !
U r ¼ U  Ub (23) the bubble shape sequences corresponding to different subcooled
temperature value. It is found clearly that increase in the subcooled
!
where Ub is bubble velocity computed from bubble mass center temperature decreases bubble life time. However, there is no linear
displacement. The relative velocity is downward in the front of relationship between bubble life time and subcooled temperature.
bubble. In lateral sides of bubble front face, relative velocity reaches Present numerical simulations in Fig. 18 indicates that the bubble
maximum value. It increases heat convection coefficient and life time has a rapidly decrease in the low subcooled degree region,

Fig. 17. Bubble shape sequences for condensation process atPsat ¼ 0.130 MPa, D0 ¼ 4 mm, DT is subcooled temperature.
130 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131

life time is almost insensitive to subcooled temperature

changes.
4 In a single vapor bubble condensation in quiescent water, mass
flux is higher in the front face of bubble than back face.
5 Thermal boundary layer in bubble swarm motion prolongs
lower bubbles condensation.

References

Alexiades, V., 1992. Mathematical Modeling of Melting and Freezing Processes. CRC
Press.
Bahreini, M., Ramiar, A., Ranjbar, A.A., 2015. Numerical simulation of bubble
behavior in subcooled flow boiling under velocity and temperature gradient.
Nucl. Eng. Des. 293, 238e248.
Berberovi c, E., Hinsberg, N.P.V., Jakirli
c, S., Roisman, I.V., Tropea, C., 2009. Drop
impact onto a liquid layer of finite thickness: dynamics of the cavity evolution.
Am. Phys. Soc. 79, 036306e036315.
Fig. 18. Bubble life time vs. subcooled temperature at Psat ¼ 0.130 MPa. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling
surface tension. J. Comput. Phys. 100, 335e354.
Brucker, G., Sparrow, E., 1977. Direct contact condensation of steam bubbles in water
at high pressure. Int. J. Heat Mass Transf. 20, 371e381.
but tends to be steady as the subcooled degree is large enough. Chen, R., Tian, W., Su, G., Qiu, S., Ishiwatari, Y., Oka, Y., 2010. Numerical investigation
Zeng et al. (Zeng et al., 2015) have been noted that once subcooled on bubble dynamics during flow boiling using moving particle semi-implicit
method. Nucl. Eng. Des. 240, 3830e3840.
degree changes, the fluid thermodynamic properties changes as
Chen, R., Tian, W., Su, G., Qiu, S., Ishiwatari, Y., Oka, Y., 2011. Numerical investigation
well, which thus leads to a non-linear relation between bubble on coalescence of bubble pairs rising in a stagnant liquid. Chem. Eng. Sci. 66,
lifetime and the subcooled degree. However, it seems this behavior 5055e5063.
returns to temporal scale of problem. The total heat capacity can be Chen, Y., Mayinger, F., 1992. Measurement of heat transfer at the phase interface of
condensing bubbles. Int. J. Multiph. Flow 18, 877e890.
transferred using conduction and convection mode in a definite Deshpande, S.S., Anumolu, L., Trujillo, M.F., 2012. Evaluating the performance of the
portion of time is limited. Hence, in low subcooled temperature, two-phase flow solver interFoam. Comput. Sci. Discov. 5, 014016.
increase in subcooled temperature increases transferred heat Gueyffier, D., Li, J., Nadim, A., Scardovelli, R., Zaleski, S., 1999. Volume-of-fluid
interface tracking with smoothed surface stress methods for three-dimensional
across the interface but in large subcooled temperature, it has no flows. J. Comput. Phys. 152, 423e456.
effect either on transferred heat or mass flux. Thus, bubble life time Guo, D.Z., Sun, D.L., Li, Z.Y., Tao, W.Q., 2011. Phase change heat transfer simulation
in large subcooled region become steady. This numerical result is for boiling bubbles arising from a vapor film by the VOSET method. Numer.
Heat. Transf. 59, 857e881.
consistent with Sudhoff correlations (Sudhoff et al., 1982) which is a Harada, T., Nagakura, H., Okawa, T., 2010. Dependence of bubble behavior in sub-
semi-empirical correlation on the basis of literature reviews. cooled boiling on surface wettability. Nucl. Eng. Des. 240, 3949e3955.
Hardt, S., Wondra, F., 2008. Evaporation model for interfacial flows based on a
continuum-field representation of the source terms. J. Comput. Phys. 227,
5871e5895.
Hoang, D.A., van Steijn, V., Portela, L.M., Kreutzer, M.T., Kleijn, C.R., 2013. Benchmark
4. Concluding remarks numerical simulations of segmented two-phase flows in microchannels using
the volume of fluid method. Comput. Fluids 86, 28e36.
Issa, R.I., Gosman, A., Watkins, A., 1986. The computation of compressible and
In present study, condensation of vapor bubble in subcooled incompressible recirculating flows by a non-iterative implicit scheme.
water is modeled using CF-VOF method in OpenFOAM solver J. Comput. Phys. 62, 66e82.
(interFoam). In order to simulate phase change process, energy Jasak, H., 1996. Error Analysis and Estimation for the Finite Volume Method with
Applications to Fluid Flows. University of London.
equation, Tanasawa mass transfer model and appropriate source
Kamei, S., Hirata, M., 1990a. Condensing phenomena of a single vapor bubble into
terms are added to base solver. Implemented code is validated with subcooled water. Exp. Heat Transf. Int. J. 3, 173e182.
analytical solution of Stefan problem, and then applied for simu- Kamei, S., Hirata, M., 1990b. Study on condensation of a single vapor bubble into
subcooled water-Part 2; experimental analysis. Heat Trans. Jpn. Res. (USA) 19.
lation of a single and multiple vapor bubbles condensation. A single
Klostermann, J., Schaake, K., Schwarze, R., 2013. Numerical simulation of a single
vapor bubble shape and life time history in condensation is verified rising bubble by VOF with surface compression. Int. J. Numer. Methods Fluids
with the experimental data (Kamei and Hirata, 1990a,b). The limi- 71, 960e982.
tations of mass transfer model and numerical method are discussed Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., Zanetti, G., 1994. Modelling
merging and fragmentation in multiphase flows with SURFER. J. Comput. Phys.
exclusively. Furthermore, some promising remarks are obtained: 113, 134e147.
Lee, H., Kharangate, C.R., Mascarenhas, N., Park, I., Mudawar, I., 2015. Experimental
1 The main limitation of bubble condensation simulation using and computational investigation of vertical downflow condensation. Int. J. Heat
Mass Transf. 85, 865e879.
interface tracking method such as VOF is spurious current. It Lee, W.H., 1980. A pressure iteration scheme for two-phase flow modeling. Multiph.
impairs numerical result by false extra heat convection near Transp. Fundam. React. Saf. Appl. 1, 407e431.
interface leading to increase in mass flux. Lucas, D., Prasser, H.-M., 2007. Steam bubble condensation in sub-cooled water in
case of co-current vertical pipe flow. Nucl. Eng. Des. 237, 497e508.
2 Despite the extensive use of Tanasawa mass transfer model for Magnini, M., Pulvirenti, B., Thome, J.R., 2013. Numerical investigation of hydrody-
simulation of phase change (Hardt and Wondra, 2008; Magnini namics and heat transfer of elongated bubbles during flow boiling in a
et al., 2013; Ranjan et al., 2011), this model is not recommended microchannel. Int. J. Heat Mass Transf. 59, 451e471.
Marek, R., Straub, J., 2001. Analysis of the evaporation coefficient and the conden-
for simulation of bubble condensation in the absence of exper- sation coefficient of water. Int. J. Heat Mass Transf. 44, 39e53.
imental data. Our study shows CFD results are highly sensitive to Pan, L.-m., Tan, Z.-w., Chen, D.-q., Xue, L.-C., 2012. Numerical investigation of vapor
mass transfer coefficient g and must be tuned up with experi- bubble condensation characteristics of subcooled flow boiling in vertical rect-
angular channel. Nucl. Eng. Des. 248, 126e136.
mental data for different working conditions.
Ranjan, R., Murthy, J.Y., Garimella, S.V., 2011. A microscale model for thin-film
3 The bubble life time is linearly dependent to bubble size and not evaporation in capillary wick structures. Int. J. Heat Mass Transf. 54, 169e179.
sensitive of initial bubble shape. It has a non linear dependence Renardy, Y., Renardy, M., 2002. PROST: a parabolic reconstruction of surface tension
to subcooled temperature. The influence of subcooled temper- for the volume-of-fluid method. J. Comput. Phys. 183, 400e421.
Rhie, C., Chow, W., 1983. Numerical study of the turbulent flow past an airfoil with
ature variation in bubble life time is higher in low subcooled trailing edge separation. AIAA J. 21, 1525e1532.
temperature. In enough large subcooled temperature, bubble Rusche, H., 2003. Computational Fluid Dynamics of Dispersed Two-phase Flows at
 N. Samkhaniani, M.R. Ansari / Progress in Nuclear Energy 89 (2016) 120e131 131

High Phase Fractions. Imperial College London (University of London). computation of thermally controlled steam bubble condensation using Moving
Samkhaniani, N., Ajami, A., Kayhani, M.H., Dari, A.S., 2012. Direct numerical simu- Particle Semi-implicit (MPS) method. Ann. Nucl. Energy 37, 5e15.
lation of single bubble rising in viscous stagnant liquid. In: International Con- Van Leer, B., 1974. Towards the ultimate conservative difference scheme. II.
ference on Merchanical, Automobile and Robotics Engineering (ICMAR'2012). Monotonicity and conservation combined in a second-order scheme. J. Comput.
Samkhaniani, N., Gharehbaghi, A., Ahmadi, Z., 2013. Numerical simulation of re- Phys. 14, 361e370.
action injection molding with polyurethane foam. J. Cell. Plastics 49, 405e421. Welch, S.W.J., Wilson, J., 2000. A volume of fluid based method for fluid flows with
Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and phase change. J. Comput. Phys. 160, 662e682.
interfacial flow. Annu. Rev. Fluid Mech. 31, 567e603. Weller, H., 2008. A New Approach to VOF-based Interface Capturing Methods for
Sideman, S., Hirsch, G., 1965. Direct contact heat transfer with change of phase: Incompressible and Compressible Flow. OpenCFD Ltd. Report TR/HGW/04.
condensation of single vapor bubbles in an immiscible liquid medium. Pre- Weller, H.G., TaboraI, G., Jasak, H., Fureby, C., 1998. A tensorial approach to
liminary studies. AIChE J. 11, 1019e1025. computational continuum mechanics using object-oriented techniques. Com-
Sudhoff, B., Plischke, M., Weinspach, P., 1982. Direct contact heat transfer with put. Phys. 12, 620e632.
change of phase-condensation or evaporation of a drobble. Ger. Chem. Eng. 5, Zalesak, S.T., 1979. Fully multidimensional flux-corrected transport algorithms for
24e43. fluids. J. Comput. Phys. 31, 335e362.
Tanasawa, I., 1991. Advances in condensation heat transfer. Adv. Heat Transf. 21, Zeng, Q., Cai, J., Yin, H., Yang, X., Watanabe, T., 2015. Numerical simulation of single
55e139. bubble condensation in subcooled flow using OpenFOAM. Prog. Nucl. Energy
Tian, W., Ishiwatari, Y., Ikejiri, S., Yamakawa, M., Oka, Y., 2010. Numerical 83, 336e346.