9086
Langmuir 2003, 19, 9086-9093
Liquid Drops on Homogeneous and Chemically Heterogeneous Surfaces: A Two-Dimensional Lattice Boltzmann Study Dai Iwahara, Hiroyuki Shinto,* Minoru Miyahara, and Ko Higashitani Department of Chemical Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8510, Japan Received March 17, 2003. In Final Form: July 21, 2003 A boundary condition of a two-phase lattice Boltzmann (LB) method is proposed to model smooth substrates of chemically heterogeneous surfaces. The present technique is applied to the LB simulations of liquid drops on homogeneous, regularly heterogeneous, and randomly heterogeneous surfaces. It is found (i) that the LB technique can mimic not only homogeneous substrates with arbitrary wettabilities but also chemically heterogeneous substrates with arbitrary heterogeneities and wettabilities and (ii) that, in the presence of surface heterogeneities, it also can capture the contact line pinning, the contact angle hysteresis, and the effects of size and distribution of the heterogeneities on them. The results indicate that the present LB method is useful for the simulation of wetting phenomena of chemically heterogeneous surfaces.
1. Introduction Liquid drops spreading on solid surfaces appear in many industrial processes, such as painting and coating, inkjet printing, lubrication, and gluing. To understand this phenomenon, a large number of experimental, theoretical, and numerical studies have been performed (see refs 1-5 and references therein). As a result, the physics of a liquid drop on flat homogeneous surfaces has been relatively well understood,6 while that on heterogeneous surfaces remains in question because of its complexity that originates from realistic situations: surface roughness, chemical contaminations or inhomogeneities in the solid surface, and solutes in the liquid (e.g., surfactants and polymers).1-5 Even with ingenious experimentation, it is still difficult to prepare solid surfaces of controlled structure and chemistry. As computers have become more powerful and less expensive, numerical simulations have turned out to be very useful for modeling liquid drops on well-defined substrates.7 Several molecular dynamics (MD) simulations were carried out to study the spreading of sessile drops on flat homogeneous surfaces,8-11 flat heterogeneous surfaces,12,13 and porous surfaces.14 MD simulations give molecular* To whom correspondence should be addressed. Phone: +81-75-383-2672. Fax: +81-75-383-2652. E-mail: shinto@ cheme.kyoto-u.ac.jp. (1) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces, 6th ed.; Wiley-Interscience: New York, 1997; Chapter X. (2) de Gennes, P. G. Rev. Mod. Phys. 1985, 57, 827. (3) Cazabat, A. M. In Liquids at Interfaces; Charvolin, J., Joanny, J. F., Zinn-Justin, J., Eds.; North-Holland: Amsterdam and New York, 1990; p 371. (4) Leger, L.; Joanny, J. F. Rep. Prog. Phys. 1992, 55, 431. (5) Blake, T. D. In Wettability; Berg, J. C., Ed.; Marcel Dekker: New York, 1993; Chapter 5. (6) De Coninck, J.; de Ruijter, M. J.; Voue´, M. Curr. Opin. Colloid Interface Sci. 2001, 6, 49. (7) Voue´, M.; De Coninck, J. Acta Mater. 2000, 48, 4405. (8) (a) De Coninck, J.; D’Ortona, U.; Koplik, J.; Banavar, J. R. Phys. Rev. Lett. 1995, 74, 928. (b) D’Ortona, U.; De Coninck, J.; Koplik, J.; Banavar, J. R. Phys. Rev. E 1996, 53, 562. (9) Bekink, S.; Karaborni, S.; Verbist, G.; Esselink, K. Phys. Rev. Lett. 1996, 76, 3766. (10) Haataja, M.; Nieminen, J. A.; Ala-Nissila, T. Phys. Rev. E 1996, 53, 5111. (11) de Ruijter, M. J.; Blake, T. D.; De Coninck, J. Langmuir 1999, 15, 7836.
level information on how the structure and chemistry of surfaces influence the wetting. Because of computational cost, however, such MD simulations are unable to access a wide range of parameters (e.g., size and distribution of physical and chemical heterogeneities of substrates, interfacial tension and viscosity of fluids, and gravitation). Another approach is the use of the lattice Boltzmann (LB) method, which is one of the mesoscopic simulation techniques. The LB method reproduces hydrodynamics following these key ideas: a fluid is approximated by a collection of imaginary particles that have a finite number of velocities and collide and propagate on lattices with local conservation of mass and momentum, and the motion of the particles is governed by particle probability distribution functions. For simulating two-phase fluid flows, there are the LB methods based on three different models: the chromodynamic model,15,16 the pseudopotential model,17,18 and the free-energy model.19,20 Wettability of solid substrates has been studied using the chromodynamic,21-24 the pseudopotential,25-29 and the free-energy models;30-33 for example, the spreading of a sessile drop21,23,26,31,33 and its displacement in gravitational29 and flow fields,30 the behavior of a captive (12) Ada˜o, M. H.; de Ruijter, M.; Voue´, M.; De Coninck, J. Phys. Rev. E 1999, 59, 746. (13) Voue´, M.; Semal, S.; De Coninck, J. Langmuir 1999, 15, 7855. (14) Seveno, D.; Ledauphin, V.; Martic, G.; Voue´, M.; De Coninck, J. Langmuir 2002, 18, 7496. (15) Gunstensen, A. K.; Rothman, D. H.; Zaleski, S.; Zanetti, G. Phys. Rev. A 1991, 43, 4320. (16) Grunau, D.; Chen, S.; Eggert, K. Phys. Fluids A 1993, 5, 2557. (17) Shan, X.; Chen, H. Phys. Rev. E 1993, 47, 1815. (18) Shan, X.; Chen, H. Phys. Rev. E 1994, 49, 2941. (19) Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Phys. Rev. Lett. 1995, 75, 830. (20) Swift, M. R.; Orlandini, E.; Osborn, W. R.; Yeomans, J. M. Phys. Rev. E 1996, 54, 5041. (21) D’Ortona, U.; Salin, D.; Cieplak, M.; Rybka, R. B.; Banavar, J. R. Phys. Rev. E 1995, 51, 3718. (22) Blake, T. D.; De Coninck, J.; D’Ortona, U. Langmuir 1995, 11, 4588. (23) van Kats, F. M.; Egberts, P. J. P. J. Colloid Interface Sci. 1998, 205, 166. (24) Hazlett, R. D.; Vaidya, R. N. J. Pet. Sci. Eng. 2002, 33, 161. (25) Martys, N. S.; Chen, H. Phys. Rev. E 1996, 53, 743. (26) Raiskinma¨ki, P.; Koponen, A.; Merikoski, J.; Timonen, J. Comput. Mater. Sci. 2000, 18, 7.
10.1021/la034456g CCC: $25.00 © 2003 American Chemical Society Published on Web 09/03/2003
LB Study of Liquid Drops on Surfaces
drop,21,22,24,25,33 the capillary rise dynamics,27 and the immiscible-fluid displacement in a capillary tube28 and a channel,16,32 where static and dynamic contact angles were evaluated. However, the effects of surface heterogeneities on the wettability have not been reported, except for ref 26, which gives a preliminary result of the surface roughness effects. The aims of the present study are to provide a model of chemically heterogeneous substrates using a two-phase LB method and then to examine their effects on wettability. In this paper, we adopt the two-phase LB method based on the free-energy model, which was recently proposed by Inamuro,34,35 and develop a model of smooth substrates with arbitrary wettabilities and chemical heterogeneities. Using the present LB method, we investigate the physics of liquid drops on homogeneous, regularly heterogeneous, and randomly heterogeneous surfaces.
Langmuir, Vol. 19, No. 21, 2003 9087
The index function φ representing an interface, the flow velocity u, and the pressure p are respectively defined as 9
φ(x,t) ) u(x,t) )
1 gi(x + ci∆x,t + ∆t) ) gi(x,t) - [gi(x,t) - geq i (x,t)] (2) τg eq where feq i and gi are the equilibrium distribution functions, τf and τg are the dimensionless single relaxation times, ∆x is a spacing of the square lattice, and ∆t ≡ (U/c)∆x is a time step during which the particles travel the lattice spacing. The effect of gravity was neglected for simplicity in the present study but can be considered by incorporating an additional term on the right-hand side of eq 2.34
(27) Raiskinma¨ki, P.; Shakib-Manesh, A.; Ja¨sberg, A.; Koponen, A.; Merikoski, J.; Timonen, J. J. Stat. Phys. 2002, 107, 143. (28) Fan, L.; Fang, H.; Lin, Z. Phys. Rev. E 2001, 63, 051603. (29) Kang, Q.; Zhang, D.; Chen, S. Phys. Fluids 2002, 14, 3203. (30) Grubert, D.; Yeomans, J. M. Comput. Phys. Commun. 1999, 121-122, 236. (31) Langaas, K.; Grubert, D. J. Pet. Sci. Eng. 1999, 24, 199. (32) Desplat, J.-C.; Pagonabarraga, I.; Bladon, P. Comput. Phys. Commun. 2001, 134, 273. (33) Briant, A. J.; Papatzacos, P.; Yeomans, J. M. Philos. Trans. R. Soc. London A 2002, 360, 485. (34) Inamuro, T. Bussei Kenkyu 2001, 77, 197. (35) Inamuro, T.; Tomita, R.; Ogino, F. Int. J. Mod. Phys. B 2003, 17, 21. (36) He, X.; Chen, S.; Zhang, R. J. Comput. Phys. 1999, 152, 642.
9
∑cigi(x,t) F(x,t)i)1
(4)
∑
(5)
with the density F
2. Simulation Methods
1 fi(x + ci∆x,t + ∆t) ) fi(x,t) - [fi(x,t) - feq i (x,t)] (1) τf
(3)
1 9 p(x,t) ) [ gi(x,t) - F(x,t)] 3 i)1
F(x,t) ) Let us here consider a liquid drop of density FL in equilibrium with its vapor of density FV. When the liquid drop is placed in contact with a flat solid substrate under negligible gravity, the capillary forces drive the interface toward the equilibrium states. 2.1. Two-Phase Liquid-Vapor Flows. Henceforth, we use nondimensional variables, which are defined by using the characteristic quantities of a length L, a flow speed U, a particle speed c, and a time scale t0 ()L/U), and a reference density F0. The two-dimensional nine-velocity (D2Q9) model was employed, where the physical space of interest is divided into a square lattice and the evolution of the particle population at each lattice node is computed. The D2Q9 model has velocity vectors ci ) (0, 0), (1, 0), (0, 1), (-1, 0), (0, -1), (1, 1), (-1, 1), (-1, -1), and (1, -1) for i ) 1, 2, ..., and 9, respectively. The two-phase LB method by Inamuro34,35 uses two particle distribution functions, fi and gi. In this sense, the method of Inamuro is similar to that of He et al.,36 but it is different from theirs in detail. The distribution functions, fi(x,t) and gi(x,t) with velocity ci at node x and at time t, evolve following
1
fi(x,t) ∑ i)1
φ(x,t) - φV (FL - FV) + FV φL - φV
(6)
where the index function φ ranges from the minimum φV to the maximum φL. The two equilibrium distribution functions are expanded as a power series of the local velocity: 2 feq i ) Hiφ + Fi(p0 - κfφ∇ φ) + 3 9 φ Eiφ 3uRciR - uRuR + uRuβciRciβ + EiκfGRβ ciRciβ (7) 2 2
(
)
geq i ) HiF + 3 9 Ei 3p + F 3uRciR - uRuR + uRuβciRciβ + 2 2 ∂F F + EiGRβ ciRciβ (8) 2FiωguR ∂xR
[
)]
(
where the appropriate coefficients
E1 ) 4/9, E2 ) ... ) E5 ) 1/9, E6 ) ... ) E9 ) 1/36 H1 ) 1, H2 ) H3 ) ... ) H9 ) 0 F1 ) -5/3, Fi ) 3Ei for i ) 2, 3, ..., 9 φ ) GRβ
F ) GRβ
[
9 ∂φ ∂φ 9 ∂φ ∂φ δ 2 ∂xR ∂xβ 4 ∂xγ ∂xγ Rβ
)]
(
}
(9)
(10)
∂F ∂F ∂F ∂F 9 κ + ωg uβ + uR 2 g ∂xR ∂xβ ∂xR ∂xβ 9 ∂F ∂F ∂F κ + 2ωguγ δ (11) 4 g ∂xγ ∂xγ ∂xγ Rβ
(
ωg )
1 1 τ - ∆x 3 g 2
(
)
)
(12)
are determined by local conservation of mass and momentum and by the constraints of Galilean invariance37 and the pressure tensor in an inhomogeneous fluid.38,39 In the above equations, the subscripts R, β, and γ denote Cartesian coordinates (R, β, γ ) x, y) and δRβ is the Kronecker delta. The summation convention is applied to R and β in eqs 7 and 8 and to γ in eqs 10 and 11. For the system without surfaces, the pressure tensor is related to the free-energy functional of the form40 (37) Inamuro, T.; Konishi, N.; Ogino, F. Comput. Phys. Commun. 2000, 129, 32. (38) Yang, A. J. M.; Fleming, P. D., III; Gibbs, J. H. J. Chem. Phys. 1976, 64, 3732. (39) Evans, R. Adv. Phys. 1979, 28, 143.
9088
Langmuir, Vol. 19, No. 21, 2003
ΨB )
∫dx
[
ψ(T,φ) +
Iwahara et al.
]
κf |∇φ|2 2
(13)
where the first term is the bulk free-energy density at a temperature T and the second term is the contribution from any gradients of index functions (or densities via eq 6). In eqs 7 and 13, κf is a constant parameter relevant to the interface width, while κg in eq 11 is one relevant to the interfacial tension.34,35 For a van der Waals fluid, the bulk free-energy density ψ has the form
ψ(T,φ) ) φT ln
(1 -φ bφ) - aφ
2
φT ∂ψ -ψ) - aφ2 ∂φ 1 - bφ
∂φ ∂xR 2
∇ φ≈
≈ 1
1
9
∑ciRφ(x + ci∆x) 6∆x i)2
(16)
9
∑φ(x + ci∆x) - 8φ(x)] 3∆x i)2 [
(17)
The kinematic viscosity ν, the pressure in the interface pint(x), and the interfacial tension γLV are respectively given by
ν)
1 1 τ - ∆x 3 g 2
(
)
pint(x) ) p(x) +
κg |∇F(x)|2 2
∫-∞∞(∂F ∂ξ)
γLV ) κg
2
dξ
FjS ≡
(18) (19) (20)
where ξ is the coordinate normal to the interface.37 The parameters in eqs 14 and 15 were chosen as a ) 1.0, b ) 6.7, and T ) 3.5 × 10-2, which resulted in φL ) 9.714 × 10-2 and φV ) 1.134 × 10-2 because of the requirements for the coexistence chemical potential at φ ) φL and φV.38 The density ratio FL/FV was set at five (i.e., FL ) 0.5 and FV ) 0.1). The other parameters used are κf ) 0.8(∆x)2, κg ) 0.01(∆x)2, τf ) 1.0, and τg ) 1.1 throughout the present simulations, unless specified otherwise. 2.2. Solid Surfaces. To model a solid substrate of a fluid-wettable surface, we assigned fluid densities to the solid-side lattice nodes that are on the solid surface and then applied the halfway bounce-back, no-slip boundary condition41 on the fluid-side lattice nodes that are next to the solid surface, where the boundary lies between two surfaces of the fluid-side and the solid-side nodes. Let us consider a chemically heterogeneous surface of the solid substrate that is represented by a collective of (40) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: Oxford, England, 1982. (41) Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Clarendon Press: Oxford, England, 2001.
(21a)
F2 ) FL - ∆Fd
(21b)
1
nS
∑F(xs,t) ) X1F1 + X2F2
nSs)1
(22)
with
X1 + X2 ) 1
(15)
which is the equation of state of the fluid.38,39 The derivatives in eqs 7, 8, 10, and 11 are calculated using the finite-difference approximations
F1 ) FV + ∆Fd
where ∆Fd ≡ (FL - FV)/2 and is a parameter ranging from 0 to 1 so that FV e F1 e F2 e FL. It is noted that the index function at a solid-side node of given density F1 or F2 is obtained using eq 6. The overall density FjS at the collective of the solid-side nodes xs is given by
(14)
where a and b are the parameters giving the critical temperature as Tc ) 8a/27b and T must be taken less than Tc. Then, p0 in eq 7 is given by
p0 ) φ
the continuous lattice nodes xs (s ) 1, 2, ..., nS) with densities F1 and F2:
(23)
where Xj is the fraction of the solid-side nodes with density Fj (j ) 1, 2). One can consider FjS as the average density of a fluid immediately next to the solid surface. The average affinity of the solid surface for the fluid, χ, is defined by
χ≡
FjS - Fd ) X1χ1 + X2χ2 ) (1 - )(2X2 - 1) (24a) FL - Fd χj ≡
Fj - Fd FL - Fd
(24b)
where Fd ≡ (FL + FV)/2 and χj corresponds to the affinity of the homogeneous surface with Fj (j ) 1, 2). The parameter χ satisfies -1 e χ e 1 and is referred to as the surface affinity in the present study. One can construct a variety of solid surfaces by tuning not only the values of and X2 (or X1) but also the distribution of F1 and F2. Note that when X2 ) 0 or 1, all the solid-side nodes have the same density F1 or F2 and the solid surface is chemically homogeneous. Given the fluid densities at solid-side nodes, the two-phase fluid near the solid surface autonomously evolves to minimize the total free energy of the system given by eq 13. Very recently, the more sophisticated approach32,33 following Cahn theory42 was proposed, where the free-energy functional of eq 13 with the different forms of ψ(T,φ) has an additional surface term which describes the interactions at the surface between the solid and the fluid. The additional term has the form ΨS ) -hφS, where φS is the value of an order parameter (e.g., the density of a one-component fluid, the composition of a fluid mixture, or the index function for them) at the wall and h is a constant parameter controlling the wettability. In contrast to the case of our approach, φS is a dynamical variable; hence, it must be determined every time step using the elaborate schemes.32,33 This approach following Cahn can model a flat homogeneous substrate with arbitrary wettabilities,32,33 but it has not been developed for the chemically heterogeneous substrates, which are our chief interest in the present study. 2.3. Simulation Details. The fluid region was represented by the 128 × 128 lattice nodes at x ) (k∆x, l∆x) for k, l ) 1, 2, ..., 128 (i.e., Lx ) Ly ) 128∆x). The nonwettable substrate of χ ) -1 was located at y ) Ly + ∆x where the solid-side lattice nodes at xs ) (s∆x, Ly + ∆x) for s ) 1, 2, ..., 128 had density FV, while (42) Cahn, J. W. J. Chem. Phys. 1977, 66, 3667.
LB Study of Liquid Drops on Surfaces
Langmuir, Vol. 19, No. 21, 2003 9089
Table 1. Systems of Lattice Boltzmann Simulations system
χ
X2
A.1 A.2 A.3 A.4 A.5 A.6 A.7 A.8 A.9 B.1 B.2 B.3 C.1 C.2 C.3 D.1 D.2 D.3 D.4 D.5 D.6 D.7 D.8 D.9
-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 -0.25 -0.25 -0.25 0.25 0.25 0.25 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
0 0.25 0.5 0.75 1 0.75 0.5 0.25 0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
surface initial shape of no. of heterogeneity liquid phase runs
0 0 0 0 0 or 1 1 1 1 1 0.25 0.25 0.25 0.75 0.75 0.75 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
regular 3-1 regular 6-2 regular 9-3 regular 1-3 regular 2-6 regular 3-9 random random random random random random random random random
circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle, rectangle circle circle circle circle circle circle circle circle circle
1, 1 1, 1 1, 1 1, 1 1, 1 1, 1 1, 1 1, 1 1, 1 4, 4 4, 4 4, 4 4, 4 4, 4 4, 4 3 3 3 3 3 3 3 3 3
the fluid-wettable substrate of surface affinity χ was placed at y ) 0 where the solid-side lattice nodes at xs ) (s∆x,0) for s ) 1, 2, ..., 128 had specified densities. The details of the fluid-wettable substrate will be described in section 4. The halfway bounce-back, no-slip boundary condition41 was imposed on the fluid-side nodes at y ) ∆x and Ly next to the solid-side nodes, while the periodic boundary condition was applied on the sides at x ) ∆x and Lx. The fluid region at time t ) 0 was composed of the liquid phase in contact with the fluid-wettable substrate and the vapor phase elsewhere, and the system was then allowed to evolve for 500 000 time steps. When the liquid phase was taken as a circular drop with radius r0 ) 16∆x, the drop spread to increase the solid/liquid contact area. Alternatively, when it was taken as a film with width 117∆x and height 7∆x, the film shrank to decrease the contact area. In the former and the latter, the advancing and the receding contact angles, θA and θR, were respectively evaluated using the final configuration of the liquid phase. The systems simulated are summarized in Table 1. It was found that under our simulation conditions the systems were almost in equilibrium after 500 000 time steps (see section 4.1). The contact angle of a two-dimensional drop is calculated by
θ ) 2 arctan(2H/D)
(25)
where H is the height of the drop and D is the length of the baseline spanned by two triple points at which vapor, liquid, and solid phases meet. Equation 25 is valid for a nearly semicircular drop, that is, for the case of negligible gravitation.1 The liquid/vapor interface was defined as the contour of the density Fd ) 0.3, and the baseline was taken as that cut out from the line of y ) 2∆x by the interface. It should be noted that the densities at the fluidside nodes of y ) ∆x were strongly affected by those at the solid-side nodes of y ) 0. 3. Preliminary Tests To validate our LB simulations, we first carried out the simulations of a two-dimensional drop of various radii in a vapor. The drop of r ) (10-50)∆x was placed at the
Figure 1. Density profile near the interface of a drop of r ) 20∆x, where ξ is the coordinate normal to the interface. The solid line is a fit of eq 26 with ξd ) 0.4∆x and δ ) 1.8∆x.
Figure 2. Pressure difference between the internal and external regions of a drop as a function of the inverse radius. The solid line is a fit of eq 28 with γLV ) (2.8 × 10-4)∆x.
center of the simulation cell of Lx ) Ly ) 128∆x, and the region excluded from the drop was filled with the vapor phase. The periodic boundary condition was applied on all the sides. The systems were allowed to evolve for 30 000 time steps. 3.1. Density Profiles. Figure 1 shows the density profile near the interface of a drop of r ) 20∆x, where ξ is the coordinate normal to the interface. The liquid/vapor interface has a width of about 10∆x, which is enough resolution to represent the interface. The interfacial tension γLV was computed by eq 20 and found to be γLV ) (2.9 × 10-4)∆x irrespective of radius of the drop. The density profile at a liquid/vapor interface is usually described by
(
Fint(ξ) ) Fd - ∆Fd tanh
)
ξ - ξd δ
(26)
where δ is a measure of the interfacial thickness, ξd is the position of the dividing surface at which Fint(ξd) ) Fd, and again Fd ≡ (FL + FV)/2 and ∆Fd ≡ (FL - FV)/2.40 Equations 20 and 26 yield
γLV )
κg(∆Fd)2 δ
∫-11(1 - χ2int) dχint )
4κg(∆Fd)2 (27) 3δ
where χint ≡ (Fint - Fd)/(FL - Fd) is the normalized density in line with eq 24. Figure 1 illustrates that the simulation data are well fitted by eq 26 with ξd ) 0.4∆x and δ ) 1.8∆x. Equation 27 with the value of δ also gives γLV ) (2.9 × 10-4)∆x. 3.2. Pressure Differences. The internal and external pressures of a drop, pin and pout, were calculated, and the difference ∆p ≡ pin - pout is plotted as a function of r-1 in
9090
Langmuir, Vol. 19, No. 21, 2003
Iwahara et al.
Figure 3. Snapshots of a drop spreading on the homogeneous surfaces of different affinities: (a) χ ) -0.25, system A.4; (b) χ ) 0.25, system A.6. Panels from top to bottom are snapshots at 2000, 10 000, and 500 000 time steps. Dark gray denotes the liquid phase, light gray indicates the vapor phase, and medium gray represents the fluid with intermediate density, where the curved black line indicates the liquid/vapor interface of Fd ) 0.3. The thick black line denotes the solid substrate.
Figure 4. Drop’s base length D as a function of time t on the homogeneous surfaces of different affinities: (4) χ ) -0.25, system A.4; (0) χ ) 0.25, system A.6. The solid lines are a fit of D ∼ tn with n ) 0.28 and 0.35 for systems A.4 and A.6, respectively.
Figure 2, where the relation ∆p ∼ r-1 appears and coincides with Laplace equation for a two-dimensional drop:
∆p ) γLV/r
(28)
The slope of Figure 2 gives γLV ) (2.8 × 10-4)∆x, which is almost equal to the aforementioned value evaluated by eq 20. These results agree with those of other studies using different two-phase LB methods.15,17-19,21,25,26,29 Thus, they clearly demonstrate that the two-phase LB method employed successfully captures the hydrostatics of a drop. 4. Results and Discussion In the three following subsections, we investigate the statics and dynamics of a drop on a smooth solid substrate of a homogeneous surface (system A) or a chemically heterogeneous surface whose heterogeneity is regular (systems B and C) or random (system D), as summarized in Table 1. 4.1. Homogeneous Surfaces. Let us now consider systems A.1-Α.9 in Table 1, where the solid surfaces are homogeneous and have different surface affinities χ defined by eq 24. 4.1.1. Spreading Dynamics. Figure 3 displays typical snapshots of the spreading process of a drop on the surfaces of χ ) -0.25 (system A.4) and χ ) 0.25 (system A.6), indicating that the spreading is faster with increasing the surface affinity χ. For quantitative discussion, Figure 4 shows the evolution of the base length D, which was measured along the line of the solid surface as in section 2.3. The base length exhibits an exponential growth D ∼ tn during a transient stage and then reaches an equilibrium. The exponent n was obtained by fitting the relation to simulation data. The results are given in Figure 5, where n increases with χ, as mentioned above. Tanner’s law for a two-dimensional thin drop predicts the slow growth D ∼ t1/7,43,44 whose exponent is smaller than those in Figure 5. This is because the exponent n strongly depends on the initial conditions of a spreading drop.26 In the present study, the drop initially had a circular shape rather than a flat semicircular shape, for which Tanner’s law is valid. Because of the fast spreading, the simulated drops exhibited deformed semicircular shapes up to the order of a few 104 time steps, as shown in Figure 3. Hence, (43) Tanner, L. H. J. Phys. D 1979, 12, 1473. (44) Safran, S. A. Statistical Thermodynamics of Surfaces, Interfaces, and Membranes; Addison-Wesley: Reading, MA, 1994.
Figure 5. Spreading exponent n as a function of the affinity χ of the homogeneous surface (system A). The solid line is a guide for the eyes.
it is difficult to compare our simulation results with the models of spreading dynamics (e.g., the models in refs 2 and 11), where the spreading drop is assumed to remain semicircular or hemispherical. Also Figure 4 demonstrates that a drop spreading on the surfaces almost reaches an equilibrium state after 500 000 time steps under our simulation conditions. As in the bottom panels of Figure 3, the equilibrium drops evidently exhibit the semicircular shapes in which eq 25 is applicable for evaluating the contact angles. The simulations of system A.6 were carried out for r0 ) 20∆x and 24∆x in addition to r0 ) 16∆x. It was found that neither the exponents of the spreading law nor the equilibrium contact angles depend on the size of the drop, as reported in MD studies of a spreading drop in a complete wetting regime.8a,9,10 4.1.2. Statics. Figure 6 displays the cross-sectional density profiles of the equilibrium fluid for system A.6 along the two lines perpendicular to the substrate (see the bottom panel in part b of Figure 3): x ) 64∆x (i.e., the center line of the drop) and x ) ∆x (i.e., being far from the center line). The curves of eq 26 with δ ) 1.8∆x and the appropriate positions of the dividing surface accurately follow the density profiles for the solid/vapor and the solid/liquid interfaces at the distances of y < 5∆x as well as the liquid/vapor interface at the distances of y ) (14-24)∆x. Hence, the interfacial tensions, γSV and γSL, are obtained in the same manner as γLV of eq 27:
γSV )
κg(∆Fd)2 (2 + 3χ - χ3) 3δ
(29)
LB Study of Liquid Drops on Surfaces
Langmuir, Vol. 19, No. 21, 2003 9091
Figure 6. Cross-sectional density profiles of the liquid drop and the vapor region in contact with the homogeneous surface of χ ) 0.25: system A.6 after 500 000 time steps (see the bottom panel in part b of Figure 3). Filled circles are the densities along the line of x ) 64∆x (i.e., the center line of the drop), and open circles are those along the line of x ) ∆x (i.e., being far from the center line). Solid lines are the fits of eq 26 with δ ) 1.8∆x and the appropriate positions of the dividing surface for solid/vapor, solid/liquid, and liquid/vapor interfaces.
Figure 7. Equilibrium contact angles of a spread drop and a shrunk drop as a function of the affinity χ of the homogeneous surface (system A). The circles denote advancing angles, and the crosses represent receding angles. The solid line represents the expectation by eq 31.
γSL )
κg(∆Fd)2 (2 - 3χ + χ3) 3δ
(30)
The equilibrium contact angle θeq is then
cos θeq )
γSV - γSL χ ) (3 - χ2) γLV 2
(31)
It should be noted that in principle eq 31 is valid for the homogeneous surfaces other than the heterogeneous surfaces, since the validity of eqs 29 and 30 for the latter is not clear. Figure 7 shows the equilibrium contact angles of a spread drop and a shrunk drop, θA and θR, as a function
of χ. We observe the relation θA ) θR ()θeq) regardless of χ, indicating that no contact angle hysteresis takes place when the substrate has a completely smooth homogeneous surface. This result agrees with the conventional postulation.1,2,4,5 Figure 7 demonstrates that the contact angles coincide with the expectations by eq 31. Thus, by changing χ we can represent the homogeneous substrate with arbitrary wettabilities: complete nonwetting (θ ) 180°) by χ ) -1; complete wetting (θ ) 0°) by χ ) 1; partial wetting (180° > θ > 0°) by -1 < χ < 1. As far as equilibrium contact angles are concerned, the results of our approach are very similar to those of the approach following Cahn (see also section 2.2).32,33 4.2. Regularly Heterogeneous Surfaces. Henceforth, we investigate how chemical heterogeneities of a smooth surface affect the drop behavior. First, the surfaces with regular distribution of two chemical groups are considered using systems B.1-B.3 of χ ) -0.25 (FjS ) 0.25) and systems C.1-C.3 of χ ) 0.25 (FjS ) 0.35) in Table 1. Figure 8 illustrates three different units of the bottom solid-side nodes for system C: (C.1) a node of F1 ) 0.2 plus three nodes of F2 ) 0.4; (C.2) the two plus six nodes; (C.3) the three plus nine nodes. In the present study, these surface arrangements of systems C.1-C.3 are respectively referred to as regular 1-3, 2-6, and 3-9 distributions, whereas those of systems B.1-B.3 are defined as regular 3-1, 6-2, and 9-3 distributions, respectively. 4.2.1. Spreading Dynamics. Figure 9 displays typical snapshots of a spreading process of a drop on the surfaces of χ ) 0.25 with 1-3 and 3-9 arrangements (systems C.1 and C.3). The snapshots revealed that the triple points are pinned and immobile within several finite intervals of a transient regime, in which we also observed a stepwise growth of the base length D (data not shown). A comparison between Figure 9 and part b of Figure 3, whose corresponding systems have the same value of χ, suggests that the stick-slip motion of the triple points results from surface heterogeneities and affects the spreading speed and the equilibrium shape of the drop. Obviously, the semicircular shape of the equilibrium drop in part b of Figure 9 for a large heterogeneity differs from those in part a of Figure 9 and part b of Figure 3 for small and zero heterogeneities, respectively. 4.2.2. Statics. The equilibrium contact angles θA and θR for different surface arrangements of χ ) -0.25 and 0.25 are shown in Figure 10. The contact angles of the heterogeneous surfaces lie around that of the corresponding homogeneous surface, θeq. The difference θA - θR remains zero for the smallest size heterogeneity and then increases with increasing size of the surface heterogeneity. Thus, the purely chemical heterogeneity in a smooth surface causes the contact angle hysteresis (θA - θR > 0).
Figure 8. Three different units of chemical heterogeneity patterns of the solid surface. The dark squares represent solid-side nodes with fixed fluid densities of F1 (open circles) and F2 (filled circles). The heterogeneity patterns from top to bottom are referred to as regular 1-3, 2-6, and 3-9 distributions.
9092
Langmuir, Vol. 19, No. 21, 2003
Iwahara et al.
Figure 9. Same as Figure 3 but for the heterogeneous surfaces of χ ) 0.25 (FjS ) 0.35) with different regular distributions of fluid densities of F1 ) 0.2 and F2 ) 0.4: (a) 1-3 distribution, system C.1; (b) 3-9 distribution, system C.3. The short- and long-dashed thick lines denote the solid-side nodes with densities F1 (white) and F2 (black). See also Figure 8.
Figure 11. Cosine of the equilibrium contact angle of a spread drop, cos θA, as a function of the fraction X2 (bottom axis) or the average affinity χ (top axis), for the heterogeneous surfaces with random distributions of fluid densities of F1 ) 0.2 and F2 ) 0.4 (system D). The dashed line represents the prediction by eq 32 with θ1 ) 133° and θ2 ) 47° (i.e., the expectations by eq 31 for systems A.3 and A.7, respectively), and the dotted line is that by eq 33. The solid line denotes the curve of eq 31 for the corresponding homogeneous surfaces of χ ) X1χ1 + X2χ2, where χ1 ) -0.5 and χ2 ) 0.5.
with that of eq 31 for the corresponding homogeneous surfaces of χ ) X1χ1 + X2χ2. Another expression48
(1 + cos θ)2 ) X1(1 + cos θ1)2 + X2(1 + cos θ2)2 (33)
Figure 10. Equilibrium contact angles of a spread/shrunk drop on the heterogeneous surfaces of χ ) -0.25 (system B, triangles) and χ ) 0.25 (system C, circles) with different regular distributions of fluid densities of F1 ) 0.2 and F2 ) 0.4. The open symbols denote advancing angles, and the filled symbols represent receding angles. The dashed lines indicate the expectations by eq 31 for the corresponding homogeneous surfaces (systems A.4 and A.6). 2,45
The results coincide with those of theoretical and numerical studies.46 4.3. Randomly Heterogeneous Surfaces. Following section 4.2, let us now consider the heterogeneous surfaces with random distribution of two chemical groups, using systems D.1-D.9 in Table 1. The fluid density, F1 ) 0.2 or F2 ) 0.4, was randomly given at each of the solid-side nodes such that the fraction of the nodes with density F2 is equal to X2 ) 0.1, 0.2, ..., 0.9 for systems D.1-D.9, respectively. For stable computation, the parameter κf ) 0.5(∆x)2 was employed where the other parameters were the same as those in systems A-C. 4.3.1. Statics. The equilibrium value of cos θA is plotted as a function of X2 in Figure 11. These simulation data lie near the line of the Cassie equation47
cos θ ) X1 cos θ1 + X2 cos θ2
(32)
where X1 and X2 are the fractions of the surface occupied by surface types having contact angles θ1 and θ2, respectively. Notice that the curve of eq 32 almost coincides (45) Joanny, J. F.; de Gennes, P. G. J. Chem. Phys. 1984, 81, 552. (46) Johnson, R. E., Jr.; Dettre, R. H. J. Phys. Chem. 1964, 68, 1744. (47) Cassie, A. B. D. Discuss. Faraday Soc. 1948, 3, 11.
considerably underestimates the variation of contact angle θA with surface fraction X2. The invalidity of eq 33 has been pointed out by experiments49 and theoretical considerations.50 Figures 10 and 11 indicate that eq 32 is suitable for flat chemically heterogeneous surfaces as long as the characteristic size of the heterogeneities is sufficiently smaller than the diameter of the drops, as reported in the MD study.12 Thus, our LB model can capture the wettability of chemically heterogeneous surfaces. 5. Conclusion We have proposed a boundary condition for the twodimensional two-phase LB method of Inamuro34,35 to model smooth substrates of chemically heterogeneous surfaces. Using the present technique, LB simulations of liquid drops on homogeneous, regularly heterogeneous, and randomly heterogeneous surfaces were performed. We found that (i) our LB technique can mimic not only homogeneous substrates with arbitrary wettabilities but also chemically heterogeneous substrates with arbitrary heterogeneities and wettabilities and (ii), in the presence of surface heterogeneities, it also can capture the contact line pinning, the contact angle hysteresis, and the effects of size and distribution of the heterogeneities on them. The results indicate that our LB technique is useful for the simulation of wetting phenomena of chemically heterogeneous surfaces. A three-dimensional two-phase LB study, in which both physically and chemically heterogeneous substrates are considered, will be reported elsewhere.51 (48) Israelachvili, J. N.; Gee, M. L. Langmuir 1989, 5, 288. (49) Laibinis, P. E.; Whitesides, G. M. J. Am. Chem. Soc. 1992, 114, 1990. (50) Johnson, R. E., Jr.; Dettre, R. H. In Wettability; Berg, J. C., Ed.; Marcel Dekker: New York, 1993; Chapter 1. (51) Shinto, H.; Iwahara, D.; Miyahara, M.; Higashitani, K. In preparation for publication.
LB Study of Liquid Drops on Surfaces
Acknowledgment. The authors thank Dr. Takaji Inamuro and co-workers for many valuable discussions. This work was partly supported by the Grantsin-Aid for Scientific Research from the Ministry of
Langmuir, Vol. 19, No. 21, 2003 9093
Education, Culture, Sports, Science and Technology in Japan. LA034456G
