7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A VOF method coupled with a dynamic contact angle model for simulation of twophase
flows with partial wetting
Bogdan A. Nichitan, Iztok Zunb and John R. Thomen
a EPFL STI IGM LTCM, ME G1 464, Station 9, CH1015 Lausanne, SWITZERLAND
b LFDT, Faculty of Mechanical Engineering, University of Ljubljana, Askerceva 6, 1000 Ljubljana, SLOVENIA
bogdan.nichita~ epfl.ch, iztok.zun ~fs.unilj.si and john.thome ~epfl.ch
Keywords: volumeoffluid, contact angle, advancing, receding, FLUENT
Abstract
This paper describes the development and the implementation of a 3D dynamic contact angle model into the VOF
method provided by the commercial CFD code FLUENT for simulations of two phase flows with wetted boundaries.
At present, FLUENT only allows user specified fixed values for the contact angle. This is physically quite limiting
and leads to nonrealistic interfacial effects for these applications. The contact angle becomes important in problems
where the surface tension forces are more \ignlilk~.lill than the inertial forces. Problems are also encountered where
dry patches are formed on channel walls, which give rise to a triple point between the wall, the liquid and the gas.
This is often the case for two phase flows inside micro and mini channels. In FLUENT, the contact angle, which the
fluid is assumed to make with the wall, is used to adjust the surface normal vector to the interface in cells near the
wall. This socalled dynamic boundary condition results in the adjustment of the curvature of the interface near the
wall. In our model, the 3D contact angle is computed based on the volume fraction and the contact line velocity, and
is limited by the experimentally available advancing and receding static contact angles. Several 2D and 3D tests are
presented, which demonstrate the accuracy of our model.
Introduction
Two phase flows have been widely studied over the last
years both numerically and experimentally. From the
point of view of numerics, a large amount of work has
been done regarding the so called "one" fluid models,
where a single set of conservation equations is solved
for the whole domain and the interface is then tracked or
captured.
Interface tracking techniques explicitly track the inter
face with marker particles like in MAC methods (Harlow
& Welch 1965), arbitrary LagrangianEulerian (ALE)
methods (Donea 1983; Hirt et al. 1997; Hughes et
al. 1991), or front tracking methods (Tryggvason et al.
2001; Unverdi & Tryggvason 1992). These interface
tracking methods are very accurate and efficient for flex
ible moving boundaries with small deformations. How
ever their main drawback is that they are difficult to use
in cases where the interface breaks up or coalesces with
another interface. Also additional remeshing is needed
when a large deformation of the interface occurs.
In interface capturing methods, an auxiliary function
is needed. These methods are very robust and have
a wide range of applicability. However, they require
higher mesh resolutions. Examples include volume of
fluid (VOF) (Hirt & Nichols 1981; Youngs 1992; Li
1995; Scardovelly & Zaleski 1999) and level set meth
ods (Sussman et al. 1994, 1998, 1999; Sussman &
Fatemi 1999).
In VOF methods, the interface is given implicitly by
a color function, which is defined as the volume fraction
of one of the fluids within each cell. From this func
tion, a reconstruction of the interface is made and the
interface is then propagated implicitly by updating the
color function. VOF methods are conservative and can
deal with topological changes of the interface. However,
VOF methods cannot accurately compute several impor
tant properties, such as curvature and the normal to the
interface. Moreover, a high order of accuracy is hard to
achieve because of the discontinuities in the color func
tion.
For two phase flows in microchannels, the surface ten
sion forces play an important role in determining the dy
namics of bubbles whereas gravitational forces are gen
erally less important, if not negligible. Also it is very
important to consider the interaction between the bound
aries and the fluids by prescribing or computing the cor
rect contact angle between the interface and the bound
ary. FLUENT allows one to input either a static value
for the contact angle, or to use User Defined Functions
(UDF) to compute the dynamic contact angle. When the
contact line starts to move, a fix static value of the con
tact angle would be inappropriate to use. Therefore the
use of a dynamic contact angle model is much more in
Vapor
Contact line o
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
dicated.
In this paper we propose a 3D dynamic contact angle
model based on volume fraction, interface reconstruc
tion, and advancing and receding static contact angles
available experimentally. To define the advancing con
tact angle, one needs to dynamically add volume to the
liquidphase to determine the maximum volume permit
ted without increasing the interfacial area between the
liquid and solid phases. For the definition of the reced
ing contact angle, one needs to remove volume from the
liquidphase to the minimum volume permitted without
decreasing the interfacial area between the liquid and
solid phases. These two angles refer to the static values.
When the contact line starts to move, the contact angles
will deviate from their equilibrium value to provide the
dynamic advancing and receding contact angles.
The rest of this paper is organized as follows: sec
tion 2 presents the dynamic contact angle model with all
the geometrical considerations, section 3 briefly outlines
the FLUENT discretization procedure and code devel
opment, section 4 presents several numerical validation
tests, and conclusions are given in section 5.
Mathematical formulation
If the surface tension is assumed constant between each
pair of the three phases (liquid, gas and solid), see Fig. 1
Young's equation reads as:
Figure 1: Equilibrium contact line and contact angle.
where in Eq. 4 p is the dynamic viscosity of the liquid, v
is the contact line velocity and a is the liquidgas surface
tension.
Once the upper and the lower limits of the contact an
gle are calculated, we evaluated the dynamic contact an
gle based on the geometrical reconstruction of the inter
face. We now give some details about the VOF method
used by FLUENT to capture the interface.
The volume of fluid method. The basic idea of the
VOF method is to consider a color function, defined as
the volume fraction of one of the fluids within each cell,
to capture the interface. This function will be unity if
the cell is filled with the gas phase, zero if the cell is
filled with the liquid phase and between zero and one if
the cell contains the interface. VOF belongs to the so
called "one" fluid family of methods, where a single set
Of COnservation equations is solved for the domain:
as **"' = GS LS
8(u V(pu~u)
dp
81
In Eq. 1, Os is the static contact angle, aLG is the liquid
gas surface tension, a~s is the solidliquid surface ten
sion, and aGs is the gassolid surface tension. The real
dynamic contact angle can take any value within the in
terval On < 0 < eA, where On and eA are the dynamic
receding and the dynamic advancing contact angles. To
compute these values, we used the correlation of Tanner
(Tanner 1979) for the receding angle, see Eq. 2, and the
correlation of Jiang et al. (Jiang et al. 1979) for the
advancing angle, see Eq. 3:
Vp+V(2pD)+Fst~pg (5)
+ V (pu) =
In Eq. 5 and Eq. 6, u is the velocity vector, p is the
density, t is the time, p is the dynamic viscosity, D
is the rate of deformation tensor with the components
Di = (.. ; + ujge), Fst is the body force due to the
surface tension, and g is the gravity vector.
In VOF, the density and the dynamic viscosity are de
scribed by the following formulae:
s ps = 72Ca
recs R
p(x, t)
p1 + (ps pf)F
coseadv s '
= tanh (4.96Cao02 (3
In Eq. 2, On is the dynamic receding contact angle and
Orec~s is the static receding contact angle. In Eq. 3, BA
is the dynamic advancing contact angle and Badv,s is the
static advancing contact angle. In both Eq. 2 and Eq. 3,
Ca is the capillarity number given by:
where the subscripts I and g denote liquid and gas
phases, respectively.
After solving Eq. 5 and Eq. 6, the color function is
advected with the velocity field by:
8F
+ u VF =0 (7)
Ca = L?
I I I I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
VOF/PLIC algorithm is divided into two parts: a recon
struction step and a propagation step. The key part of
the reconstruction step is the determination of the ori
entation of the segment. This is equivalent to the de
termination of the unit normal vector n to the segment.
The normal vector n and the volume fraction F then
uniquely determine the straight line. The second step
of the VOF algorithm is the propagation. Usually, the
propagation step is done using a fractional step or an
operator split method, which updates the volume frac
tion F by advecting the interface along one spatial di
rection at a time. Intermediate F values are calculated
during this process, and the final F field is obtained only
after advection of the interface along all coordinate di
rections. Although they are more complicated, unsplit
methods exist as well.
The surface tension term in Eq 5 is approximated us
ing the Continuum Surface Force model of Brackbill et
al. (Brackbill et al. 1992) given by the following rela
tion:
pkVF
Fs t = (8)
where in Eq. 8 p is the discontinuous density, k is the
curvature of the interface, and pi and p, are the liquid
and the gas densities respectively.
The contact angle model. A first step in determining
the contact angle between the interface and boundary is
similar to the reconstruction step from the advection al
gorithm of the volume fraction F. Given the normal vec
tor to the interface m = [ml, m2, m3], the equation of
the interface, which is a plane in 3D and a line in 2D, is
given by:
mlJx + m2y + m3z = a (9)
where in Eq. 9 a is the plane constant (line in 2D). Scar
dovelli and Zaleski (Scardovelly & Zaleski 1999) pre
sented analytical relations for determining a. Once the
plane constant a is determined, we used the same re
lation which connects a with the volume fraction F to
express the contact angle, and the following notations
like as in (Fang et al. 2008) g m2/ml, h a/~ml,
f ms/ml. For 3D we get:
In general, a VOF algorithm solves the problem of up
dating the volume fraction field F given the fixed grid,
the velocity field u and the field F at the previous step.
In two dimensions, the interface is considered to be a
continuous piecewise smooth line, and the problem of
its reconstruction is that of finding the approximation to
the section of the interface in each cell, by knowing only
the volume fraction F in that cell and in the neighboring
cells. The simplest type of VOF methods are the Simple
Line Interface Calculation (SLIC) (Noh & Woodward
1976) or the SOLAVOF algorithm (Hirt & Nichols
1981). They are first order accurate in the reconstruction
of the interface. Typically, the reconstructed interface is
made up of a sequence of segments aligned with the grid.
Fig. 2 shows the exact color VOF function for a smooth
circular arc, whilst Fig. (3(a)) shows the interface recon
struction using SLIC method. More accurate VOF tech
niques attempt to fit the interface with piecewise linear
segments. These are known as Piecewise Linear Inter
face Calculation (PLIC) (Li 1995); Fig. 3(b) shows an
interface reconstructed with VOF/PLIC.
0.3
0.01
Figure 2: The exact VOF color function for a circular
are over a square grid.
6gF~ f'3~ (i= my ~( i
mi>
01
(a) SLIC
e.1 oo
0
(b) PLIC
where in Eq. 10 H is the discontinuous Heaviside func
tion. For 2D Eq. 10 becomes:
1 2 2 mii mii2
2 A Amia & i miai1 (11)
Figure 3: Interface reconstruction schemes,
Because F in Eq. 7 is not a continuous function, the
0.01  0
0.
Now observing that tan(0) = Jm m/ms, Eq. 10
becomes:
tan 6 F H Am A
6g
m4Axt *
~i=~~1 92
J1+9"a,,miass
tan0 mi
m am a
mi (12)
a m , 2 tan 6 AM
For the 2D case we have tan(0) = ml/m, so Eq. 11
becomes:
tan(0) 2 2 mi i mi
Eq. 12 and Eq. 13 are solved iteratively for a with a New
ton iterative method. The derivative of Eq. 12 with re
spect to a is given by:
A+B+C
In Eq. 14 A, B and C are given by the following rela
tions:
A = AHAm
6g COS92 9~~[l~ ~~ i mlaz
( 
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
for the cell containing the interface, compute the
plane (line) constant a;
by knowing the position at the current time step and
at the previous time step, a contact line velocity can
be rapidly computed;
with the contact line velocity known, the dynamic
advancing and dynamic receding contact angles can
be computed using Eq. 3 and Eq. 2 using the input
values of static advancing and static receding con
tact angles;
the dynamic contact angle is then computed using
Eq. 12 (or Eq. 13 for the 2D case) and if the value
falls between the receding and advancing angles,
then this value will be returned to the FLUENT
solver, otherwise the advancing or the receding an
gle is returned to the FLUENT solver;
move to the next time step.
FLUENT discretization and code development
FLUENT is a finite volume code, so one needs to solve
the integral form of Eq. 7 which in conservative form
reads as:
H A An gAZ2 tan8 ug
i=, ( a~ a, J+2
mI )
~ ~T~ ~O (puF') =
,& 2
tan0
C=2g sin2 9
tan 8 3
For the unsteady term, FLUENT employs first order Eu
ler discretization. For the convective term, the Green
Gauss theorem is applied and the volume integral is
transformed into a surface integral.
FLUENT uses a midpoint rule integration of the sur
face integral which is secondorder accurate. It also
prOvides five schemes to interpolate the face values
Ff, namely: first order upwind, second order upwind,
power law, QUICK, and MUSCL. Because the MUSCL
(Monotone Upwind Schemes for Conservative Laws)
scheme is less diffusive than the others, it was used in
all our numerical applications to interpolate face val
ues in Eq. 7 and Eq. 5. In Eq. 5, if the pressure field
is known, one can solve for the velocity field. How
eVer, the pressure field is not known a priori and must
be obtained as part of the solution. FLUENT offers sev
eral pressurevelocity coupling algorithms like SIMPLE
(SemiImplicit Pressure Linked Equations), SIMPLEC
(SIMPLE Consistent), and PISO (Pressure Implicit with
Split of Operators). Since FLUENT uses a collocated
grid where pressure and velocity are stored at cell
centers, an interpolation procedure is needed to interpo
late the pressure face values from the cellcenter values.
For twophase flow, FLUENT provides two schemes,
For the 2D case A, B and C in Eq. 14 have the following
form:
miAzi & miAzi 2
m1 01
tans a 1
B= H( ) (6
sin2 8 tan 8
C = 0.0 (16
To summarize, the following algorithm was followed:
'" 1 )`
H (A f) A 
tan 8 / J
2
H (A g f) A g tan 6
H (A 1 f) A 1 
tan 8
A = 1 & A
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and Scardovelli (Manservisi & Scardovelli 2008) and
experimentally by Lim et al. (Lim et al. 2008) and
Grundke et al. (Grundke et al. 2008). Grundke et
al. (Grundke et al. 2008) measured a contact angle
between a sessile water drop and a silicone surface be
tween 70" and 12(7 The error Ef between succesive
mesh refinements is computed using the Li norm, which
is well suited for the interfacial flow problems:
El  H(,) H(Fe da (19)
where H(F) is the discontinuous Heaviside function,
and F, and Fe are solutions from, respectively, re
fined and coarse grids. For the numerical integration of
Eq. 19, a domain of 1000x1000 mesh points was used
and the solution was interpolated from the actual do
main to the computational one. The results obtained are
shown in Table 2. As can be seen, increasing the resolu
tion of the grid leads to a decrease in computation error.
(a) t)s
(c) t=().(4s
namely bodyforceweighted which computes the pres
sure by assuming that the normal gradient of the differ
ence between pressure and body forces is constant, and
PRESTO! (PREssure STaggering Option) which uses
the discrete continuity balance for a "staggered" control
volume about the face to compute the face pressure. In
all our simulations the PISO scheme was used for the
velocitypressure coupling, and the PRESTO! scheme
for the pressure interpolation.
In FLUENT the dynamic contact angle model was im
plemented using User Defined Functions. The values of
the contact angle and other parameters are stored for the
current and the previous time step, for handling the sit
uation when the contact line is shifting from one cell
to another. Once the contact angle is computed and re
turned to the FLUENT solver, the surface normal at the
live cell next to the wall is:
n = n, cos 9 + nt sin 9 (18)
where in Eq. 18 n, and n, are the unit vectors normal
and tangential to the wall, respectively.
Numerical validations
In this section several 2D and 3D simulations are pre
sented for airwater and a coated silicon wafer surface.
2D examples. In the present study, a semicircular
droplet of 2cm diameter is placed on a wall of a rectan
gular domain of 8x2.5cm at (4.0,0.5) cm with a gravity
vector (0.0, 9.81) m/s The fluids used are air and
water, see Table 1 for their physical properties. Fang et
Table 1: Physical properties of the fluids.
Property Density Viscosity Surface Tension
ky/n?" Ns/m2 N/m
liquid 1000 0.001 0.072
gas 1.225 0.0000178
al. (Fang et al. 2008) have measured the contact angle
for airwater on a coated silicon wafer surface and found
that the advancing contact angle is around 105" and that
the receding contact angle is around 65 Consequently,
in our simulations we used the above mentioned values
as limiters for the dynamic contact angle. A no slip wall
boundary condition was used for all boundaries. Be
sides this, for the bottom boundary the dynamic con
tact angle was computed and prescribed. Three differ
ent grids were used: 80x25, 160x50 and 320x100. Fig 4
shows the bubble contour at different times using a grid
of 80x25. Fig. 5 and Fig. 6 show the bubble contour on
a grid of 160x50 and on a grid of 320x100, respectively.
Similar shapes were found numerically by Manservisi
(b) t=().(2s
(d) t=().(6s
Figure 4: Snapshots of the 2D water drop for a grid of
80x25 and gravity vector (0.0,9.81).
Table 2: Error between successive mesh refinements for
2D airwater simulations with gravity vector (0,9.81)
and (9.81,0).
Grid Case 1, Ei, Case 2, Ei,
(0, 9.81) (9.81, 0)
80x25 N.A. N.A.
160x50 1.22e2 7.049e3
320x100 3.359e3 2.829e3
The second simulated case was similar to the above,
however the gravity vector is (9.81,0)nt/s2. Thus we
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) t)s
(c) t=().(4s
(e) t().(8s
Figure 8: Snapshots of the 2D bubble contour for a grid
of 160x50 and gravity vector (9.81,0).
3D examples. In this section we extend the 2D ex
amples shown above to the 3D case. For the first
case a semispherical droplet of 2cm diameter is placed
on a wall of a parallelipipedic domain of 8x2.5x8
cm at (4.0,0.5,4.0) cm. Based on the convergence
study carried out for the 2D case, the grid chosen was
160x50x160. A no slip wall boundary condition was
used for all boundaries. The dynamic contact angle, lim
ited by 65" and 105 was computed and prescribed at
the bottom boundary. Fig. 10 shows the bubble contour
in the XY plane for Z=0.04m at different times. Again
good agreement was found when compared with numer
(e) t=().(8s
Figure 6: Snapshots of the 2D water drop for a grid of
320x100 and gravity vector (0.0,9.81).
expect now to see different contact angles at the two
droplet extremities, using the same boundary conditions
as in the previous case and three different grids: 80x25,
160x50 and 320x100. Once again the error decreased
with increasing mesh densities. Fig 7 shows the droplet
contour at different times using a grid of 80x25. As can
be seen, the drop moves from right to left under the ac
tion of gravity. Fig. 8 and Fig. 9 show the bubble con
tour on grids of 160x50 and on a grid of 320x100, re
spectively. The error Ef between succesive mesh refine
ments is computed as is explained above, and the results
are shown in Table 2.
(a) t)s
(c) t().(4s
(b) t=().(2s
(d) t=().(6s
(a) t)s
(c) t=().(4s
(b) t=().(2s
(d) t=().(6s
(e) t=().(8s
Figure 5: Snapshots of the 2D water drop for a grid of
160x50 and gravity vector (0.0,9.81).
(e) t().(8s
Figure 7: Snapshots of the 2D bubble contour for a grid
of 80x25 and gravity vector (9.81,0).
(a) t)s
(c) t().(4s
(b) t=().(2s
(d) t=().(6s
(b) t=().(2s
(d) t=().(6s
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Fig. 11 shows how the shape and position of the droplet
changed with time.
Figure 11: Snapshots of the 3D bubble contour for a
grid of 160x50x160 and gravity vector (9.81,0.0,0.0) in
plane z=0.04m.
Conclusions
The commercial CFD code FLUENT allows one to use
a static contact angle between the interface and the wall.
In dynamic conditions (i.e. when the contact line starts
to move), this choice of a static value is not physically
appropriate. Therefore, a dynamic contact angle model
was developed and implemented via UDF into the com
mercial CFD code FLUENT for two phase flows with
wetted boundaries.
As a first step, the interface is reconstructed based on
the volume fractions. After the interface position is de
termined, the contact angle is determined based on the
volume fractions and the normal to the interface. The
model has as its lower limit for the contact angle, the re
ceding contact angle, and as its upper himit, the advanc
ing contact angle. The advancing and receding contact
angles are available experimentally for several pairs of
fluids/surfaces, or can be estimated using several pub
lished correlations.
Several 2D and 3D tests were performed for air/water
on a coated silicon wafer surface with different gravity
vectors, which proved the accuracy of our model when
compared to both numerically and experimentally data
available in the literature.
References
J. Brackbill, D. Kothe, C. Zemach, A Continuum
Method for Modeling Surface Tension, Journal of Com
putational Physics 100 (1992) 335354.
Figure 9: Snapshots of the 2D bubble contour for a grid
of 320x100 and gravity vector (9.81,0).
ical results of Manservisi and Scardovelli (Manservisi &
Scardovelli 2008), and when compared with experimen
tal data of Lim et al. (Lim et al. 2008) and Grundke et
al. (Grundke et al. 2008).
(b) t=0.02s
(d) t=0.06s
(a) t0s
(c) t0.04s
(b) t=0.02s
(d) t=0.06s
(a) t0s
(c) t=0.04s
(b) t=0.02s
(d) t=0.06s
(e) t=0.08s
(e) t0.08s
(a) t0s
(c) t0.04s
(e) t=0.08s
Figure 10l: Snapshots of the 3D bubble contour for a
grid of 160x50x160 and gravity vector (0,9.81,0.0) in
plane z=0.04m.
As for the 2D case, the gravity vector was then
changed from (0.0,9.81,0.0) to (9.81,0.0,0.0). A semi
spherical droplet of 2cm diameter was placed on a
wall of a parallelipipedic domain of 8x2.5x8 cm at
(4.0,0.0,4.0) cm using a grid of 160x50x160. A no slip
wall boundary condition was used for all boundaries.
The dynamic contact angle, limited by 650 and 1050,
was prescribed at the bottom boundary. Again different
contact angles at the droplet extremities were expected.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
R. Scardovelli, S. Zaleski, Direct numerical simulation
of free surface and interfacial flow, Annual Review of
Fluid Mechanics 31 (1999) 657603.
R. Scardovelli, S. Zaleski, Analytical Relations Con
necting Linear Interfaces and Volume Fractions in Rect
angular Grids, Joumnal of Computational Physics 164
(2000) 228237.
M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H.
Howell, M. L. Welcome, An Adaptive Level Set Ap
proach for Incompressible Two Phase Flows, Joumnal of
Computational Physics 148 (1999) 81124.
M. Sussman, E. Fatemi, An efficient interface
preserving level set redistancing algorithm and its ap
plication to interfacial incompressible fluid flow, SIAM
Journal on Scientific Computing 20 (4) (1999) 1165
1191.
M. Sussman, E. Fatemi, P. Smereka, S. Osher, An im
proved level set method for incompressible twophase
flows, Computers & Fluids 27 (56) (1998) 663680.
M. Sussman, P. Smereka, S. Osher, A level set approach
for computing solutions to incompressible twophase
flow, Joumnal of Computational Physics 114 (1994) 146
159.
L. H. Tanner, The spreading of silicone oil drops on hor
izontal surfaces, Joumnal of Physics D: Applied Physics
12 (1979) 14731484.
G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al
Rawahi, W. Tauber, J. Han, S. Nas, Y.J. Jan, A Front
Tracking Method for the Computations of Multiphase
Flow, Joumnal of Computational Physics 169 (2001)
708759.
S. O. Unverdi, G. Tryggvason, A Front Tracking Method
for Viscous, Incompressible, Multi fluid Flows, Jour
nal of Computational Physics 100 (1992) 2537.
D. Youngs, Time dependent multimaterial flow with
large fluid distortion, in: K. Morton, M. Baines (eds.),
Numerical Methods for Fluid Dynamics, Academic
Press, New York, 1982, pp. 273285.
J. Donea, Arbitrary LagrangianEulerian finite element
methods, Computational Methods for Transient Anal
isys 1 (1983) 473516.
C. Fang, C. Hidrovo, F. Wang, J. Eaton, K. Goodson,
3D numerical simulation of contact angle hysteresis for
microscale two phase flow, Intemnational Joumnal of Mul
tiphase Flow 34 (2008) 690705.
K. Grundke, S. Michel, G. Knispel, A. Grundler, Wet
tability of silicone and polyether impression materials:
Characterization by surface tension and contact angles
measurements, Colloids and Surfaces A: Physicochemi
cal and Engineering Aspects 317 (2008) 598609.
D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski,
VolumeofFluid Interface Tracking with Smoothed Sur
face Stress Methods for ThreeDimensional Flows, Jour
nal of Computational Physics 152 (1999) 423456.
F. Harlow, J. Welch, Numerical Calculation of Time
Dependent Viscous Incompressible Flow of Fluid with
Free Surface, Physics of Fluids 8 (1965) 21822189.
C. Hirt, A. Amsden, J. Cook, An arbitrary Lagrangian
Eulerian computing method for all flow speeds, Joumnal
of Computational Physcs 135 (1997) 203216.
C. Hirt, B. Nichols, Volume of fluid (VOF) method for
the dynamics of free boundaries, Joumnal of Computa
tional Physics 39 (1981) 201225.
T. Hughes, W. Liu, T. Zimmermann, Lagrangian Eule
rian finite element formulation for incompressible vis
cous flow, Computer Methods in Applied Mechanics
and Engineering 29 (1981) 239349.
T.S. Jiang, S.G. Oh, J. C. Slattery, Correlation for Dy
namic Contact Angle, Joumnal of Colloid and Interface
Science 69 (1979) 7477.
J. Li, Piecewise Linear Interface Calculation, Comptes
Rendus de l'Academie des Sciences Serie II. Fascicule
B Mecanique 320 (1995) 391396.
T. Lim, S. Han, J. Chung, J. T. Chung, S. Ko, Exper
imental study on spreading and evaporation of inkjet
printed picoliter droplet on a heated substrate, Interna
tional Joumnal of Heat and Mass Transfer.
S. Manservisi, R. Scardovelli, A variational approach to
the contact angle dynamics of spreading droplets, Com
puters & Fluids.
W. F. Noh, P. Woodward, SLIC (Simple Line Interface
Calculation), in: Proceeding of the Fifth Intemnational
Conference on Numerical Methods on Fluid Dynamics,
vol. 59 of Lecture Notes in Physics, Springer, 1976.
