7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Numerical study of single bubble dynamics with dynamic contact angle using a
conservative level set method
Muhammad Sajid and Rachid Bennacer
L2MGC, Universit6 de CergyPontoise,
33 boulevard du Port, CergyPontoise, 95011, France
Muhamaad.Sajid@ucergy.fr and Rachid.Bennacer@ucergy.fr
Keywords: Dynamic contact angle, Bubble, Boiling, Level set method
Abstract
During nucleate pool boiling bubbles nucleate at the heated surface and stay attached to it through out their growth period. The
shape of a bubble growing on a heated surface depends on the material properties of the liquid surrounding the bubble, the
vapor inside it and the surface on which it is resting. This dependency manifests as the contact angle which has been described
as a function of the interfacial tensions using Youngs equation for a static equilibrium contact angle. However, in real systems
a spectrum of contact angles are observed during the bubble growth process. The variation of the contact angles influences
bubble departure diameters and growth rates. Therefore, for the numerical simulation of bubble dynamics, an appropriate
contact angle model is essential for a correct prediction of the bubble formation. In this study a forcebalance based approach
is to modelling dynamic contact angle is applied. A conservative level set method coupled with the Navier Stokes equations is
used to represent the liquid vapor interface. The results are compared with two approaches, the static mean contact angle
approach and contact line velocity based dynamic contact angle approach as well as experimental results available in literature.
The contribution of the dynamic contact angle in bubble dynamics and the associated heat transfer is highlighted.
Introduction
A thorough understanding of the process of boiling heat
transfer in pool boiling requires an investigation of the
thermodynamics of bubbles and the hydrodynamics of the
flow resulting from the departure of the bubbles from the
heater surface. Consequently for mechanistic modelling, it
is necessary to be able to predict data for a single bubble
nucleation which involves the mechanics of multiphase
fluids as well as the dynamics of phase transition. Contact
angles play an important role in the bubble growth cycle.
It has been experimentally observed that the contact angle at
the bubble base varies, ranging from the advancing contact
angle, OACA, to the receding contact angle, ORCA. The
experimental data typically obtained for variation of contact
angle, e.g. Dussan (1979) is shown in Figure1. The upper
limit of the spectrum is the advancing contact angle OACA,
which is the contact angle found at the advancing edge of
the liquid. The lower limit is the receding contact angle ORCA,
which is the contact angle found at the receding edge. In
accordance with these observations contact angles have
typically been modelled on the basis of the contact line
velocity.
However, the fundamental physics governing the contact
angle are based on the forces acting locally in the contact
line region. During bubble growth the contact line force,
which depends on the contact angle pulls the bubble
outwards while the buoyancy force, which depends on the
bubble volume, pulls it upwards. The balance between them
determines the dynamics of the growing bubble and
influences the contact angle, bubble growth rate, heat
transfer, bubble volume etc. The present numerical study is
carried out to evaluate an approach to modelling the
dynamic contact angle based on the balance of the local
forces, i.e. the contact line force and buoyancy force.
Nomenclature
Cp specific heat at constant pressure
d distance function (m)
g gravitational constant (ms1)
hfg Latent heat of evaporation (Jkg1)
th Mass transfer rate at interface (kgm2s1)
P pressure (Nm2)
u horizontal velocity (ms1)
v Vertical velocity (ms1)
Greek letters
PT coefficient of thermal expansion (K1)
F interface
6 Dirac delta
E half of interface thickness (m)
0 angle (o)
K curvature (m2)
k heat transfer coefficient (W m 2)
I viscosity (Pas)
p density (kgm 3)
o surface tension (Nm1)
0 level set function
Subscripts
ACA advancing contact angle
CA contact angle
CL
1
max
min
R
RCA
sat
v
w
contact line
liquid
maximum
minimum
characteristic
receding contact angle
saturation
vapor
wall
Literature Review
In the numerical simulation of bubble dynamics two types
of contact angles have been considered, static and dynamic.
Because of the large uncertainties associated with
determination of advancing and receding contact angles,
researchers have often used a constant value for the
apparent contact angle, i.e. static contact angle. Son et al
(1999) assumed a static contact angle at the bubble base
during their numerical simulation of a single bubble on a
horizontal surface during nucleate pool boiling. Their results
showed that the departing bubble became larger with
increase in contact angle. Abarajith and Dhir (2002) fixed
the contact angle throughout the bubble growth and
departure process while relating it to the Hamakar constant
which was found to change with surface wettability. A
careful study of contact angle hysteresis, using drops of
hydrocarbon chains was conducted by Lam et al.(2002).
Their work shows that as the advancing and receding
contact angles (OACA and ORCA) increase the hysteresis
between them decreases, such that at OACA and ORCA close to
90" liquids will have a constant contact angle. Thus a static
contact angle approach to simulating bubble dynamics is
more accurate at high values of contact angle.
0ACA
0RCA
UCL,max 0 UCL,max
Figure 1: Characteristic experimental results for contact
angle measurements against contact line velocity.
A generalized representation of the dynamic contact angle in
comparison to the contact line velocity is shown in Fig 1.
The curve is characteristic of experimental observations and
shows a quasilinear variation of contact angle between
limiting interface velocities. Here UCL,.m denotes the
limiting value of the contact line velocity at which a
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
measurement of the contact angle was made. Based on these
experimental observations the dynamic contact angle has
been typically modelled as a function of the interface speed.
Ramanujapu and Dhir (1999) derived a contact angle model
based on their experiments using a silicon wafer as the test
surface with micromachined cavities for nucleation, during
nucleate pool boiling. They calculated the interface velocity
and measured the bubble base diameter as a function of time.
They observed the dynamic contact angle as function of the
interface velocity and obtained a figure similar to Fig2.
They found that the contact angle varied during different
stages of bubble growth and was weakly dependent on the
interface velocity and concluded that the contact angle
could be determined based on the sign of the interface
velocity.
More recently, Mukherjee and Kandlikar (2007) studied the
effects of dynamic contact angles on the formation of single
bubbles, without considering the evaporation in the
microlayer, using two kinds of contact line models. In the
first approach they use a constant advancing and constant
receding contact angle depending on the sign of the contact
line velocity, uCL which was approximated as the rate of
change of bubble base diameter. In the second approach the
dynamic contact angle is a function of the rate of change of
bubble base diameter and thus contact line velocity based
on Fig1 with a smooth transition from receding contact
angle to advancing contact angle.
ACA if UCL< ax
CA AC ACA A *
C2 x 2
f CA L>U 2ax
Francois and Shyy (2003) used a dynamic contact angle
model based on contact line velocity to simulate droplet
spreading. Huang et al. (2004) computed the motion of a
twodimensional drop/bubble attached to a plane surface.
For the motion of contact line and the determination of the
dynamic contact angle, they used ad hoc models based on a
relationship between moving contact angle and contact line
speed. Smith et al. (2005) employed a contact line velocity
based model to study the contact line dynamics of an
interface between immiscible fluids.
One prominent aspect of the contact line velocity approach
is stickslip behaviour of the liquid vapour interface at the
contact point. The change in the contact angle affects the
balance of forces by altering the surface tension which acts
through the contact angle. The sudden disturbance in the
competing forces, the buoyancy force that pulls the bubble
up and the contact line force that keeps the bubble down,
causes perturbations in the liquid vapour interface, which in
turn affects the velocity of contact line. This cyclic process
is at the root of the stickslip behaviour of the liquid vapour
interface.
Numerical Model
The nature of the interface between two fluids has been the
subject many investigations. In the case where the interface
is a plane of discontinuity that separates regions of different
phases, tracking it as a function of time becomes a
challenging task, especially if it develops topological
singularities. To overcome such problems one typically
resorts to level set methods. A complete numerical
simulation of a growing and departing bubble during partial
nucleate boiling was performed by Son et al. (1999) by
employing the level set method (LSM). Since then the level
set method has been used by various researchers for
studying boiling related phenomena.
A persistent problem with the level set method is that it
leaks phase from the typically less dense phase to the more
dense phase. Olsson and Kreiss (2005) and Olsson, Kreiss
and Zahedi (2007) introduced a conservative level set
method which was applied to simulate single bubble
dynamics in nucleate pool boiling using a commercial
software package COMSOL, in Sajid and Bennacer (2009).
In this work we extend the model to include the contribution
of the dynamic contact angle between surface tension forces
and buoyancy forces. It is a hypothesis that has not yet been
justified in literature. The prescribed values of the apparent
receding and the apparent advancing contact angles (RCA
and ACA) define the limits of the instantaneous apparent
contact angle. The transition between receding contact angle
towards the advancing contact angle takes place based on
the balance of buoyancy forces and surface tension forces.
In this approach the liquidvapor interface has the receding
contact angle for a buoyancy force less than or equal to the
minimum contact force, i.e. contact force in effect at
receding contact angle. While the interface forms advancing
contact angle with the surface when the buoyancy force is
greater than or equal to the maximum contact force, i.e.
contact force in effect at advancing contact angle, with
intermediate values of OCA for buoyancy force in between
the minimum and maximum values. We model the contact
angle as a hyperbolic tangent function which will be valid
through out the simulation.
SACA+ Rca + tanh nF g0 ACA R
A 2 F,, l ,Y 2)
Where, Fg is a function of the equivalent bubble diameter
and the gravitational constant, g. While Fo,mmi and Fo,ma are
functions of the bubble base diameter and the known values
of the minimum receding contact angle and maximum
advancing contact angle respectively.
Governing Equations
The stabilized level set advection equation, from Sajid and
Bennacer (2009).
kt+V.[+l(0)F (V0)]=
at
The density and viscosity of the liquid and vapor phases is
defined by scalar fields that use the level set function to
distinguish between phases.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
S= p, + (p, p~ 4
=U + (ut,U )0 5
The initial level set phase function is defined as.
0
I= tanh +1d(x)
1
for d < 
for e
for e
Where 0 = 0 in vapor phase and $ = 1 in liquid phase with
intermediate values for the liquid vapor interface. As
described in Sajid and Bennacer (2009) Solving
successively for Eq. 3 accomplishes the transport of the 0 =
0.5 isocontour, while preserving the shape of the profile
and ensuring the conservation of 0 without reinitializing the
level set function at every time step.
The interface normal vector and interface curvature are
obtained from the level set function gradient and the
divergence of the gradient.
r rv
F Vo
KC = V nF = V K1 8
The level set equation is coupled to the Navier Stokes
equation through the level set function and the velocity
field.
P + Vu =V p+ pB+F,
p, (T T,, ) + V (Vi + V a )
Where F, is the surface tension force at the interface
evaluated using the continuum surface force method of
Brackbill et al (1992).
F, = WCF
Where, 6 is the Dirac delta function localizing the surface
tension force to the interface. The continuity equation is
derived by including the effects of volume expansion due to
liquidvapor phase change at the interface.
V. = ep .Vp
p
Where the vapor velocity at the interface due to evaporation
is defined in terms of the vapor mass flux
where VT
where =
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
implementation of contact angle OCA, as follows:
n,, i, = cos(CA)
And the energy equation is;
pCP +i.VT = V .VT +Q 13
The source term Q is used to absorb the latent heat and
retain temperature in the vapor phase at Tsat. The conductive
and convective terms of the energy equation are active only
in the liquid phase.
Characteristic Scales
The governing equations are implemented in
nondimensional form using the following characteristic Cas
scales. cont
valu
Cas
IR= 14 dyn
g(pp,9 con
Cas
S15 smo
tR 15 betv
g betv
(Eq
UR = g 16
CaseI:
The nondimensional temperature is defined as:
The bul
from the
* sat 17 on a hea
T T shown i
0.003
Computational Domain
An axially symmetric domain is considered.
05 o l o.
Figure 2: Computational Domain
1 1.
0.002
0
E
Axial symmetry conditions are applied on the first boundary,
0Q1. For the second boundary (T__, the component of the
velocity normal to the wall, u has to be zero. The noslip
condition is relaxed to navierslip, to facilitate
Symmetry conditions are applied on the right boundary
while the upper boundary has open boundary conditions.
The domain is initially at saturation conditions for pure
water with a spherical bubble centered at the origin with a
diameter of 0.5mm. And the wall superheat, AT = 7 K. The
receding and advancing contact angles are 48 and 61
respectively.
Results
Three approaches to contact angle modelling were studied.
eI: Static contact angle approach, here the apparent
act angle OCA is taken as the equivalent of the mean
e of the advancing and receding contact angles
eII: Contact line velocity approach, OCA is a
amic contact angle that varies as a function of
act line velocity uCL (Eq 1).
eIII: Force balance approach, OCA is modelled as a
oth hyperbolic tangent function of the force balance
veen the contact line force and the buoyancy force
2).
Static contact angle approach
ble diameter and bubble base diameter obtained
Numerical simulation of a single bubble growing
ited surface with a static contact angle of 54.5o are
i figure 3.
O Bubble diameter
A Bubble base diameter
0 A
0
0 &
0 A
A
0000
0000
ooo
oooo
A
,,, ,,,,I,,, I, ,, ,, I,,
0 1 2 3 4 5 6 7 8 9 10
Time [t/tO]
Figure 3: Bubble growth for static contact angle 54.5o
The contact line velocity is plotted in figure 4. Because the
interface is endowed with a thickness in the level set method
exact values of the contact line velocity can not be obtained.
The approximate value of the interface velocity is obtained
by integrating the velocity term over the interface near the
contact point and dividing by the length of the interface
covered.
uevp =
P,
 I
8 92 "824
Ls? W42
ci
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Um,. The limiting value of the interface velocity, Uma is
chosen to be 0.015m/s after considering the evolution of the
interface velocity for a static contact angle of 54.5"
0.1 1
0.05
.2
0
>
A
A a A
0.05
A
A Aa
a
A A~
a A
A
U. 
0 2 4 6 8 10
Time [titO]
Figure 4: Contact line velocity for static contact angle
54.5
The fluctuations in the contact line velocity are caused by
discrete approach which can not accurately capture
continuous phenomena. There fore it would be difficult to
model the dynamic contact angle on the basis of the contact
line velocity.
9E05
A Contaact line force
o Buoyancyforce
6E05 
0
L
3E05
o
0
A A 0
A C
A
A
A
0I
0 1 2 3 4 5 6 7 8 9 10
Time [t/tO]
Figure 5: Force balance for static contact angle 54.5"
Comparing the contact line velocity from Fig 4 with the
force balance in Fig5, we observe that the absolute value of
the interface velocity is proportional to the absolute
difference between the buoyancy force and contact force.
We also note that as the buoyancy force overtakes the
contact force the interface velocity changes signs. Thus,
there is a strong relation between the balance of forces and
the contact line velocity. To be more precise, the contact line
velocity is inversely proportional to Fg F,.
CaseII: Contact line velocity approach
In this case the instantaneous apparent contact angle, OCA is
a function of the interface velocity that changes linearly
from the receding contact angle value, ORCA to the advancing
contact angle value, OACA as the interface velocity goes from
a chosen positive limiting value, um, to a negative one,
A
a
a a
^ AA^ &A A
A
aa & a&Z & &a
a~ naa a
A A
0.1 F
SI, I I ,
0" 2 4 6 8 10
Time [t/tO]
Figure 6: Contact line velocity for dynamic contact angle
48"61" with UcL,max = 0.015m/s
Since the contact angle is a function of the contact line
velocity alone, it is sensitive to changes in the interface
velocity. The interface velocity depends on the rate of
expansion or contraction of bubble. Which in turn depend
on many factors, including the balance of forces i.e. surface
tension, bubble momentum etc. Thus, the interface
behaviour is caught in a cyclic process where the interface
velocity affects the contact angle and the contact angle
affects the interface velocity by altering the surface tension
and bubble momentum. The result of this cyclic process is
visible in the fluctuation of the dynamic contact angle
between ORCA and OACA shown below.
aaaa
"AA A AA A
A AAnAAAAA
' A
AA
aA A a A
AA AA A A
00 1 2 3 4 5 6 7 8 9 10 11
Time [t/tO]
Figure 7: Apparent contact angle for a velocity based
dynamic contact angle approach with UcL,max = 0.015m/s
The force balance near the contact point influences the
ability of the bubble base to expand or contract, inducing
the slip/stick behaviour described in the previous case.
Dussan (1979) explains slip/stick behaviour as when the
motion of a contact line appears to be unsteady and
spasmatic when the contact angle is changing below a
certain smallest reported speed. Beyond this range the
0.15 r
I. I .... I .. ..... I . .. .. . . I . . I
contact line movement ought to be smooth. However we
note that the contact angle is oscillating continuously
leading to a solution that is computationally expensive and
shows pronounced slip/stick behaviour of the contact line.
CaseIII: Force balance approach
Numerical calculations are carried out with the apparent
instantaneous contact angle at the bubble base obtained as a
function of the force balance between buoyancy force and
the contact line force, which is shown in figure 8.
0.00012
8E05
o
IE
4E05
a Contaact line force
O Buoyancyforce
a o
A o
A o
OR?,.
AO
a
oo
oo
C
6A A. "
4 Al
I,,,, I,,,,I,,,, I ,,,, I ,,,, I ,,,, I ,,,,I,
a
,,I,,,, I
0 1 2 3 4 5 6 7 8 9 10 11 12
Time [t/t0]
Figure 8: Force balance for a force balance dynamic
contact angle approach.
Comparing Fg and Fo between the CaseI and CaseII, we
note that the buoyancy force, Fg overtakes the contact force
Fo at a slightly later time for CaseIII. This is because as the
buoyancy force tries to cross the contact line force the
contact angle begins to change from minimum receding
value towards maximum advancing value (Fig 9).
70 
30 
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
delay in the buoyancy force crossing over the contact line
force. As a side effect the there is more evaporation and
hence greater bubble growth.
We not that as the contact angle changes in a smooth fashion
from ORCA to OACA the slip/stick behaviour, exhibited in the
fluctuations of the contact line force at the beginning and
intermediate stages of bubble growth (Fig 8) stops.
Discussion
Ramanujapu and Dhir (1999) experimentally studied the
dynamic contact angle at the base of a vapor bubble
growing on a horizontal surface during nucleate pool boiling
of deionized water. The contact angle was found to vary
between 4861 for AT = 7K, which corresponds to our test
case. Table1 summarizes the results of the three studies in
comparison with the experimental results of Ramanujapu
and Dhir (1999) of the dynamic contact angle and its effects
on the dynamics of a single bubble growing on a heated
surface during nucleate pool boiling.
There is excellent agreement between the experimentally
observed bubble departure diameter and the bubble
departure diameter obtained by numerical simulation using
a force balance model for the dynamic contact angle,
CaseIII.
Diameter Time
[mm] [s]
Case I 2.56 0.151
Case II 2.65 0.167
Case III 2.80 0.178
Experimental 2.82 0.042
Results
Table 1: Bubble characteristic at departure
for various test cases.
However there is considerable difference in bubble
departure time between the experimental and numerical
studies. This is due to the absence of a microlayer
evaporation model as well as a thermal boundary layer in
the initial condition.
The heat absorbed by the system during the transformation
from the liquid to the vapor phase for the three cases of
dynamic contact angle approaches is shown in the following
figure.
In the case of the static mean contact angle the latent heat is
high at the beginning since the static mean contact angle of
2 4 6 8 10 12 54.5o is greater than the receding contact angle 48 and the
Time [t/tO] heat flux available for evaporation is proportional to the
9: Apparent contact angle for a force balance contact angle. Similarly at the end of the simulation the heat
ynamic contact angle approach. transfer for the static contact angle is less as the mean
contact angle is less than the advancing contact angle, 68.
turn tries to pull the bubble base outwards causing However, the velocity base approach doesn't respect this
ase in the bubble base diameter. This causes an criterion. This is because of the excessive stick/slip
in the contact line force which is responsible for the behaviour of the interface which limits the amount of heat
10
a
Figure
based c
This in
an incre
increase
000
o0
flux by causing the contact angle to fluctuate.
I . I I I o I I I''000000
0.01250 0.05 0.1 0.15 0.2
Time [s]
Figure 10: Apparent contact angle for a force balance
based dynamic contact angle approach.
Conclusions
The influence of a dynamic apparent contact angle on the
dynamics of a single vapor bubble nucleating on a heated
wall were included in the numerical simulation using a
conservative level set method that liberates from the
necessity of reinitialization while considering different
approaches of contact angle at the bubble base from which
the following conclusions can be drawn.
1. The static mean contact angle model is the simplest
to implement but under predicts the bubble
departure characteristics.
2. The contact line velocity based approach is
susceptible to discretization errors while estimating
contact line velocity and excessive stick/slip
behavior leading to loss of accuracy in the
determination of bubble departure characteristics.
3. The force balance approach provides an excellent
alternative to modeling dynamic contact angles that
is based on the fundamental physics of contact
angle variation. This approach provides the most
accurate results for bubble departure characteristics
in comparison with experimentally measured data.
Inorder to completely predict bubble dynamics over time the
microlayer evaporation needs to be taken into account
Acknowledgements
The work was conducted in L2MGC at the University of
CergyPontoise with funding from the Higher Education
Commission (HEC) of the government of Pakistan whose
support is gratefully acknowledged.
References
1. Abarajith, H. S., and Dhir, V. K. Effect of Contact
Angle on the Dynamics of a Single Bubble During
Pool Boiling. Proceedings of the ASME IMECE,
^ Casel
o Case II
O Case II
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
New Orleans, LA, (2002), IMECE200233876.
2. Dussan, V.E.B. Spreading of liquids on solid
surfaces static and dynamic contact lines. Ann.
Rev. Fluid Mech. 11, (1979), 371400.
3. Francois, M., and Shyy, W., Computations of drop
dynamics with the immersed boundary method Part
2: drop impact and heat transfer. Numer. Heat
Transfer B 44, (2003), 119143.
4. Huang, H., Liang, D., and Wetton, B., Computation
of a moving drop/bubble on a solid surface using a
fronttracking method. Commun. Math. Sci. 2,
(2004), 535552.
5. Lam, C.N.C., Wu, R., Li, D., Hair, M.L., and
Neumann, A.W., Study of the advancing and
receding contact angles: liquid sorption as a cause
of contact angle hysteresis. Adv Colloid Interface
Sci 96, (2002), 169191.
6. Mukherjee, A., and Kandlikar, S.G., Numerical
study of single bubbles with dynamic contact angle
during nucleate pool boiling. Int. J. Heat Mass
Transfer 50, (2007), 127138.
7. Olsson E., Kreiss G., A conservative level set
method for two phase flow, J. Comput. Phys. 210
(2005) 225246.
8. Olsson E., Kreiss G., Zahedi S, A conservative
level set method for two phase flow II, J. Comput.
Phys. 225 (2007) 785807
9. Ramanujapu, N., and Dhir, V.K., Dynamics of
Contact Angle During Growth and Detachment of
a Vapor Bubble at a Single Nucleation Site. In;
Proceedings of 5th ASME/JSME Joint Thermal
Engineering Conference [CDROM] American
Society of Mechanical Engineers, San Diego, CA,
(1999).
10. Sajid M., and Bennacer R., Simulation of single
bubble dynamics in nucleate pool boiling using a
conservative level set method, Proceedings of the
ASME 2009 Heat Transfer Summer Conference, San
Francisco, California, July 2009
11. Smith, K.A., Ottino, J.M., and Warren, P.B. Simple
representation of contactline dynamics in a
levelset model of an immiscible fluid interface.
Ind. Eng. Chem. Res. 44, (2005), 11941198.
12. Son G., Dhir V.K., Ramanujapu N., Dynamics and
heat transfer associated with a single bubble during
nucleate boiling on a horizontal surface, J. Heat
Transfer 121 (1999) 623631.
