Wetting Dynamics of Drop Spreading. New Evidence for the

pubs.acs.org/Langmuir © 2010 American Chemical Society

Wetting Dynamics of Drop Spreading. New Evidence for the Microscopic Validity of the Molecular-Kinetic Theory D. Seveno,* N. Dinter, and J. De Coninck Laboratory of Surface and Interfacial Physics (LPSI), Universit
e de Mons, Place du Parc, 20. 7000 Mons, Belgium Received June 14, 2010. Revised Manuscript Received July 27, 2010 We study the spontaneous wetting of liquid drops on FCC solid substrates using large-scale molecular dynamics simulations. By varying the solid lattice parameter, five different drop/solid dynamic systems are investigated. It is shown that the results are in agreement with the molecular-kinetic theory (MKT) describing the dynamics of wetting. Moreover, it is established that the microscopic parameters resulting from fits using the MKT, the so-called molecular jump frequency at equilibrium and the jump length, correspond to the values that can be estimated directly from the simulations. This agreement strongly supports the validity of the MKT at the microscopic scale.

I. Introduction Wetting dynamics has been the subject of much research over the past few decades, both theoretically and experimentally.1-11 In particular, the wetting of liquid drops on solid surfaces has been the focus of much attention. Whenever a liquid drop touches a solid surface, it starts to wet. Since the drop shape changes with time, some energy has to be dissipated by the system. Several channels of dissipation can be considered: viscous dissipation in the wedge of the contact angle, frictional processes at the threephase contact line, and so on. To be able to predict the dynamics of wetting, we have to identify the correct dissipation channels to describe the link between the wetting speed v(t) and the dynamic contact angle θ(t): v(t) as a function of θ(t). Usually, one is limited to fitting data to one model or another. This approach sometimes leads to strange values for the fitted parameters such as fractions of angstroms for an atomic cutoff length or 20 nm for the average length between adsorption sites on a solid surface. We suggest that these nonphysical values result from fitting data with inappropriate theories through incomplete understanding of the mechanisms at work. However, this claim makes sense only if we can first prove that the theories themselves are meaningful at the appropriate scale. Among the few existing theories, there is one which can be tested directly, without fitting, by molecular dynamics (MD) simulations. Following the work of Blake and Haynes,2-4 the dynamics of wetting can be described via the molecular-kinetic theory (MKT), which expresses the triple-line (TL) velocity v(t) as a function of the dynamic contact angle θ(t), via equilibrium parameters: K0, λ, and θ0, being, respectively, the molecular jump (1) de Gennes, P.-G. Rev. Mod. Phys. 1985, 57, 827. (2) Blake, T. D.; Haynes., J. M. J. Colloid Interface Sci. 1969, 30, 421. (3) Blake, T. D. J. Colloid Interface Sci. 2006, 299, 1. (4) Blake, T.D. In Wettability; Berg, J., Ed.; Marcel Decker: New York, 1993; Chapter 5. (5) Blake, T. D.; De Coninck, J. Adv. Colloid Interface Sci. 2002, 96, 21. (6) Hayes, R.; Ralston, J. Colloids Surf., A 1994, 93, 15. (7) De Coninck, J.; de Ruijter, M.; Vou
e, M. Curr. Opin. Colloid Interface Sci. 2001, 6, 49. (8) Bonn, D.; Egers, J.; Indekeu, J.; Meunier, J.; Rolley, E. Rev. Modern Phys. 2009, 81, 739. (9) Brochart-Wyart, F.; de Gennes, P. G. Adv. Colloid Interface Sci. 1992, 39, 1. (10) de Ruijter, M.; De Coninck, J.; Oshanin, G. Langmuir 1999, 15, 2209. (11) Petrov, P; Petrov, J Langmuir 1992, 8, 1762.

14642 DOI: 10.1021/la102412s

frequency, the jump length, and the equilibrium contact angle. The relation between these parameters and v(t) is given by eq 1, with n being the number of interaction sites per unit area of solid surface (usually, n is taken as λ-2), γ the liquid surface tension, kB the Boltzmann constant, and T the temperature: ! γðcos θ0 - cos θðtÞÞ 0 ð1Þ vðtÞ ¼ 2K λ sinh 2nkB T The remarkable aspect of this equation is that the dynamics of wetting is totally controlled by the equilibrium properties K 0, λ, and θ0. The displacement of the TL is a result of the collective motion of the liquid molecules and their component atoms. Atoms from the layer in contact with the solid surface move perpendicular and parallel to the surface, leaving holes. These holes are repopulated by other atoms. The energy lost in this process is considered to be the source of contact-line friction. This equation has been shown to describe rather well several sets of data.12-28 Various MD simulations (refs 15, 29-31 and references (12) Saiz, E.; Tomsia, T. Nat. Mater. 2004, 3, 903. (13) Saiz, E.; Cannon, E.; Tomsia, A. Ann. Rev. Mater. Res. 2008, 38, 197. (14) Webb, E. B., III; Hoyt, J. J.; Grest., G. S. Curr. Opin. Solid State Mater. Sci. 2005, 9, 174. (15) de Ruijter, M. J.; Blake, T. D.; De Coninck, J. Langmuir 1999, 15, 7836. (16) de Ruijter, M. J.; De Coninck, J.; Blake, T. D.; Clarke, A.; Rankin, A. Langmuir 1997, 13, 7293. (17) de Ruijter, M. J.; Charlot, M.; Vou
e, M.; De Coninck, J. Langmuir 2000, 16, 2363. (18) Petrov, J. G.; Ralston, J.; Schneemilch, M.; Hayes, R. J. Phys. Chem. B 2003, 107, 2003. (19) Brooks, C. F.; Grillet, A. M.; Emerson, J. A. Langmuir 2006, 22, 9928. (20) Vega, M. S.; Seveno, D.; Lemaur, G.; Ad~ao, M. H.; De Coninck, J. Langmuir 2005, 21, 9584. (21) Reddy, S.; Schunk, P. R.; Bonnecaze, R. T. Phys. Fluids 2005, 17, 122104. (22) Phan, C. M.; Nguyen, A. V.; Evans, G. M. Langmuir 2003, 19, 6796. (23) Le Grand, N.; Daerr, A.; Limat, L. J. Fluid Mech. 2005, 541, 293. (24) Mc Hale, G; Newton, M. I. Colloids Surf., A 2002, 206, 193. (25) Ray, S.; Sedev, R.; Priest, C.; Ralston, J. Langmuir 2008, 24, 13007. (26) Ranabothu, S. R.; Karnezis, C.; Dai, L. L. J. Colloid Interface Sci. 2005, 288, 213. (27) Stalcup, E. J.; Seemann, R.; Herminghaus, S.; Law, B. M. J. Colloid Interface Sci. 2009, 338, 523. (28) Prevost, A.; Rolley, E.; Guthmann, C. Phys. Rev. Lett. 1999, 83, 348. (29) De Coninck, J.; Blake, T. Ann. Rev. Mater. Res. 2008, 38, 1. (30) Bertrand, E.; Blake, T. D.; De Coninck, J. J. Phys.: Condens. Matter 2009, 21, 464124. (31) Benhassine, M.; Saiz, E.; Tomsia, A. P.; De Coninck, J. Langmuir 2009, 25, 11450.

Published on Web 08/18/2010

Langmuir 2010, 26(18), 14642–14647

Seveno et al.

Article

Table 1. FCC Parameters Used in This Study relative lattice minimum interatomic density parameter (A˚) distance (A˚) 1 0.75 0.7 0.6 0.5

5.55 6.11 6.25 6.59 7

atomic density (10-2atom/A˚3)

3.93 4.32 4.42 4.66 4.95

2.332 1.749 1.632 1.399 1.166

therein) have shown that the perpendicular jump frequencies are smaller than the horizontal ones due to the attraction of the solid surface. The limiting factor for repopulation is thus the perpendicular jump frequency. These have been shown to be compatible with the values of K 0 obtained by fitting eq 1 to several sets of MD data.15,29-31 However, the length of the jump λ has usually been assumed to equal to the average distance between adsorption sites on the solid surface. In the one case30 where the value was directly measured, it was shown that this assumption was broadly correct. This point is rather critical in eq 1 since one has to obtain the parameters 1/n and 2K0λ by fitting the data and assuming n = λ-2. A bad value for λ will thus affect the estimate of K 0. The object of this Article is to use MD to explore this crucial aspect of the problem in detail by investigating several solid surfaces using the same liquid. By changing the substrate geometry but not the nature of the atoms (which is only possible using simulations), we will show that the equilibrium jump length indeed corresponds closely to the distance between adsorption sites on the solid. We also tested our data upon the hydrodynamic (HD) approach describing the viscous dissipation1 and predicting that " #   L  -1 γ 3 3 θt - θ0 ln vðtÞ ¼ 9η Ls

ð2Þ

with η being the liquid viscosity (Pa 3 s), L a characteristic length scale of the drop, and Ls the slip length. There is only one free parameter which is the logarithm of the ratio of these lengths. Its value is usually found to be about 10-15 as L is on the order of a millimeter and Ls of a molecular scale.

II. Simulation Details The systems we consider are two-dimensional periodic cells extending in the surface directions containing FCC substrates and the liquid droplet. The substrates are orientated with the [100] axis along the z-axis. They are composed of 40 500 atoms with varying lattice parameters. All the parameters remain constant for the simulations except the FCC lattice parameter. These take the values 5.55, 6.11, 6.25, 6.59, and 7 A˚. In terms of relative solid density, if d denotes the ratio of substrate volume to the most close packed FCC, we get, respectively, d = 1, d = 0.75, d = 0.7, d = 0.6, and d = 0.5. Table 1 lists the FCC solid parameters. Liquid droplets are composed of 25 600 atoms of 16 carbon-like chains. The surface expanses prevent interaction between the left part of the drop and its right part; in addition, the substrates are thick enough to ensure that atoms located in the bottom layer do not influence drop wetting. A first equilibrium step is conducted with a thermostat (velocity rescaling) to thermalize the drop and to ensure that the interatomic distances within the drop are equilibrated. The drop is then dragged close enough to the surface (5 A˚) to be influenced by the surface and to wet during the MD simulations. In order to simulate realistic experiments, the solid temperature is fixed via the thermostat whereas the temperature of the drop is free to evolve during the wetting dynamics. Langmuir 2010, 26(18), 14642–14647

Figure 1. Contour map of the potential energy due to Lennard-Jones interactions for a single liquid atom interacting with the solid at distance σ from the solid for d = 1. The white areas show the position of the solid atoms and the dark areas the location of the adsorption sites.

The potentials between the atoms, both solid and liquid, are described by standard pairwise Lennard-Jones 12-6 interactions: Vij ðrÞ ¼ 4εij

 12  6 ! σ ij σ ij Cij - Dij r r

ð3Þ

where r is the distance between any pair of atoms i and j. The parameters εij and σij are related, respectively, to the depth of the potential well and an effective molecular diameter. For both solid and liquid atoms, the Lennard-Jones parameters are εij = 0.267  103 J/mol and σij = 3.5 A˚. The pair potential is set to zero for r g 2.5 σij (the cutoff length). The repulsive and attractive coupling parameters Cij and Dij are set equal, with a value of 1.0 for solid/ solid and liquid/liquid interactions and 0.9 for the solid/liquid interaction, ensuring partial wetting.29 This model is rather simplistic, but it contains all the basic ingredients to describe the details of wetting.15,29,30 In addition, we consider a confining potential between nearest neighbors to maintain a constant distance between any two adjacent atoms within a given molecule potential, V conf ¼ Dconf r6

ð4Þ

where Dconf = 1.0 and r is the distance of separation. Moreover, the solid atoms vibrate thermally around their initial equilibrium position according to a harmonic potential. The MD time step is fixed to 5 fs, and the temperature is fixed by the equality T = kBεij (= 33.33 K) by a velocity scaling scheme. Figure 1 illustrates the energy landscape experienced by a liquid atom at a distance σ from the solid. The sites where the attraction is the greatest correspond to the energy well where a liquid atom can jump (minimum interatomic distance). Two types of jumps are therefore possible and should dominate: first between two closest sites and then between two sites in diagonal. We assumed that the ratio between these types of jump can vary from 1:2 to 2:1 so that the jump lengths induced by the geometry of the solid surface should be λm = 4.75 ( 0.16, 5.22 ( 0.18, 5.34 ( 0.18, 5.63 ( 0.19, and 5.98 ( 0.20 A˚ for, respectively, d = 1, 0.75, 0.7, 0.6, and 0.5.

III. Measurements Let us first describe the equilibrium situation of our systems. The system of molecules, even at equilibrium, is subjected to DOI: 10.1021/la102412s

14643

Article

Figure 2. Atom density versus the distance from the solid: (0) d = 1, (O) d = 0.75, (4) d = 0.7, (3) d = 0.6, and (g) d = 0.5.

fluctuations in its density. However, if we average in time, the density becomes reproducible. In Figure 2, the vertical atom density profiles are given at equilibrium. Close to the solid, we see strong changes in the profile because of the presence of the substrate, while at larger distance the atom density reaches a constant value corresponding to the liquid bulk (see inset). Indeed, near the solid, the liquid is much more ordered, and a few layers can be distinguished. This is in agreement with previous simulations.13-17 At a distance larger than 4 layers away from the substrate, the liquid density is uniform (Figure 2, inset) corresponding to approximately 15 A˚. To describe the dynamics of wetting, we need to know the speed and the contact angle versus time. To achieve this, we record the position of the edge of the drop throughout the simulation. As usual, the atom density parallel to the surface is used to locate this edge. First, we subdivide the liquid droplet into many horizontal sheets of arbitrary thickness. The constraint on the number of sheets is provided by the need to maximize the number of sheets while ensuring that each sheet contains enough atoms to give a uniform density. For each sheet, we locate its center by symmetry and compute the density of particles as a function of the distance from the center. We then locate the extremity of the layer as the distance where the density falls below some cutoff value (half the liquid bulk density). To check the consistency of the method, different sheet thicknesses and cutoff values were considered and were found to give almost identical results. These methods enable us to construct the complete profile of the drop and to determine how it evolves with time. We approximate the drop shape by a spherical cap, and the best circular fits are determined throughout the simulations; see Figure 3. At the very start of the simulation, when the drop just comes into contact with the surface, the spherical cap approximation is not valid, so these configurations are omitted from our data. However, after short times, corresponding typically to the first 1.0 ns, the spherical cap gives a good approximation of the drop shape as illustrated in Figure 4 where the standard deviation fit 2 1/2 ([ΣN i=1(zi - zi ) /N ] , where N is the number of data points in the droplet profile, zi is a point of the profile, and zfit i is its associated point obtained from the best circular fit) between the profile of the droplet and the best circular fit is plotted versus time for d = 0.75. The simulated drops subsequently retain their spherical form during wetting, except for the first few layers at the liquid/solid interface. 14644 DOI: 10.1021/la102412s

Seveno et al.

Figure 3. Spherical cap approximation for a sessile droplet (in gray) on a FCC substrate (black dots).

Figure 4. Standard deviation between the profile of the drop and the best circular fit versus time for d = 0.75.

Figure 5. Dynamic contact angle measured during the MD simulations: (0) d = 1, (O) d = 0.75, (4) d = 0.7, (3) d = 0.6, and (g) d = 0.5.

Here, we expect the profile to be perturbed for energetic and entropic reasons.32 To avoid this problem, we investigated the profile as a function of the number of layers used, from top to bottom. The tangent to the circular profile at the intersection point with the substrate defines the contact angle. Using this (32) De Coninck, J.; Dunlop, F.; Menu, F. Phys. Rev. E 1993, 47, 1820.

Langmuir 2010, 26(18), 14642–14647

Seveno et al.

procedure, we are able to measure the contact angle θ(t) as a function of the number of time steps during our simulations. Typical results are presented in Figure 5 for the five different solids. For clarity, the error bars (typically 5
) are not shown. Since the volume of the drop is constant, we can easily measure the variation of the drop radius r(t). To determine the speed of the three-phase contact line, we use a Levenberg-Marquardt procedure to fit the radius to a robust function.33 The function used is a ratio of polynomials of order less or equal to 10, which can

Figure 6. TL velocity versus dynamic contact angle: (0) d = 1, (O) d = 0.75, (4) d = 0.7, (3) d = 0.6, and (g) d = 0.5. The lines are the best fits found by G-Dyna.

Article

accurately reproduce the observed evolution of the contact radius for the range of times considered here, n P

rðtÞ ¼

i¼0

1þ

ai ti

n P

i¼1

ð5Þ bi t i

The order of the polynomial has to be fixed to maximize the quality of the fit and, at the same time, to avoid negative speed values. This equation provides parameters with some physical meaning as the time approaches 0 or infinity. As t f þ ¥, the drop radius reaches some finite value (an/bn), and at the first point measured (t = 0) the radius has already some finite value (a0). To estimate the error on the speed from the radius data, we used the so-called bootstrap method. With this powerful method, the original data are used as the basis for a Monte Carlo simulation to build a new set of data from the original data. Some points are randomly chosen from the original data and are replaced by new data chosen according to a normal distribution function, with the original data as the average and the expected error on each datum point as the standard deviation. In this way, we replace the original set of data with a new one, for which the corresponding parameters are then fitted with the method described above. Successive cycles result in a simulated set of values for each parameter. If enough cycles are used (typically 1000), these sets turn out to be normally distributed, providing us with the mean value and standard deviation for the corresponding speed. The expected errors on each data point were fixed at 5
on angles yielding to radius error from 4 to 10%. This methodology has been implemented in the freeware G-Dyna33 which has been

Figure 7. Distribution of (K0, λ) and ln(L/Ls) (top right) parameters obtained from G-Dyna for d = 0.75. Langmuir 2010, 26(18), 14642–14647

DOI: 10.1021/la102412s

14645

Article

Seveno et al. Table 2. K 0, λ, ln(L/Ls), and θ0 Fitted by G-Dyna and K 0m, λm, and θ0m Measured Directly from the Simulationsa

d

K 0 (GHz)

1 1.87 (0.30) 0.75 0.86 (0.18) 0.7 0.62 (0.30) 0.6 0.46 (0.08) 0.5 0.44 (0.65) a Errors shown in parentheses.

K 0m (GHz)

λ (A˚)

λm (A˚)

ln(L/Ls)

θ0 (
)

θ0m (
)

1.03 (1.1) 0.79 (1.1) 0.83 (1.10) 0.37 (0.63) 0.20 (0.35)

3.39 (0.46) 5.23 (0.79) 5.35 (1.06) 6.06 (0.91) 5.97 (1.82)

4.75 (0.16) 5.22 (0.18) 5.34 (0.18) 5.63 (0.19) 5.98 (0.20)

81.44 (0.45) 59.08 (1.80) 84.47 (2.59) 80.77 (2.62) 111.50 (3.10)

41.65 (0.26) 72.18 (0.30) 76.21 (0.29) 84.10 (0.30) 92.30 (0.30)

46.59 (1.37) 73.23 (1.03) 77.95 (1.95) 86.00 (2.01) 94.30 (1.73)

designed to yield in particular the general MKT parameters θ0, λ, and K 0, together with their distribution and relative error. These data were used for the wetting analysis. In Figure 6, the measured TL velocities versus the dynamic contact angle are compared with the best MKT fit (the HD fits are of equal qualities) found with G-Dyna (the symbols represent every fifth data point). The distribution of the free parameters is shown in Figure 7 for d = 0.75 (the same type of graphs is obtained for the others densities). There is an obvious correlation between the parameters: a low K0 involves a high λ and vice versa, but the distributions are narrow enough so the calculation of an average value and standard deviation is meaningful. The corresponding MKT parameters, HD parameters, and equilibrium contact angles θ0 are listed in Table 2 using γ = 10.73 ( 2.7 mN/m for the surface tension.34 The error associated to every value takes into account the standard deviation of the parameters and the error on the surface tension. Let us first point out that the fits of the hydrodynamic model lead to a microscopic cutoff length which cannot be explained in terms of atomic length and is therefore not a model appropriate to describe the dynamics of our droplets. For comparison purposes, we also made direct measurements of atomic jump frequencies at the equilibrium, K 0m.30 In order to measure the frequencies, we count all the perpendicular jumps between the first and second layers. Since among these transitions thermal fluctuations do not have to be counted, we counted jumps only when a minimum perpendicular length was exceeded, fixed here to the adsorption site spacing. The average perpendicular jump frequencies are given by the ratio of the number of effective jumps divided by the time corresponding to the measurement; this ratio is then divided by the mean number of atoms within the first layer to obtain the average perpendicular jump frequency for an individual atom. The corresponding measurements are listed in Table 2.

IV. Results and Discussion The results reveal that θ0m is a monotonous function of the relative density d. We can compare these values to that predicted by the Cassie and Baxter relationship35,36 for the equilibrium contact angle between a liquid and a two-component surface: 0m cos θ0m ¼ R cos θ0m 1 þ ð1 - RÞcos θ2

ð6Þ

with R being the concentration of chemical species 1, and θ0m 1 and θ0m 2 the equilibrium contact angles on surfaces composed of pure species 1 and 2, respectively. In the present case, we can consider our surfaces to be composed of pure d = 1 surface with 0m θ0m 1 = 46.59
and vacuum with θ2 = 180
. When this is done, (33) Seveno., D.; Vaillant., A.; Rioboo., R.; Helena., H.; Conti, J.; De Coninck, J. Langmuir 2009, 25, 13034. (34) Seveno, D.; Ogonowski, G.; De Coninck, J. Langmuir 2004, 20, 8385. (35) Cassie, A. B. D.; Baxter, S. Discuss. Faraday Soc. 1948, 3, 11. (36) De Coninck, J.; Ruiz, J.; Miracle-Sol
e, S. Phys. Rev. E 2002, 65, 036139.

14646 DOI: 10.1021/la102412s

Figure 8. Cosine of contact angle with a 5
error bar versus d (black squares). The straight line corresponds to the Cassie-Baxter prediction.

the linear Cassie and Baxter relationship is recovered with R = d, and the resulting angle differs from our results by less than 5
; see Figure 8. We can see that there is a broad agreement between the measured and fitted frequencies. Due to the relative large width of the frequency distributions, the error bars (standard deviations) are large and then do not allow a detailed analysis. Yet, the trends are similar and therefore confirm that the MKT can be used to model our data. The values found can be interpreted as follows: when the surface density of solid atoms in contact with the liquid decreases, the liquid atoms feel a lower effective coupling. Previous studies30,31 have shown that for a given surface (with the distance between adsorption sites fixed) the jump frequency increases as the coupling decreases. In the present case, the distance between adsorption sites changes, and as expected the jump length λ is a decreasing function of d; that is, it increases with atom spacing. However, using the same liquid at the same temperature, we observe that increasing the distance between adsorption sites decreases the perpendicular jump frequency. To explain such a result, we must keep in mind that here the solid geometry is changing. Taking the same liquid in the same conditions, increasing the distance to jump must reduce its frequency. Indeed, the change in density but not in coupling will allow the attractive forces to become effectively more short-range, so that although the integral adhesion decreases, the short-range adhesion of the first layer increases. The liquid atoms will thus go closer to the solid surface as can be seen in Figure 2. As a result, the perpendicular jump frequency decreases. The agreement is also good between the jump length predicted by the MKT and λm except for d = 1 where λ is found to be too short compared to λm. Nevertheless, the comparison between fitted and measured values implies that taking the jump length as the adsorption site spacing makes sense. Langmuir 2010, 26(18), 14642–14647

Seveno et al.

Article

V. Conclusion In this Article, we have used molecular dynamics simulations involving a large number of atoms to model the wetting of a liquid drop on a solid surface. We have shown that the evolution of the drop may be modeled using the molecular-kinetic theory of wetting. Furthermore, we have demonstrated that the effective molecular parameters in the theory are consistent with those arising directly from the simulation, and, in particular, the jump distance closely follows the lattice spacing of the solid, which determines

Langmuir 2010, 26(18), 14642–14647

the adsorption sites in the simulations. These results provide strong support for the underlying validity of the molecular-kinetic theory of wetting. They also amply demonstrate the power of molecular dynamics simulations in relating the macroscopic dynamics to the microscopic properties of the system. Acknowledgment. This research is partly supported by the Fonds National pour la Recherche Scientifique, the European Community, and the R
egion Wallone.

DOI: 10.1021/la102412s

14647