7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Dynamic Contact Angle of a Spreading Drop Represented by Fronttracking Simulation
Yasufumi Yamamoto and Tomomasa Uemura
Department of Mechanical Engineering, Kansai University
3335 Yamatecho, Suita, Osaka 5648680, JAPAN
yamayasu @ kansaiu.ac.jp
Keywords: Dynamic Contact Angle, Fronttracking, StickSlip, Velocity Dependency, Slip Length, Rolling Belt
Abstract
In this study, we present a fronttracking method containing an accurate representation of interfacial tensions on solid surface.
A drop growing and shrinking on a flat plate is simulated employing two moving contact line models, the classical Navier slip
condition and a rolling belt model. Numerical results with different parameters are compared, then the relations between
contact angle dynamics and the parameters are examined. We can find that contact angle hysteresis and stickslip motion are
naturally represented by our fronttracking simulations with the rolling belt model and the velocity dependency of contact
angle can be represented by the slip length model.
Introduction
"Wetting" liquidgas interface contacts solid surface is
very important phenomenon in many industrial processes.
Because the wettability of solid surface affects many
phenomena such as boiling, condensation, coating,
lubrication and so on, understanding of the wetting
mechanism and modeling of the wettability are very
important issues (e.g. Dussan V 1979, de Gennes 1985).
For gasliquid two phase flows, nowadays, numerical
simulations using volume of fluid, levelset, diffused
interface etc. have been frequently performed. In some of
them, the wettability of solid surface is taken into account
(Renardy, Renardy & Li 2001, Himeno et al. 2004, Spelt
2005, Mukherjee & Kandlikar 2007, Takada et al. 2008).
However, in most cases, the contact angle is fixed as a given
static contact angle (Renardy, Renardy & Li 2001, Himeno
et al. 2004, Takada et al. 2008), or is given by an empirical
correlation with the contact line velocity (Spelt 2005,
Mukherjee & Kandlikar 2007). In those implicit interface
methods, it is natural because the interfacial tension effect is
evaluated by the curvature calculated from spatial
derivatives of phase indicator function like CSF model
(Brackbill, Kothe & Zemach 1992). Such models need a
boundary condition of gradient of indicator, namely
interfacial shape need to be prescribed.
On the other hand, in fronttracking method (Unverdi &
Tryggvason 1992, Yamamoto, Yamauchi & Uemura 2008)
representing the interface explicitly by connected marker
points, the interfacial tension effect can be directly
evaluated as pull strength between the marker connected
elements. Such representation can naturally treat the
interfacial tension effect even at three phase contact location
(Yamamoto & Uemura 2008). Thus, as the computational
boundary condition, if not interfacial shape but interfacial
tensions between solidliquid and solidgas are given,
dynamic contact angle may emerge spontaneously by
balance of viscous stress, pressure and the tensions.
In the present study, we propose a contact line movement
model in order to represent real dynamic contact conditions
like contact angle hysteresis and stickslip motion. Then a
drop growing and shrinking on a flat plate is simulated by a
fronttracking method. The new model for contact line
movement based on the original characteristics of
fronttracking method is proposed. This new model and the
classical Navier slip model are employed, then the relations
between contact angle dynamics and those model
parameters are examined.
Nomenclature
Ca capillary number (= Usi,p/o) ()
D weighting function for distribution (m2)
F force due to interfacial tension on fixed grid
(Nm 1)
f interfacial tension force per unit interface length
(Nm2)
G gradient of phase indicator function (m1)
g gravitational acceleration (m s 2)
h height of drop (m)
I phase indicator function ()
k marker number ()
'lap capillary length (m)
n unit normal ()
p pressure (Pa)
Q flow rate from nozzle (m3 s 1)
s coordinate along interface (m)
t unit tangent ()
t time (s)
Us,p slip velocity of contact line (m s1)
u velocity vector (m s 1)
u radial (horizontal) component of velocity (m s 1)
V volume of liquid drop (m3)
v axial (vertical) component of velocity (m s ')
x position vector of fixed grid point (m)
x, position vector of front marker or element (m)
x radial (horizontal) coordinate (m)
y axial (vertical) coordinate (m)
Greek letters
A grid spacing (m)
A, length of front element (m)
0 approximated delta function (n1)
A1 slip length (m)
A2 threshold distance for rolling belt model (m)
p viscosity (Pa s)
0 dynamic contact angle
O, static contact angle
p density (kg m 3)
a surface tension of liquid (N n1)
s,g surface tension of solid (N n1)
oa interfacial tension between liquid and solid
(N n1)
Numerical Method and Modeling
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
viscosity of the multiphase fluid are given by,
p = IP, + (l I)Pg = IPu + (1 I)P g'
where subscript 1 and g represent liquid and gas,
respectively.
Gradient of the indicator function G is related to the unit
normal n(x,) calculated on the marker elements (line
segments) through the delta function (1),
G(x) VI(x)= D(x x,)n(x,)A,
x,
where A, is the element length. Taking divergence of
equation (4), Poisson equation of is obtained as,
V21=V.G.
By solving this equation by a successive over relaxation
(SOR) method, the smoothed Heaviside function I is
obtained on the fixed grid points.
Furthermore, force due to interfacial tension, F, described
later is considered, then basic equations are as follows,
In this study, axisymmetric field is treated, so description
hereafter is on axisymmetric twodimensional coordinates. x
represents a radial (horizontal) component and y represents
an axial (vertical) component.
Fronttracking method
In fronttracking method, liquidgas interfaces are tracked
by marker points, and the interfaces are represented by the
connected elements of the marker points. A multiphase flow
field in the object system is treated as one fluid on a fixed
rectangular grid system with assistance of the marker
elements.
Information interchange between the fixed grid system and
the marker system is done through a weighting function,
which is a smoothed delta function approximated by,
D(x) = q(x)q(y), (1)
S+cos2 x<2A
I 4A x2A 2
0 x> 2A
where A is grid spacing. For example, velocity of marker
point, u(x,), is calculated by velocity at fixed grid, u(x), and
equation (1), as,
u(x,)= D(x x,)u(x)A2. (2)
x
The marker point x, is updated by integration of this
velocity u(x,). After that, marker distribution on the
interface is rearranged to keep the resolution uniformity by
addition and deletion of marker points.
The phase indicator function I is considered on every grid
points. When I is 0 in gas and is 1 in liquid, the interfaces
are represented by 1=0.5 contour lines. The density and the
Vu= 0,
p + V.uu
( a~t
Vp + V (Vu + Vu)+ F + pg (7)
In this study, equations (6) and (7) are solved by a
simplified marker and cell (SMAC) method on a staggered
grid system. The spatial derivatives are approximated by a
secondorder central finite difference. The advection term is
treated by a second order AdamsBashforth scheme and the
viscous term is by the secondorder CrankNicolson
method.
Interfacial tension evaluation and the treatment at the
contact line
The force F due to interfacial tension is given by
distribution of the tensile forces f pulling each line segments
tangentially. The unit tangent at marker point is calculated
by analytical derivative of an interpolation polynomial using
3 points. Using the unit tangent, the force acting the line
element is evaluated by,
fluid B
fluid A
fluid C
Figure 1: Evaluation of interfacial tension force at 3phase
contact point in fronttracking method. k means marker
number.
= ot (ot)+ (ot)
fA s = a = (ot)A (ot) (8)
A cgs A,
where s is a coordinate along the interface, superscripts +
and are the both edges of the element, respectively. The
effect of curvature around the symmetry axis is calculated
by the radial distance and the unit normal and added to
equation (8) separately.
We have been confirmed that the interpolation polynomial
constructed by onesided 3 points is enough accurate when a
side of the line segment is the three phase contact point
(Yamamoto & Uemura 2008). Using the unit tangent and
each interfacial tensions, the force acting the element
adjacent to the three fluid contact point is evaluated by (see
Figure 1),
fAs = (at)+ (ot) = crAt + CrCtB,O + C CAt,o. (9)
In case of flat solid plate like in this study (for example,
fluid C in Figure 1 is assumed as solid phase and the surface
is flat), unit tangents of gassolid interface and liquidsolid
interface are known. Then, the interfacial tension condition
on the solid wall is given by one parameter, 0,, the static
contact angle, using Young's equation,
UBC CA = 'sg (7s = UcOS O .
This relation only gives the interfacial tension condition on
the solid wall not decides the instantaneous dynamic contact
angle 0. Then, horizontal component of equation (9) is
represented by,
(fAs)x = o(cos0s cos9). (11)
Models for contact line movement
In order to evaluate the effect of interfacial tensions by
equation (11), which marker is the contact line must be
known explicitly. However, if a marker is set on the wall
and nonslip boundary condition is strictly applied, the
contact line marker can never move. In this study, a popular
slip model and a new model are used in order to achieve the
movable contact line marker.
The first model is based on the Navier condition (e.g. Huh
& Scriven 1971) considering molecular scale slip by the
velocity gradient and the slip length i,
Us =' (12)
If the contact line particle moves very slowly, a fluid
element may climb over and pass the previous contact line.
For example, Dussan V and Davis (1974) observed the
"rolling" motion of surface of a drop, which slides down an
inclined plate. This rolling motion of liquidgas and
liquidsolid interfaces is considered to be natural and
desired to be modeled. In this study, the interfaces of
liquidgas and liquidsolid are modeled to form a rolling belt
like Caterpillar track, which replaces contact point by one
after another. In this rolling belt model, if a marker near the
contact line approaches to the solid wall within the threshold
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
y*
C axisymmetric
ot hemisphere at t=Os
SRm=3mm=48Ax
i gas(B)
) liquid
(A) x
ozze=0mmsolid(C)
Y IwNozzle=0.5mm
Figure 2: Configuration of the simulated model.
value A2, the marker is replaced as the contact line one. This
unique modeling is just based on the fronttracking method.
In this study, these two model parameters with length scale,
2i and L2 are examined. However, we found that those
parameters should be used in nondimensional form divided
by the grid spacing. The reason is explained as follows. The
interfacial tension force at the contact line is distributed
through equation (1) as,
F(x)= CD(x x,)f(x,)A (13)
so F is proportional to 1/A2 due to the weighting function. It
may be supposed that the viscous force and the interfacial
tension at the contact line are almost balanced, then,
lu U, 1 2 (coso cosO),
y 2 A 2 A
where u,,1/2 is a velocity component nearest to the contact
line discretized on the staggered grid. It must be noticed that
velocity component nearest to the contact line does not
depend on the grid spacing. So the slip velocity is
approximated by,
u i 1 = 1 0 o (cos O cost).
A/2 A c
Thus, the slip velocity depends on the grid spacing A, so
1L1/A rather than .i should be treated as the model parameter.
For the rolling belt model parameter .2, the same treatment
is suitable because of equation (14).
Simulation Condition
As the target of the present simulation, a drop growing and
shrinking on a flat plate by liquid ejection and suction
through a nozzle shown in figure 2 is considered. The
density and the viscosity of the liquid are pi=925 kg/m3 and
Pu =9.25x10 2Pa s, respectively. And those of the
surrounding gas are 1.19kg/m3 and 1.85x10 Pa s. The
interfacial tension a between liquid and gas is
1.69x102N/m. The capillary length l/p (=(Uc/(g))1 2) is
calculated as l/ap =1.36 mm. The interfacial tension
condition for the solid wall is given by the static contact
angle 0 ,=90.
At initial state, all the phase is stationary and the shape of
the drop is set to be hemispherical with radius of 3mm. A
0 2 4 6 8 10
time [s]
0 2 4 6 8 10
time [s]
4 4.2 4.4 4.6 4.8 5
time [s]
Figure 3: Time history of contact line position. (a) by the
slip length model with z2=0, (b) by the rolling belt model
with i =0, (c) partial magnification of (b).
nozzle with 0.5mm radius is set on the symmetry axis. The
liquid is ejected in seconds after 1 sec settling time, and
then sucked in 4 seconds. The flow rate of ejection and
suction is set to 50mm3/s.
As conditions for numerical computation, the number of
grid points inside the initial drop radius is set to 48 and the
number of all grid points is 256x128. So the grid spacing is
about 0.05/cap. The time increment is set to 1x10 s.
Element length for the fronttracking is arranged to be
between 0.60.8 times of the fixed grid spacing. However,
in case using the rolling belt model, the length of the
element having the contact line is not controlled.
Results and Discussion
Contact line movement
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
slip stick slip
Sfor 60ms
(a) 4.212s (b) 4.213s (c) 4.273s (d) 4.274s
stick ms
s lip stick slip
for 42mfor95ms
(e) 4.316s
(f) 4.317s
(g) 4.412s
(h) 4.413s
o or S
stick slip stick slip
for 18ms s 7for 28ms
(i) 4.431s (j) 4.432s (k) 4.460s (1) 4.461s
Figure 4: Instantaneous shape of "advancing" interface for
case of /2/A=0.2, slip motion observation is done by
animation with ims frame period.
" 1
S1
'b
0 2 4 6
time [s]
A/A:
"_1
0
C5
0 2 4 6
time [s]
8 10
8 10
Figure 5: Time history of contact line velocity. (a) by the
slip length model with A2=0, (b) by the rolling belt model
with 1 =0.
Figure 3 shows the time histories of the contact line position,
(a) only with the slip length model and (b) only with the
rolling belt model. In figure 3(a), it is found that the contact
lines move smoothly by the slip length model and
differences by the length parameters cannot be found. On
the other hand, the contact lines by the rolling belt model
move intermittently and the profiles at shrinking depend on
the parameter 2z. By figure 3(c), which is the partial
magnification of figure 3(b) at the range from 4 to 5 seconds,
intermittency of the moving patterns can be found clearly.
What you should pay attention to is not periodical jump
related with element interval is occurred but irregular
repeating of moving and unmoving, socalled "stickslip"
motion is observed. Figure 4 shows the instantaneous
profiles of the marker connected elements for case of
A2/A=0.2. Comparing the marker positions on the wall and
stationary mesh, sticking and slipping can be recognized.
After a slip from (a) to (b), the contact line marker sticks to
the wall for 60 milliseconds long. Then, stickslip motions
with irregular periods are repeated.
Figure 5 shows the time histories of the contact line velocity,
(a) only with the slip length model and (b) only with the
rolling belt model. In order to obtain the locally time
averaged velocity while observing stickslip motion for (b),
velocity is calculated by time derivative of 5th order
polynomials of the profiles shown in figure (3), separately
treated for advancing and receding. We can find that the
contact line movement models and parameters do not affect
the drop expanding velocity. In this figure, velocity of a
radialexpanding cylinder with a constant height is also
shown. In such case, the expanding front position x and the
front velocity u are calculated as follows,
F,, +Qt*
x= d
STh
1 Q
2 h(V, +Qt)
Where, Q is the flow rate from nozzle, h, the constant height,
V*,,, the volume at the starting time of growing or
shrinking, t the time from such starting times. The profiles
with h=2mm shown in figure 5 agrees with all the results
very well. In case of shrinking with the rolling belt models,
velocity profiles are just shifted due to the differences of
shrinking start times but the model parameter does not
change the velocity directly. Contact line velocity may be
decided by the model parameters of contact line movement
in case of wettability driven flow like liquid column in a
capillary tube. On the other hand, these parameters in this
study do not affect the contact line velocity essentially in the
case that contact line is moved by forced flow field. Further
exanimation must be needed and is left for future work.
Dynamic Contact Angle
Figure 6 shows the relation between the dynamic contact
angle and the contact line velocity. For (b), velocity data are
plotted only when the derivatives of polynomial are
available. The contact angles at the onset of advancing and
receding are also plotted. Evaluation of the contact angle is
from unit tangent of the markerconnected element at the
contact line.
At zero velocity by the slip length model, the advancing
contact angles coincide with the receding ones. So, contact
angle hysteresis cannot be represented by the slip length
model. On the other hand, by the rolling belt model, the
contact angle at the onset of advancing differs from the one
at the onset of receding. So, the rolling belt model can
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
120
12 /A=0.025 (a)
0.05
x 0.1
o 0.2
90 ........
60
2 1 0 1 2
contact line velocity [mm/s]
180
,150
120
j 90
S60
30
0
2 1 0 1 2
contact line velocity [mm/s]
Figure 6: Velocity dependency of dynamic contact angle.
(a) by the slip length model with 22=0, (b) by the rolling
belt model with 2i =0.
5
5
o
10
0
1000Ca []
Figure 7: Relation between nondimensional parameters
in case using the slip length model.
represent contact angle hysteresis. It is found that the
magnitude of the hysteresis can be controlled by the
parameter 2z/A.
Although the contact angles by the rolling belt model do not
show the velocity dependency, the ones by the slip length
model show the velocity dependency. This dependency
becomes smaller with larger slip length. This dependency
can be described from equation (15), by using
nondimensional parameters as,
Ca= C (cos cos ),
A
k L/A=0.05 (b)
A 0.1
X 0.2
0 0.4 1
: ................................
an
where C is a nondimensional constant. Figure 7 shows the
relations between Ca and 0 with different Ai/A. Almost all
the cases lie on one line roughly. By all the numerical
methods that convert the effect of interfacial tension to the
body force at a cell as well as the present method, the
velocity component at the first point from the wall near the
contact line is obtained independently from the mesh size.
As the result, the length parameter Ai based on physical
assumption becomes unphysical, just numerical parameter.
Thus, the numerical slip length Ai/A, the nondimensional
slip velocity Ca and the dynamic contact angle 0 are related
with the simple equation (18) as shown in figure 7. By using
this relation, the parameter Ai/A well representing
experimental observations may be found uniquely from the
relation between Ca and 0.
Conclusions
An axisymmetric drop growing and shrinking on a solid
plate was simulated by the fronttracking method. A suitable
treatment of interfacial tension at the three phase contact
position in the fronttracking method was presented. And a
rolling belt model for representing the contact line
movement utilizing the originality of fronttracking was
introduced. Two model parameters of slip length and rolling
belt replacement threshold were examined. By the rolling
belt model, contact angle hysteresis and stickslip motion
were reproduced. For reproducing the contact angle's
dependency on contact line velocity, the slip length model
was suitable. However, by all the numerical methods that
convert interfacial tension to body force at a cell, the slip
length become not physical but just numerical parameter,
and the parameter can be decided uniquely by the relation
between the contact angle and the contact line velocity in
experiments. Combination of these two models and detail
examination of stickslip motion are left for future work.
References
Brackbill, J.U., Kothe, D.B. and Zemach, C., A continuum
method for modeling surface tension. J Comput. Phys., 100,
335354, (1992).
de Gennes, P.G, Wetting: statics and dynamics. Rev Mod.
Phys. 57, 827863, (1985).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Dussan V, E.B. On the spreading of liquids on solid
surfaces: static and dynamic contact lines. Annu. Rev. Fluid
Mech. 11, 371400, (1979).
Dussan V, E.B. and Davis, S.H., On the motion of a
fluidfluid interface along a solid surface. J. Fluid Mech., 65,
7195, (1974).
Himeno, T, Watanabe, T., Nonaka, S., Naruo, Y. and Inatani,
Y., Numerical analysis of sloshing and wave breaking in a
small vessel by CIPLSM. JSME Int. J B, 47, 709715,
(2'"'4).
Huh, C. and Scriven, L.E., Hydrodynamic model of steady
movement of a solid/liquid/fluid contact line. J. Colloid
Interface Sci., 35, 85101, (1971).
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, 127138,
(2007).
Renardy, M., Renardy, Y. and Li J., Numerical simulation of
moving contact line problems using a volumeoffluid
method. J. Comput. Phys., 171, 243263 (2001).
Spelt, P.D.M, A levelset approach for simulations of flows
with multiple moving contact lines with hysteresis. J.
Comput. Phys., 207, 389404, (2005).
Takada, N., Matsumoto, J., Matsumoto, S. and Ichikawa, N.,
Application of a phasefield method to the numerical
analysis of motions of a twophase fluid with high density
ratio on a solid surface. J. Comput. Sci. Technol., 2, 318329,
(2008).
Unverdi,S.O. and Tryggvason,G., A fronttracking method
for viscous, incompressible, multifluid flows. J Comput.
Phys., 100, 2537, (1992).
Yamamoto, Y. and Uemura, T., Numerical experiment of
drop spreading by fronttracking method accurate
representation of interfacial tension at contact line. J. Japan.
Soc. Exp. Mech., 8, Special, 4348 (2008).
Yamamoto, Y., Yamauchi, M. and Uemura, T., Numerical
simulation of a contaminated water drop sinking in an oil by
a fronttracking method. J Comput. Sci. Technol., 2,
246257, (2008).
