Langmuir 2004, 20, 8137-8141
8137
Lattice Boltzmann Study on the Contact Angle and Contact Line Dynamics of Liquid-Vapor Interfaces Junfeng Zhang and Daniel Y. Kwok* Nanoscale Technology and Engineering Laboratory, Department of Mechanical Engineering, University of Alberta, Edmonton, Alberta T6G 2G8, Canada Received March 18, 2004. In Final Form: June 23, 2004
The moving contact line problem of liquid-vapor interfaces was studied using a mean-field free-energy lattice Boltzmann method recently proposed [Phys. Rev. E 2004, 69, 032602]. We have examined the static and dynamic interfacial behaviors by means of the bubble and capillary wave tests and found that both the Laplace equation of capillarity and the dispersion relation were satisfied. Dynamic contact angles followed the general trend of contact line velocity observed experimentally and can be described by Blake’s theory. The velocity fields near the interface were also obtained and are in good agreement with fluid mechanics and molecular dynamics studies. Our simulations demonstrated that incorporating interfacial effects into the lattice Boltzmann model can be a valuable and powerful alternative in interfacial studies.
I. Introduction The moving contact line (CL) problem is important to a wide range of phenomena and industrial processes such as coating and oil recovery.1,2 Classical fluid dynamics treatment on movement of the CL leads to a stress singularity and a multi-valued velocity field near the CL.3,4 To overcome such difficulties, several theoretical models had been proposed;3-6 these methods, however, typically require the dynamic contact angle (CA) as a boundary condition (BC). Experimental studies have shown that the dynamic CA can deviate significantly from the static value when the CL is in motion.7-11 Yet, the nature of this phenomenon is not well understood. The general behavior of dynamic CAs is that they increase with advancing CL motion and decrease when the CL recedes. To understand this phenomenon, several relations had been presented;4,7,12,13 the fact is that these relations always involve some unknown parameters, which can only be determined by nonlinear fittings with experimental data. On a microscopic level, molecular dynamics (MD) simulations have also been employed to study this problem.14-19 Most of them were performed over a two* To whom correspondence should be addressed. Tel.: (780) 4922791. Fax: (780) 492-2200. E-mail:
[email protected]. URL: http://nanotech.ualberta.ca. (1) de Gennes, P. G. Rev. Mod. Phys. 1985, 57, 827. (2) de Coninck, J.; de Ruijter, M. J.; Voue, M. Curr. Opin. Colloid Interface Sci. 2001, 6, 49. (3) Seppecher, P. Int. J. Eng. Sci. 1996, 34 (9), 977. (4) Sciffer, S. Chem. Eng. Sci. 2000, 55, 5933. (5) Jacqmin, D. J. Fluid Mech. 2000, 402, 57. (6) Fan, L.; Fang, H.; Lin, Z. Phys. Rev. E 2001, 63, 051603. (7) Blake, T. D. In Wettability; Berg, J. C., Ed.; Marcel Dekker: New York, 1993; pp 251-309. (8) Rose, W.; Heinz, R. W. J. Colloid Interface Sci. 1962, 17, 39. (9) Tanner, L. J. Phys. D: Appl. Phys. 1979, 12, 1473. (10) Fermigier, G. M.; Jenffer, P. Ann. Phys. (Paris) 1988, 13, 37. (11) Hoffmann, R. L. J. Colloid Interface Sci. 1975, 50, 228. (12) Cox, R. G. J. Fluid Mech. 1986, 168, 169. (13) Shikhmurzaev, Y. D. Int. J. Multiphase Flow 1993, 19, 589. (14) Koplik, J.; Banavar, J. R.; Willemsen, J. F. Phys. Rev. Lett. 1988, 60 (13), 1282. (15) Thompson, P. A.; Robbins, M. O. Phys. Rev. Lett. 1989, 63 (7), 766. (16) Hadjicostantinou, N. G. Phys. Rev. Lett. 1999, 59 (2), 2475. (17) Koplik, J.; Banavar, J. R. Phys. Rev. Lett. 2000, 84 (19), 4401.
component system, and the velocity fields are usually difficult to resolve.17 As an alternative numerical approach on a mesoscopic level between the traditional computational fluid dynamics and MD methods, the lattice Boltzmann method (LBM) has experienced significant development recently in simulating fluid behaviors.20,21 One of the attractions of the LBM is that the numerical algorithm can be easily implemented with complex solid or free boundaries for multicomponent/multiphase systems. In the bulk fluid, LBM is in fact a Navier-Stokes solver; however, at the fluid-solid interface, the mesoscopic nature of this method becomes manifest as BCs are imposed on the particle distributions rather than directly on the fluid quantities such as velocity.22 For these reasons, the LBM has been widely utilized in interfacial studies.6,23-27 For the moving CL problem, Blake et al.27 and Fan et al.6 had conducted LBM simulations for twocomponent systems to study dynamic CA behaviors. However, no direct relationship between the CA and the CL velocity relative to the substrate was given in ref 27 while the reported changes of the CA in ref 6 were limited and increased by only about 3.5°. In this paper, we present a LBM study on this CL problem by focusing on the dynamic CAs (advancing and receding) and the velocity field near the liquid-vapor interface. We have selected a newly presented mean-field free-energy approach to the LBM for its better representation of fluid-solid interfaces.28 Our results were com(18) Qian, T.; Wang, X.-P.; Sheng, P. Phys. Rev. Lett. 2003, 68, 016306. (19) Gentner, F.; Ogonowski, G.; De Coninck, J. Langmuir 2003, 19, 3996-4003. (20) Succi, S. The Lattice Boltzmann Equation; Oxford University Press: Oxford, 2001. (21) Chen, S.; Doolen, G. D. Annu. Rev. Fluid Mech. 1998, 30, 329. (22) Succi, S. Phys. Rev. Lett. 2002, 89 (6), 064502. (23) Raiskinmaki, P.; Koponen, A.; Merikoski, J.; Timonen, J. Computat. Mater. Sci. 2000, 18, 7. (24) Yang, Z. L.; Dinh, T. N.; Nourgaliev, R. R.; Sehgal, B. R. Int. J. Heat Mass Transfer 2001, 44, 195. (25) Kang, Q.; Zhang, D.; Chen, S. Phys. Fluids 2002, 14 (9), 3203. (26) Raiskinmaki, P.; Shakib-Manesh, A.; Jasberg, A.; Koponen, A.; Merikoski, J.; Timonen, J. J. Stat. Phys. 2002, 107, 143. (27) Blake, T. D.; De Coninck, J.; D’Ortona, U. Langmuir 1995, 11, 4588-4592. (28) Zhang, J.; Li, B.; Kwok, D. Y. Phys. Rev. E 2004, 69, 032602.
10.1021/la049293q CCC: $27.50 © 2004 American Chemical Society Published on Web 08/14/2004
8138
Langmuir, Vol. 20, No. 19, 2004
Zhang and Kwok
pared with those from other theoretical models, and good agreement was found. II. Mean-Field Free-Energy LBM Scheme
of particles at rest.36,37 The macroscopic density F and velocity u can be calculated from the distribution function fi given by
F)
According to the mean-field version of van der Waals’ theory,29-33 the total free-energy function for a fluid system can be expressed as
F)
∫
{
1 dr ψ[F(r)] + F(r) 2
∫
Fu )
}
where ψ(F) is a local free energy with respect to the bulk density F. The second term is a nonlocal term taking into account the free-energy cost of variations in density; φff(r′ - r) is the interaction potential between two fluid particles locating at r′ and r. This term can be reduced to that of a square-gradient approximation when the local density varies slightly.30,34 The third term represents the contribution of the external potential energy V(r) to the free energy F. Both integrations are taken over the entire space. With this expression of free energy, we followed the procedures described in ref 35 and defined a nonlocal pressure as
1 P(r) ) F(r)ψ′[F(r)] - ψ[F(r)] + F(r) 2
∫ dr′ φff(r′ - r) [F(r′) - F(r)] (2)
For a bulk fluid with uniform density, the nonlocal integral term disappears and eq 2 reverts to the equation of state of the fluid
P ) Fψ′(F) - ψ(F)
(3)
Here, we describe the implementation of these results into a LBM algorithm. In general, after discretization in time and space, the lattice Boltzmann equation with a Bhatnagar-Gross-Krook collision term can be written as
1 fi(x + ei, t + 1) - fi(x, t) ) - [fi(x, t) - feq i (x, t)] (4) τ where the distribution function fi(x, t) denotes the particle population moving in the direction of ei at a lattice site x and at a time step t; τ is a relaxation time and feq i (x, t) is a prescribed equilibrium distribution function with
[
]
1 2 u , c2 1 - d0 D(D + 2) D + 2 ei‚u + (ei‚u)2 feq i ) F l cl 2c4l
[
(6)
∑i fiei
(7)
and
dr′ φff(r′ - r)[F(r′) F(r)] + F(r) V(r) (1)
feq 0 ) F d0 -
∑i fi
]
However, if an external force F(x, t) exists, we can modify eq 7 to reflect the momentum change as
Fu )
∑i fiei + τF
(8)
and employ the u produced here to calculate the equi37,38 Redefining librium distribution function feq i in eq 5. the fluid momentum Fv to be an average of the momentum before collision ∑ifiei and that after ∑ifiei + F yields
Fv )
1
∑i fiei + 2F
(9)
Following the Chapman-Enskog procedure, a NavierStokes equation with the equation of state
P)
c2(1 - d0) F+Φ D
(10)
can be obtained, where Φ is the potential energy field related to F by
F(x, t) ) -∇Φ(x, t)
(11)
To obtain the Navier-Stokes equation with a pressure term similar to that given by eq 2, we set an artificial Φ as follows:
∫
1 Φ(x, t) ) F(x)ψ′[F(x)] - ψ[F(x)] + F(x) dx′ φff(x′ 2 c2(1 - d0) F(x) (12) x)[F(x′) - F(x)] D The above equations set up a complete LBM scheme with the mean-field free-energy function implemented. III. Examinations of Liquid-Vapor Interface Behaviors
for a D-dimensional lattice with l links where the particle speed is |ei| ) c. The constant d0 is the equilibrium fraction
Before applying the above method to study the moving CL problem, examination of the liquid-vapor interface behaviors is necessary. In this section, we perform bubble tests and capillary wave tests, which have been commonly employed by different multiphase/multicomponent LBM models.36,37,39-41 The bubble tests were carried out over a
(29) Rowlinson, J.; Widom, B. Molecular Theory of Capillarity; Claredon: Oxford, 1982. (30) Sullivan, D. E. J. Chem. Phys. 1981, 74 (4), 2604. (31) van Giessen, A. E.; Bukman, D. J.; Widom, B. J. Colloid Interface Sci. 1997, 192, 257. (32) Zhang, J.; Kwok, D. Y. J. Phys. Chem. B 2002, 106, 12594. (33) Zhang, J.; Kwok, D. Y. Langmuir 2003, 19, 4666. (34) Widom, B. J. Stat. Phys. 1978, 19, 563. (35) Yang, A. J. M.; Fleming, P. D.; Gibbs, J. H. J. Chem. Phys. 1976, 64 (9), 3732.
(36) Shan, X.; Chen, H. Phys. Rev. Lett. 1993, 47, 1815. (37) Shan, X.; Chen, H. Phys. Rev. Lett. 1994, 49, 2941. (38) Buick, J. M. Lattice Boltzmann Methods in Interfacial Wave Modelling. Ph.D. Thesis, The University of Edinburgh, Edinburgh, U.K., 1997. (39) Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Phys. Rev. Lett. 1995, 75 (5), 830. (40) Swift, M. R.; Orlandmi, E.; Osborn, W. R.; Yeomans, J. M. Phys. Rev. Lett. 1996, 54 (5), 5041. (41) Appert, C.; Rothman, D. H.; Zaleski, S. Physica D 1991, 47, 85.
D 2 u , i ) 1, ..., l (5) 2lc2
CA and CL Dynamics of Liquid-Vapor Interfaces
Langmuir, Vol. 20, No. 19, 2004 8139
Figure 1. Pressure differences across the interface for different bubble sizes from LBM simulations (symbols). The solid line is a curve fit according to Laplace’s equation, eq 15.
Figure 2. Dispersion relation of capillary waves on a liquidvapor interface. The symbols are LBM results while the solid line is a curve fit according to fluid mechanics theory, eq 16.
128 × 128 D2Q7 lattice domain with periodic BCs in both directions, having D ) 2, l ) 6, d0 ) 1/2, c ) 1, e0 ) [0, 0], and ei ) [cos(i - 1)π/3, sin(i - 1)π/3] (i ) 1, ..., 6) in eq 5. Following refs 28, 39, 40, and 42, we employed a van der Waals fluid model to express the free-energy of the bulk fluid:
ψ(F) ) FkBT ln
F - aF2 1 - bF
(13)
where a and b are the van der Waals constants, kB is the Boltzmann constant, and T is the absolute temperature. Here, we have selected a ) 9/49, b ) 2/21, and the scaled temperature kBT ) 0.52. The fluid-fluid interaction was simplified by considering only that between nearest neighbors:36,37
{
K, |x′ - x| ) c φff(x′ - x) ) 0, |x′ - x| * c
We also calculated the dispersion relation of capillary waves. In the absence of gravity, this relation for a fluidfluid interface is given by fluid dynamics as
ω2 ∼ k3 (14)
which K measures the interaction strength among the nearest neighboring particles and the nonlocal integral term in eq 12 can be replaced by a summation over the neighbors of a site x. In this work, K was set to -0.01. A negative value implies an attraction between two fluid particles, and its magnitude changes the liquid-vapor interfacial tension. The initial density distribution was assigned with a small higher density region in the center. After about 10 000 steps, the system reached equilibrium and the drop radius and densities inside and outside were measured. The pressure difference across the interface can be obtained through the equation of state eq 3. The results were plotted in Figure 1 together with a fitted line according to Laplace’s equation
Pin - Pout ) γ/R
Figure 3. Schematic diagram (not to scale) of the simulation framework. θa and θr in the figure are the advancing CA with contact line velocity U and the receding CA with -U, respectively.
(16)
where ω is frequency and k is the wavenumber of the capillary wave.39,43 Our simulations were performed over W × H rectangular lattice domains with periodic BCs in the horizontal direction and bounce-back BCs on the top and bottom layers. The height H was selected as 10W to overcome the effects of finite fluid depth. Initially, the density in the upper part of the domain was set to be higher while that of the lower was set to be less. After equilibrium has been established, a planar liquid-vapor interface was formed in the middle and we introduced a single sinusoidal disturbance with a wavelength W on the interface. The wave amplitude was thereafter measured to obtain the frequency ω and wavenumber k ) 2π/W. Our results (symbols) are displayed in Figure 2 with a nonlinear fitting (solid line) of ω ) Rkβ. The fitted line provides a value of β as 1.48 and is nearly identical to the theoretical value of 3/2 in eq 16, validating the dynamics of our LBM for liquid-vapor interfaces.
(15)
IV. Dynamic CAs and Velocity Fields around CLs
where γ is the interfacial tension. Clearly, the symbols from our LBM results follow the Laplace equation well, and the interfacial tension γ was found to be 0.026 from the slope of the fitted line. By direct integration of the excess free energy across a planar interface, the interfacial tension had been found to be 0.029,28 which is in principle similar to the value obtained here by bubble test simulations.
Similar to the procedures described in ref 6, simulations of moving CLs were conducted on a two-dimensional channel. By taking advantage of symmetry about the centerline, only half of the channel was considered here with a mirror BC applied along the center line (Figure 3). The simulation domain was 1000 × 20 with periodic BCs at the left and right ends. A general bounce-back no-slip BC was applied along the channel wall. The solid-fluid interaction was also simplified as described above and
(42) Angelopoulos, A. D.; Paunov, V. N.; Burganos, V. N.; Payatakes, A. C. Phys. Rev. Lett. 1998, 57 (3), 3237.
(43) Landau, L. D.; Lifshitz, E. M. Fluid Mechanics; Pergamon Press: New York, 1959.
8140
Langmuir, Vol. 20, No. 19, 2004
Zhang and Kwok
Figure 4. Dynamic CAs obtained at different CL velocities (symbols), with nonlinear fittings (solid lines) according to Blake’s theory, eq 18.
was implemented as an attraction force between a solid (xs) and a fluid (xf) site as
Fs )
{
KwF(xf)(xs - xf), |xs - xf| ) c |xs - xf| * c 0,
(17)
Figure 5. Flow fields near the advancing (a) and receding (b) interfaces. The curves are constant density contours.
where θ0 is the static CA (U ) 0). Others are parameters related to the adsorption/desorption process: n is the number of adsorption sites per unit area, λ is the average length of each particle displacement, and κw0 is the
equilibrium frequency of particle displacements. Typically, these parameters can only be determined by curve-fitting.7 It can be seen in Figure 4 that our simulation results follow Blake’s relation eq 18 very closely, except for large velocities. The difference may be caused by the fact that, when the velocity becomes large, the liquid-vapor interface shape deviates from a circular shape and the CA obtained from the geometry relation becomes less accurate. However, manual measurements of such angles are difficult to be consistent and could cause a much larger uncertainty. Nevertheless, our LBM results demonstrate the general behaviors of dynamic CAs. The fitting parameter values are θ0 ) 126.0°, 2κw0λ ) 12.62, and 2nkBT/γ ) 0.23 for the Kw ) 0.04 surface and θ0 ) 55.1°, 2κw0λ ) 52.06, and 2nkBT/γ ) 0.67 for the Kw ) 0.08 surface; the latter implies more adsorption sites and a higher displacement frequency on a surface with stronger attractions to the fluid particles. Another aspect that we were interested in is the flow in the interfacial region. Figure 5 displays the velocity fields around the advancing and receding ends. The hydrophobic channel wall was modeled using Kw ) 0.04 with an interfacial velocity of U ) 0.014. Focusing on the three-phase contact points, the flow fields are similar to that of the picture using Blake’s adsorption/desorption process: when fluid 1 moves forward to displace fluid 2 (the two fluids may be either both liquids or a liquid and a vapor), the fluid 2 particles adsorbed on the solid surface would be pushed away (desorption) and the fluid 1 particles would be attached onto the surface site left behind by the particles of fluid 2 (adsorption).7 In moving CL studies, the frame of reference is typically selected for a stationary CL and an interface.3,5,16,44,46 Thus, we plotted the flow streamlines (from Figure 5) in Figure 6 with the contours of constant density (circular curves in the center) in such a reference system. However, because the flow field near interfaces is complicated, errors can be introduced when integrating fluid velocities across the interfaces to produce streamlines. We note that, far away from the interfaces, the flow displays the same features as the Poiseuille flow with a parabolic velocity profile.
(44) Sheng, P.; Zhou, M. Phys. Rev. A 1992, 45 (8), 5694. (45) Blake, T. D.; Haynes, J. M. J. Colloid Interface Sci. 1969, 30, 421.
(46) Hadjiconstantinou, N. G.; Patera, A. T. Int. J. Numer. Methods Fluids 2000, 34, 711.
The coefficient Kw is positive for attraction forces, and by adjusting its magnitude, solid surfaces with different wettability (CAs) can be easily modeled, from wetting to nonwetting.28 In this paper, two kinds of channel surfaces, one hydrophobic (Kw ) 0.04) and another hydrophilic (Kw ) 0.08), were studied. Here, we label a hydrophobic surface to be low energy (CA > 90°) and the hydrophilic one to be high energy (CA < 90°). The density distribution in the channel at the beginning was set to the density of the liquid phase, except for a segment in the center to that of the vapor phase, simulating a vapor bubble enclosed in the channel by two liquid columns with two liquid-vapor interfaces. The densities of the liquid and vapor phases were obtained as described in the previous section. To cause the interfaces to move, we applied a body force on the fluid particles in the direction along the channel and various velocities can be easily obtained by adjusting the force magnitude. Such velocities were calculated by measuring the interface positions at different time steps. The dynamic CAs were also measured at both ends of the enclosed bubble, corresponding to the advancing and receding CAs, respectively. Our CA calculation utilizes a geometry relation which is the same as that given in refs 6 and 44, which assumed the liquid-vapor interface as a circular arc. The liquid-vapor interface was treated as the isoline of the average density. The dynamic CAs at different velocities are displayed in Figure 4 where the solid lines were determined by nonlinear fittings using Blake’s adsorption/desorption model due to its good agreement with experimental results.7,45 According to Blake’s theory, The relation between the dynamic CA θ and CL velocity U is given as
U ) 2κw0λ sinh
[
]
γ(cos θ0 - cos θ) 2nkBT
(18)
CA and CL Dynamics of Liquid-Vapor Interfaces
Langmuir, Vol. 20, No. 19, 2004 8141
transfer was explained as an evaporation/condensation process to overcome the CL velocity singularity. V. Conclusions
Figure 6. Streamlines and density contours (circular curves in the center) in the vicinities of the advancing (a) and receding (b) interfaces. The reference frame is moving at a constant velocity with the interfaces.
The general streamline patterns near the interfaces are, however, quite similar to the results from fluid dynamics3,5,44,46 and MD simulation.16 The flow agrees in character with the wedge flow solution.47 We also observed small vortexes near the CL and beside the interface found in refs 3 and 44. The reason of the mass transfer across the interface at the advancing end is not clear. It may be due to spurious current in the interfacial region,48 fluid diffusion reported in ref 3, or simply the integration errors mentioned above. More recently, Briant et al.49 studied droplet deformation in a shear flow field and the mass (47) Huh, E.; Scriven, L. E. J. Colloid Interface Sci. 1971, 35, 85. (48) Hou, S.; Shan, X.; Zou, Q.; Doolen, G. D.; Soll, W. E. J. Comput. Phys. 1997, 138, 695.
In this paper, we have examined the liquid-vapor interface behaviors from a recently proposed mean-field free-energy LBM model. Results from our LBM satisfies the Laplace law and dispersion relation of capillary waves. We also studied the dynamic advancing/receding CAs at different velocities. The simulation results agree well with that from Blake’s theory. The flow field near the threephase contact points displays a similar picture of the adsorption/desorption mechanism of Blake’s model. The general flow patterns also follow a wedge flow and are consistent with the results from hydrodynamics and MD studies. In reality, the moving CL problem is complex and its physics is still not well understood. The fluid properties as well as the slipping over a solid-liquid interface may all play a crucial role. Nevertheless, our work has shown that our new LBM approach for solid-liquid interfaces could be a valuable alternative besides fluid mechanics and MD simulations. Acknowledgment. This work was supported, in part, by the Alberta Ingenuity Establishment Fund, Canada Research Chair (CRC) Program, Canada Foundation for Innovation (CFI), and Natural Sciences and Engineering Research Council of Canada (NSERC). J.Z. acknowledges the financial support from the Alberta Ingenuity through a studentship fund. LA049293Q (49) Briant, A. J.; Wagner, A. J.; Yeomans, J. M. Phys. Rev. E 2004, 69, 031602.
