1194
Ind. Eng. Chem. Res. 2005, 44, 1194-1198
Simple Representation of Contact-Line Dynamics in a Level-Set Model of an Immiscible Fluid Interface Kurt A. Smith* and Julio M. Ottino Department of Chemical Engineering, Northwestern University, Evanston, Illinois 60208
Patrick B. Warren Unilever R&D Port Sunlight, Quarry Road East, Bebington CH63 3JW, U.K.
We describe a macroscopic level-set method for handling contact-line motion at a solid interface and extract the relationship between an apparent contact angle and the contact-line slip velocity in a simplified contact-line scenario, confirming a theoretical argument for the force balance in the vicinity of the contact line. This approach offers a way to incorporate contact-line dynamics within level-set methods. The motion of a contact line leads to a force singularity in the Navier-Stokes equations if a no-slip boundary condition at the solid surface is assumed.1 This fact has led to the proposal of various microscopic mechanisms (such as molecular diffusion or evaporation-condensation) by which contact-line motion could occur.2 Recently, several numerical techniques have been used to simulate contact-line dynamics. These range from molecular dynamics3 to mesoscopic models4-7 to macroscopic sharp interface methods.8,9 Each of these methods carries advantages and disadvantages. Molecular-scale simulations capture the most physical detail, but their computational intensity prohibits simulations of large length and time scales. Mesoscopic models offer a means to reduce the level of detail while still giving rise to contact-line motion implicitly. Macroscopic techniques face the most difficulty because they contain no mechanism for producing contact-line motion. Instead, a slip velocity must be explicitly specified. However, these methods are widely used to study practical systems containing immiscible fluids. They have the advantage of providing an unambiguous description of the velocity field and interface location. Thus, in many applications it would be valuable to be able to treat contact-line motion accurately within a macroscopic model. In this paper, we develop a treatment of contact-line motion for a finite-difference Navier-Stokes solver paired with a level-set method to track the evolution of the interface (the Navier-Stokes level-set method or NSLSM). The level-set method as proposed by Osher and Sethian10 is a novel method for capturing the evolution of complex interfaces. In this method, the interface is represented as the zero-level surface of a level-set function of one higher dimension. The evolution of the interface is then embedded in an evolution equation for the level-set function itself, which is solved on a fixed Eulerian mesh. The method has many wellestablished advantages over competing front-tracking algorithms. These include a natural handling of topological changes in the interface, easy extension to any number of spatial dimensions, and a fundamental connection to numerical methods for hyperbolic conser* To whom correspondence should be addressed. Present address: LPMCN, Universite Claude Bernard Lyon 1, 69622 Villeurbanne, France.
vation laws, which enables the stable and accurate evolution of corners and cusps in the interface. In the case at hand, the local curvature and slope of the interface at the contact line plays an important role in driving slip. The level-set method is particularly suited to this point. The reason is that the interface position is determined from a smoothly varying levelset function φ. In particular, at a given point x, it is customary to set |φ(x)| to the shortest distance from x to the interface. Thus, the shape of the interface near the contact line can be determined with good precision near the contact line from the values of φ at neighboring grid points. In contrast, in the volume-of-fluid method, φ(x) is defined as the volume fraction of a given fluid in each cell. In any cell that is not intersected by the interface, φ can only have values of 0 (in one fluid) and 1 (in the other fluid). Thus, information about the interface is only contained on a limited number of grid points. This reduces the precision with which one can measure the local interface position. Level sets have been shown to handle a range of free surface flows accurately (see refs 11 and 12), including the case of a fluid/fluid/fluid triple line.13-15 The aim of the present study is to further extend the applicability of level sets to problems involving a contact line. A macroscopic numerical technique, such as NSLSM, is limited in what it can tell us about the nature of slip at the contact line. However, there are open questions regarding the relationship of contact-line dynamics to that of the bulk fluid. We take the approach that is possible to incorporate a description of contact-line motion based on, say, molecular considerations into NSLSM as a boundary condition for motion of the interface. In the present work, we consider a simple relationship between the slip velocity and contact angle, but the method is easily extended to other definitions of slip. Our aim is to show that NSLSM can be used to study the relationship between contact-line and bulk dynamics and also to handle more general interfacial flows with a proper treatment of the contact line. The level-set representation of an interface is as follows. Consider a curve Γ in the plane propagating normal to itself with speed F. We can embed the initial position of the front as the zero-level set of a function φ defined on the plane, the level-set function, and then identify the evolution of φ with the propagation of the
10.1021/ie0498605 CCC: $30.25 © 2005 American Chemical Society Published on Web 07/16/2004
Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1195
front itself through a time-dependent initial value problem. The zero-level set of the evolving function φ defines the interface so that if x is a position on the interface at time t, then
φ(x,t) ) 0
(1)
The interface normal n, its curvature κ, and an interfacial Dirac delta function δ are given by
∇φ |∇φ| κ ) ∇‚n
n)
δ )
{
(2) (3)
[1 + cos(πx/)]/2 0
if |x| < otherwise
(4)
The equations governing the motion of an unsteady, viscous, incompressible multiphase fluid system are the Navier-Stokes equations. Following ref 16, interfacial tension can be expressed as a body force S localized at the interface. Thus, the governing equations can be written in dimensionless form as
∂u ) -∇p - u‚∇u + r∇‚(2D) + S ∂t
(5)
∇‚u ) 0
(6) Figure 1. Shear-driven interface. A contact line is circled.
where u is the velocity, p is the pressure, µ is the viscosity, F is the density, and r ) 1/Re ) µ/Ful. D is the rate of deformation tension, defined as Dij ) (∂jui + ∂iuj)/2. In addition, S ) κδn/We, where the Weber number, We ) Flu2/σ, is the ratio of inertia to interfacial tension, and σ is the interfacial tension coefficient. The advection of the interface is given by
∂φi ) -u‚∇φi ∂t
(7)
The evolution of the interface and velocity field is determined by eqs 5-7. The numerical implementation of the governing equations is as follows. A regular finitedifference grid is used with φ and p defined at the cell centers, labeled (i, j). The horizontal velocity u is defined at the cell faces (i ( 1/2, j), and the vertical velocity v is defined at the cell faces (i, j ( 1/2). Derivatives for the pressure, surface tension, and viscous terms on the right-hand side of eq 5 are calculated by standard central differencing. For the advection terms in eqs 5 and 7, an essentially nonoscillatory scheme, described in ref 17, is used to calculate derivatives. Explicit time discretization is used except, for the viscous term, r∇‚ (2D) in eq 5, which is calculated by a semi-implicit scheme developed in ref 18. This removes the viscous limit on ∆t, which becomes prohibitive at low Re. To satisfy the incompressibility condition (eq 6), a projection method19 is used to determine the pressure field. The first step is to determine an intermediate velocity, which is updated by all terms except the pressure gradient, i.e.
u* ) un + ∆t[-u‚∇u + r∇‚(2D) + S]
(8)
Then u* is not, in general, divergence-free. In the second step, u* is corrected by the pressure term
(∇pF )
un+1 ) u* - ∆t
(9)
Given the requirement that un+1 is divergence-free, the pressure field can be obtained from the Poisson equation
∇‚
(∇pF ) ) - ∇‚u* ∆t
(10)
Further details of the level-set method are given in refs 11 and 12. We consider two immiscible fluids occupying a twodimensional region between parallel plates located at y ) 0 and L such that the interface separating the two fluids forms a contact line at either plate (Figure 1). Boundaries in the x direction are periodic. The plates are moved in opposite directions at constant velocities (Uwall. Assuming that inertial effects are small, there are three forces acting on the system: viscous shear stresses, surface tension, and a “slip force” located at the contact line. Contact-line motion is a nonequilibrium process. It is driven by a force arising from the fact that the system is out of equilibrium, in other words, that the contact angle differs from its equilibrium value. We thus consider the slip velocity as a function of the contact angle, Uslip ) f(θeq - θ), and for small values, we consider only the first term in the Taylor expansion:20
Uslip ) k(θeq - θ)/θeq
(11)
where θeq is the equilibrium contact angle given by the Young-Laplace condition and k is a coefficient with units of velocity. We implement this slip velocity numerically in the following manner. The position of the interface is updated by eq 7. Away from the walls, φ is advected with the continuum fluid velocity u. The x component of velocity of the interface at the wall is Uwall + Uslip, while the fluid velocity is u ) Uwall at the wall. Thus, in solving eq 7, we prescribe a velocity of u ) (Uwall + Uslip, 0) at points lying on the wall. The contact angle θ is
1196
Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 Table 1. Grid Refinement Studya 1/∆x
θa (deg)
displacement
128 256 512
74.3 71.2 69.1
0.129 0.164 0.181
a Comparison of the apparent contact angle θ and the displacea ment of the contact line from its initial position due to shear for β ) 400 and Ca ) 0.05.
Figure 2. Streamlines (a) and velocity field (b) for steady-state flow with p ) 1/2, Ca ) 0.05, and β ) 400. The grid size is 128 × 128.
determined by measuring the slope from the location of the interface at the wall and one cell length away from the wall. The balance of forces in the system is characterized by a capillary number Ca ) µUwall/σ and a dimensionless slip velocity β ) k/Uwall. In order for a steady state to be reached, both Ca and β-1 must be sufficiently small. Thus, steady states are only possible for a restricted area in the (Ca, β-1) plane. Both slip and surface tension forces play a role in opposing shear. Finally, the aspect ratio is p ) X/L, where X is the distance between adjacent interfaces. We consider a symmetrical system where θeq ) π/2 and the viscosities in both phases are equal. A typical steady-state flow, for p ) 1/2, is shown in Figure 2. The box length is 1, and the resolution is ∆x ) 1/128. We note that the flow resembles a driven cavity flow. At steady state, the fluid interfaces play a role similar to that of the side walls in a lid-driven cavity,21
with the main difference being that the tangential component of the velocity can be nonzero on the interface. There are two main circulating cells in a fluid phase, one at the top and one at the bottom. In Figure 2, the midpoint M of an interface is a hyperbolic point, whose stable manifolds originate at the contact lines and whose unstable manifolds are directed toward the center of each phase. This produces the secondary cells seen in Figure 2b, which do not exist in a driven cavity. For p ) 1 (not shown), the two cells merge, except near the center of the system, so that over most of the fluid there is a single circulating flow. This is similar to the change in driven cavity flow when the aspect ratio is varied. To check the stability and accuracy of the method, we performed a grid refinement study. The apparent contact angle, defined by the slope of the straight portion of the interface, and the displacement of the contact line from its initial position are listed in Table 1. The evolution of the flow to steady state is shown in Figure 3. As the interface moves, the flow field develops a “bump”, which grows along the interface. Flow normal to the interface is largest at the tip of the bump. In the present system, the evolution of this bump is limited by the box size. Eventually, the bumps of the upper and lower contact lines intersect at the center of the box and the interface attains a steady shape. In the final (steady state) image of Figure 3, the righthand fluid undergoes a “caterpillar” motion, described by Dussan,22 in which fluid elements lying on the wall approach the contact line, where they are then lifted off the wall and move along the fluid/fluid interface. A necessary condition for this flow is the existence of an injected streamline in the left-hand fluid, which flows toward the contact line, as is clearly seen in the figure. The resulting flow is similar to that in Figure 8 of ref 22. Consider now the balance of forces in the contact region. Our objective is to derive an approximate relation between an apparent contact angle and the parameters Ca and β. For the contact line (at the upper wall) to be stationary in the system frame of reference, it is necessary that Uslip + Uwall ) 0. Thus, from eq 11 the dynamic contact angle depends on β as
θd/θeq ) 1 - β-1
(12)
We first consider the limit β f ∞ so that θd/θeq ) 1. In the outer region, away from the contact line, the interface is nearly straight. The slope of the interface in this region defines an apparent contact angle θa. The value of θa depends on surface tension forces near the contact line, where the interface has large curvature. We approximate the interface as two line segments so that it has an angle of θa outside of a region of distance h from the wall as shown in Figure 4. The curvature is thus concentrated at the kink in the interface, and the contribution of the surface tension forces can be treated as a point force located at a height h. From geometrical
Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005 1197
Figure 3. Evolution to the steady-state flow in Figure 2 for a contact line subjected to sudden shear.
Figure 4. Approximation of curvature in the contact region as a point force. θd ) π/2 at β f ∞.
arguments, the components of the force vector parallel and perpendicular to the wall are
F|surf ) σ cos θaF⊥surf ) σ(sinθa - 1)
(13)
We suppose that the normal component of the force in eq 13 is balanced by nonhydrodynamic effects. The scenario at the contact line can thus be reduced to a calculation of the velocity field as a result of a point force in a fluid at a height h above a moving-plane wall. For the point to be stationary, the horizontal component of the fluid velocity at the point must vanish. Therefore, in the system frame of reference, the force F|surf must be balanced by a hydrodynamic drag force and thus σ cos θa ∝ µUwall or
cos θa ∝ Ca
Figure 5. [θeq - θ]/[θeq - θa] vs y for β ) 20 and Ca ) 0.01 (solid line), 0.015 (dotted line), and 0.02 (dashed line).
(14)
For the more general case of finite β, the horizontal component is
F|surf ) σ(cos θa - cos θd)
(15) Figure 6. [θeq - θ]/[θeq - θa] vs y for β ) 400. Results for Ca ) 0.2, 0.3, 0.4, and 0.5 are shown. Data collapse onto a single curve.
and the appropriate expression for θa is
cos θa ) cos([1 - β-1]θeq) + KCa
(16)
where K is a constant. Because cos θa < 1, this analysis furnishes an explicit expression for the region in the (Ca, β-1) plane where steady-state solutions exist:
cos([1 - β-1]θeq) + KCa < 1
(17)
There is a critical value, βcr ) 1, such that for β-1 > βcr-1 there is no steady state for any value of Ca. This also follows trivially from eq 12 because it is the necessary condition that θd > 0. To verify these predictions numerically, a series of interface shapes were computed for several values of β and Ca. Interfacial profiles are shown in Figures 5 and 6. In all cases, the interface is nearly flat away from
the walls, with θ ) θa. For small β, the shape of the curve depends on Ca. This is because θd ) θ(y)0) is determined by β, while θa ) maxy θ(y) decreases with Ca. Thus, the curves in Figure 5 diverge at the walls. In the large β limit, θd ≈ θeq and thus the normalized data in Figure 6 lie on a single curve that approaches zero at the walls. The dependence of θa on Ca is shown in Figure 7. At the largest values of β, the data are well described by eq 14, while for lower β, there is an offset from the origin, as predicted by eq 16. For all values of β, the data are essentially linear with a slope of K ≈ 5.3. The physical significance of length h should be addressed. It can be roughly thought of as the length scale over which wall interactions are significant. Numerically, it scales with the grid resolution because the value
1198
Ind. Eng. Chem. Res., Vol. 44, No. 5, 2005
Figure 7. cos(θa) vs Ca for several values of β.
of θ used in eq 11 is determined from information within one cell length from the wall. As the resolution is increased, h should decrease. Ultimately, molecularscale phenomena near the contact line are beyond the reach of a macroscopic model. However, our results show that this approach provides a means to incorporate these phenomena into a macroscopic description of the fluid. These results allow for an approach in which the motion of the contact lines within NSLSM is related to molecular-scale simulations or theories by treating the relationship between the apparent contact angle and the wall velocity as an extra constitutive relationship. The desired relation should first be established by mesoscale or microscopic simulation. Then the free parameters in NSLSM are chosen so that eq 16 matches the desired relation. If necessary, we can elaborate the algorithm or generalize eq 11 to improve the match. This step works because NSLSM is constructed to have the right physics (fluid mechanics and surface tension) away from the contact line, and the remaining problem is to improve the algorithm such that it accurately reproduces the physics in the vicinity of the contact line. The value of incorporating contact-line motion into a macroscopic method is that a range of flows can be easily studied and fluid properties, such as viscosity and density, can be varied. A logical next step would be to compare the model against experiments, for instance, for spreading droplets.23 Acknowledgment K.A.S. was supported by NSF Igert Grant 9987577 and by NSF Grants DMR-0076097 and DMR-0109610. Literature Cited (1) Huh, C.; Scriven, L. E. Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. J. Colloid Interface Sci. 1971, 35, 85-101.
(2) De Gennes, P. G. Wettingsstatics and dynamics. Rev. Mod. Phys. 1985, 57, 827-863. (3) Barrat, J. L.; Bocquet, L. Large slip effect at a nonwetting fluid-solid interface. Phys. Rev. Lett. 1999, 82, 4671-4674. (4) Briant, A. J.; Papatzacos, P.; Yeomans, J. M. Lattice Boltzmann simulations of contact line motion in a liquid-gas system. Philos. Trans. R. Soc. London, Ser. A 2002, 360, 485495. (5) Jacqmin, D. Contact-line dynamics of a diffuse fluid interface. J. Fluid Mech. 2000, 402, 57-88. (6) Chen, H. Y.; Jasnow, D.; Vinals, J. Interface and contact line motion in a two phase fluid under shear flow. Phys. Rev. Lett. 2000, 85, 1686-1689. (7) Jones, J. L.; Lal, M.; Ruddock, J. N.; Spenley, N. A. Dynamics of a drop at a liquid/solid interface in simple shear fields: A mesoscopic simulation study. Faraday Discuss. 1999, 112, 129-142. (8) Lowndes, J. The numerical simulation of the steady movement of a fluid meniscus in a capillary tube. J. Fluid Mech. 1980, 101, 631-646. (9) Renardy, M.; Renardy, Y.; Li, J. Numerical simulation of moving contact line problems using a volume-of-fluid method. J. Comput. Phys. 2001, 171, 243-263. (10) Osher, S.; Sethian, J. A. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulation. J. Comput. Phys. 1988, 79, 12-49. (11) Sethian, J. A.; Smereka, P. Level set methods for fluid interfaces. Annu. Rev. Fluid Mech. 2003, 35, 341-372. (12) Sussman, M.; Smereka, P.; Osher, S. A. level set approach for computing solutions to incompressible 2-phase flow. J. Comput. Phys. 1994, 114, 146-159. (13) Smith, K. A.; Solis, F. J.; Tao, T.; Thornton, K.; Olvera de la Cruz, M. Domain growth in ternary fluids: a level set approach. Phys. Rev. Lett. 2000, 84, 91-94. (14) Smith, K. A.; Solis, F. J.; Chopp, D. L. A projection method for motion of triple junctions by level sets. Interface Free Boundaries 2002, 4, 263-276. (15) Smith, K. A.; Ottino, J. M.; Olvera de la Cruz, M. Dynamics of a drop at a fluid interface under shear. Phys. Rev. E 2004, 69, 046302. (16) Brackbill, J. U.; Kothe, D. B.; Zemach, C. A continuum method for modeling surface tension. J. Comput. Phys. 1992, 100, 335-354. (17) Chang, Y. C.; Hou, T. Y.; Merriman, B.; Osher, S. A level set formulation of eulerian interface capturing methods for incompressible fluid flows. J. Comput. Phys. 1996, 124, 449-464. (18) Li, J.; Renardy, Y.; Renardy, M. Numerical simulation of breakup of a viscous drop in simple shear flow through a volumeof-fluid method. Phys. Fluids 2000, 12, 269-282. (19) Chorin, A. J. Numerical Solution of the Navier-Stokes Equations. Math. Comput. 1968, 22, 745. (20) Pomeau, Y. Recent progress in the moving contact line problem: a review. Comptes Rendus Mecanique 2002, 330, 207222. (21) Leong, C. W.; Ottino, J. M. Experiments on mixing due to chaotic advection in a cavity. J. Fluid Mech. 1989, 209, 463-480. (22) Dussan, E. B. Spreading of liquids on solid surfacessstatic and dynamic contact lines. Annu. Rev. Fluid Mech. 1979, 11, 371400. (23) Foister, R. T. The kinetics of displacement wetting in liquid/liquid/solid systems. J. Colloid Interface Sci. 1990, 136, 266-282.
Received for review February 19, 2004 Revised manuscript received June 8, 2004 Accepted June 9, 2004 IE0498605
