Predicting the Wetting Dynamics of a Two-Liquid ... - ACS Publications

ARTICLE pubs.acs.org/Langmuir

Predicting the Wetting Dynamics of a Two-Liquid System D. Seveno,* T. D. Blake, S. Goossens, and J. De Coninck Laboratory of Surface and Interfacial Physics, Universite de Mons, 20 place du parc, 7000 Mons, Belgium ABSTRACT: We propose a new theoretical model of dynamic wetting for systems comprising two immiscible liquids, in which one liquid displaces another from the surface of a solid. Such systems are important in many industrial processes and the natural world. The new model is an extension of the molecular-kinetic theory of wetting and offers a way to predict the dynamics of a two-liquid system from the individual wetting dynamics of its parent liquids. We also present the results of large-scale molecular dynamics simulations for one- and two-liquid systems and show them to be in good agreement with the new model. Finally, we show that the new model is consistent with the limited data currently available from experiment.

’ INTRODUCTION When we talk about dynamic wetting, the image we have might be that of a raindrop splashing onto the windscreen of a car or rolling over the leaf of a plant. We might also think of the formation of a film of oil in the pan when we cook. In everyday life, wetting is mostly seen as the interaction between three phases: a solid, a liquid, and a gas (usually air), i.e., a one-liquid system. Many industrial applications also involve a single liquid wetting a solid, e.g., coating, mineral flotation, plant protection, paper printing, or soldering. Dynamic wetting of this type has been studied intensively for more than 50 years. Over this period, several models have been developed that attempt to explain the observed behavior. These include the molecular-kinetic theory (MKT),1,2 hydrodynamic models based on lubrication theory (HD),3
5 and the more recent hydrodynamic interface-formation model (IFM) of Shikhmurzaev.6,7 Other approaches based on the diffuse-interface model are also being investigated.8
10 The seminal review of de Gennes11 and the more recent ones by Blake12 and Bonn et al.13 show that the subject remains an active and challenging field of research, with new areas of interest constantly emerging, e.g., the wetting of structured or reactive surfaces, superhydrophobicity, and wetting by complex liquids. However, these reviews say little specifically about the wetting dynamics of two-liquid systems, i.e., those in which one liquid displaces another immiscible liquid. This does not mean that the subject is of little interest. Indeed, from a theoretical point of view, all the current models take into account the potential presence of a second liquid. These models have been tested against experimental data for liquid
liquid displacement in capillaries by Blake and Haynes1 for the MKT, by Mumley et al.14,15 and Fermigier and Jenffer16,17 for the HD model, and by Shikhmurzaev6 for the IFM. The agreement between experiment and prediction is good except for the HD model, with which the contact angles are generally larger than predicted.15,17,18 Foister18 and Fetzer et al.19 considered another geometry: the spontaneous displacement of a liquid during the spreading of an immiscible liquid drop. This is also the approach we have adopted. r 2011 American Chemical Society

Many natural and industrial processes are in fact dependent on the displacement of one liquid by another from the surface of a solid. Relevant topics include biosurfaces,20
23 oil recovery,24
29 groundwater cleanup,30 water
oil filtration,31
34 emulsification,35
39 microfluidics,39
42 fiber wetting,43 metal treatment,44 surface characterization,45,46 and detergency/cleaning processes.47
51 Though not exhaustive, this list highlights the interest and importance of two-liquid dynamic wetting. The addition of a second liquid, displaced by the first, not only involves an obvious change in the interfacial tensions, but also may enhance effects that depend, for example, on the viscosity ratio of the liquids. Furthermore, factors involving thin film rupture or the mutual miscibility of the liquids may come into play, especially in near-critical systems. Thus, it becomes clear that the problem is complex and requires careful study in its own right. For this reason, the strategy explored in this paper is first to investigate whether we can establish a link between one-liquid systems, which are now fairly well understood and easier to study, and two-liquid systems based on the same liquids and solids. It is generally accepted that the energy dissipation in the vicinity of the moving contact line controls the dynamics of wetting.11,52,53 Since several channels of dissipation are possible, it may be necessary to consider several different theoretical models to interpret the observed behavior. The MKT and HD model, or some combination thereof,52
54 are the most widely used. All these approaches attempt to explain the evolution of the dynamic contact angle θ as a function of the spreading velocity V, the material properties of the system, such as the liquid viscosity η and surface tension γ, and various model-specific parameters. Experimentally, it is known that advancing angles increase while receding angles decrease with increasing rates of contact-line displacement, and despite considerable efforts by many research groups, it appears that the data can often be fitted as well by one Received: September 6, 2011 Revised: October 28, 2011 Published: October 31, 2011 14958

dx.doi.org/10.1021/la2034998 | Langmuir 2011, 27, 14958–14967

Langmuir model as by another.12,55
57 In the present paper, we use the MKT to interpret our data. This is a first step, and other models may be considered in later publications. However, as it turns out, the MKT seems to be very effective in linking dynamic wetting in one- and two-liquid systems. A characteristic parameter that can be extracted from the MKT is the contact-line friction ζ.58 For a given system, this provides an index of the dissipation due to the movement of the three-phase contact line across the surface of the solid as the solid
liquid interface is either created or destroyed. With respect to the two-liquid case, this suggests an interesting question: is the contact-line friction for a two-liquid system a function of the contact-line frictions of the individual parent liquids spreading in air? Here, we present the results of numerical simulations designed to investigate this possibility. Various numerical simulation techniques have already been used to study the displacement of two-liquid systems confined between walls, including computational fluid dynamics (CFD)49
51,59 and lattice Boltzmann techniques.60,61 Molecular dynamics (MD) simulations have also been performed with this type of geometry to investigate slip between the liquid and the wall62
64 or the formation of a precursor film.65 MD has a number of advantages.66 The interactions between the liquids and with the solid surfaces are all defined by their relative affinities, which can be adjusted at will. With MD, unlike, for example, CFD, it is not necessary to prescribe any parameter such as the static contact angle or the slip length in the vicinity of the contact line. Such quantities appear directly as outputs from the simulation. Its nanoscopic length scale is also well adapted to the molecular description of the wetting process used in the MKT and so provides a direct way of evaluating the underlying model. Our approach, therefore, is to use large-scale MD simulations of droplet spreading for both one- and two-liquid systems and apply the MKT to investigate a potential link between the contact-line frictions obtained in each case. We also test our conclusions against experimental results obtained using liquids of different viscosities.67 That contact-line friction is a linear function of viscosity has recently been confirmed.68
70 It is, therefore, important to understand how this translates to the two-liquid case. Throughout our analysis, we make use of a systematic statistical tool, G-Dyna,57 to evaluate the dynamic contact angle data. Having such a method is valuable, as we can make the most of hard-won data to improve our understanding of the physical mechanism of wetting. The new work is an extension of the recent studies carried out in this laboratory by Bertrand et al.71 (MD simulations of one-liquid dynamic wetting) and Goossens et al.67 (experiments with one- and two-liquid systems). The paper is organized as follows. The MKT is first briefly reviewed, and a new model is proposed that links the contact-line friction of a two-liquid system to the frictions of its parent oneliquid systems. Then the simulation details are presented, followed by the results and analysis section, where the utility of the new model is discussed.

’ THEORETICAL MODELS When wetting equilibrium is disturbed, the contact angle deviates from its equilibrium value θ0. This leads to an out-ofbalance interfacial tension force F = γ(cos θ0
cos θ), which works to restore equilibrium. As demonstrated by de Gennes,11 this work balances three distinct channels of dissipation: viscous dissipation, dissipation in the close vicinity of the solid near the

ARTICLE

Figure 1. A liquid drop (L1) partially wetting a solid surface (S) showing a schematic illustration of the contact zone.

contact line, and in the complete wetting case (θ0 = 0) dissipation within a precursor film. For partial wetting (θ0 > 0), only the viscous losses and those associated with the contact line need be considered. According to this view, the viscous dissipation may be described by an HD model and dissipation near the contact line by the MKT. Voinov3 was the first to establish a hydrodynamic relation between the apparent dynamic contact angle θ measured at a mesoscopic height hm on the liquid surface, the local microscopic contact angle θm, and the contact-line velocity V. Voinov did not consider any phenomena occurring below hm, including any losses due to frictional processes between the liquid and the solid, though he specifically did not rule them out. Cox5 and others, including Dussan,4 assumed that the liquid slips on the solid in a region of length Ls (the slip length (m)) near the contact line and that θm is independent of V and can be set equal to the apparent equilibrium contact angle θ0. These theories have been widely applied in the literature with minor variations. More recently, Shikhmurzaev6 introduced a more comprehensive hydrodynamic model (the IFM, mentioned above) in which material transfer associated with the transformation of the liquid surface into a solid
liquid interface (and vice versa) enters directly into the equations of flow. In this model, the local microscopic angle varies with the flow and is equal to the observed dynamic angle. The first model to take account of the dissipation process occurring in the close vicinity of the contact line (and the first detailed model of dynamic wetting) was the molecular-kinetic theory of Blake and Haynes,1,2 which was based on the Frenkel
Eyring view of liquid transport as a stress-modified molecular rate process. According to their approach, the macroscopic behavior of the moving contact line is determined by the statistics of molecular displacements (jumps) occurring within the contact zone, i.e., the contact line when viewed on a molecular scale. This is illustrated schematically in Figure 1. If each jump has a length λ (m) and an equilibrium frequency k0 (Hz), then the velocity of the contact line is given by γðcos θ0
cos θÞ V ¼ 2k0 λ sinh 2nkB T

! ð1Þ

with kB the Boltzmann constant, T the temperature (K), and n the density of absorption sites on the solid surface. When data are being fitted, the latter is usually set equal to λ
2. This equation has been found to be effective in modeling the velocity 14959

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

dependence of the contact angle in a wide range of experimental systems2,12,55,67
69,72,73 and in MD simulations.66,71,74
77 When the argument of the sinh function is small, e.g., near equilibrium or when γ is small or n large, eq 1 can be simplified to γ ð2Þ V ¼ ðcos θ0
cos θÞ ζ with ζ (Pa 3 s) the friction per unit length of the contact line given by ζ ¼ kB T=k0 λ3

ð3Þ

Since k0 may be written in terms of the activation free energy of wetting (per mole) ΔGW* as   kB T
ΔGW  0 exp k ¼ ð4Þ h NA kB T we have

  h ΔGW  ζ ¼ 3 exp NA kB T λ

ð5Þ

where h and NA are the Planck constant and Avogadro number, respectively. If we consider ΔGW* to have contributions arising from both viscous interactions and interactions with the solid surface, ΔGV* and ΔGS*, respectively, such that ΔGW* = ΔGV* + ΔGS*, as suggested by Blake,2 then   ην ΔGS  ζ ¼ 3 exp ð6Þ NA k B T λ where ν is the molecular flow volume of the liquid, usually approximated by the molecular volume for simple liquids. Moreover, if we define ΔgS* = nΔGS*/NA as the specific activation free energy of wetting due to surface interactions (per unit area) and approximate it by the reversible work of adhesion between the liquid and solid Wa0 = γ(1 + cos θ0) (see Blake and De Coninck78), we finally obtain ! ην γð1 þ cos θ0 Þ ζ ≈ 3 exp ð7Þ nkB T λ Thus, ζ is proportional to the viscosity of the liquid and exponentially dependent on the work of adhesion. This relationship has been verified for one-liquid systems both experimentally68,69,72 and in simulations.71,79 The agreement is sufficient to suggest that the MKT, while comparatively simple and schematic, has underlying validity. Equation 1 seems to work equally well for one-liquid and twoliquid systems; indeed, it was first derived to explain the variation in the dynamic contact angle observed for benzene
water displacements in hydrophobic capillary tubes.1 An attempt was later made to extend the equation more explicitly to two-liquid systems2 by supposing that the viscous contributions to the activation free energy of wetting from each liquid were simply additive, i.e., ΔGW* = ΔGvis1* + ΔGvis2* + ΔGS*, where subscripts 1 and 2 distinguish the two liquids. This approach looked promising, but was never verified experimentally. Furthermore, in the two-liquid case it becomes very difficult to identify the surface contribution ΔGS* with any specific work of adhesion or combination thereof. As we show below, a fresh approach is needed. Our objective is to establish a relationship between key parameters such as k0 and ζ that can be measured for both one- and

Figure 2. Schematic illustration of the molecular process within the moving three-phase zone as a drop of L1 spreads on solid surface S immersed in a bath of L2.

two-liquid systems. To avoid ambiguity, as well as using subscripts 1 and 2 to specify liquid 1 or liquid 2 (L1 or L2) in contact with its vapor (or vapor and air), we use subscripts 1,2 and 2,1 to denote some property of L1 or L2 when in contact with the other liquid, L2 or L1. Thus, for example, θ01 refers to the equilibrium contact angle of L1 alone, while θ01,2 identifies the equilibrium contact angle measured through L1 in contact with L2. The interfacial tension between L1 and L2 is denoted γ12 in the usual manner. To proceed, we need to examine more closely the molecular mechanism by which the contact line advances or recedes in the two cases. In the one-liquid case, we can assume that the liquid
vapor interface is characterized more or less by the sigmoidal density distribution that we see in our MD simulations. 71 At the contact line, this merges with the more layered liquid structure found at the solid surface to form the contact zone. In the absence of a wetting film, the density of liquid will be low at the leading edge of this zone, so that molecular jumps in the forward direction will be comparatively unhindered and even backward jumps will find many vacant sites. In contrast, when a second liquid is present, instead of this low-density region, there will be a comparatively dense threephase contact zone (TPZ) in which the properties of one liquid will give way to those of the other over some short distance commensurate with the thickness of the L1
L2 interface. Therefore, unless there is significant depletion due, say, to very strong net repulsion between the two liquids, most of the sites within the TPZ will already be occupied by molecules of one liquid or the other. This will be enhanced by any higher density layering near the solid surface. Thus, a jump in the forward direction by a molecule of, say, L1 will require the prior creation of a vacancy into which it can jump. Consider a drop of L1 spreading on the surface of solid S immersed in a bath of immiscible L2 as depicted in Figure 2. At the leading edge of the drop we suppose L1 (black) is separated from L2 (white) by a TPZ comprising 50% L1 and 50% L2 molecules (gray). The detail shows the various steps required to translate the TPZ one adsorption site to the right according to the mechanism outlined below. The TPZ is arbitrarily shown as being about three molecules thick, but this does not affect our argument. At equilibrium, molecular displacements between L1 and L2 must occur to the left and right with the same overall 14960

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

frequency ( k012, with local fluctuations in the density and position of the TPZ. For unit molecular displacement of the TPZ from, say, left to right, a new molecule of L1 must arrive and a molecule of L2 must depart. The overall process of translating the TPZ by one molecular distance consists, therefore, of two steps. This differs from the one-liquid case, where sites to the right are unoccupied and we need to consider only the time required for an L1 molecule to move forward or backward. Intuitively, this two-step process suggests that the overall contact-line friction in the two-liquid case should be the sum of the respective frictions for each liquid on its own, i.e., ζ12 = ζ1 + ζ2. However, this neglects the fact that, for the TPZ to move, empty sites must be created within it as there are very few there initially. For unit displacement by one site from left to right, one vacancy must be created into which a molecule of L1 can move and another will be created when a molecule of L2 leaves. The energy required to empty these sites is in addition to that required to move a molecule in the one-liquid case (where there are already many vacancies) and must therefore be included in our formulation. With reference to Figure 2, let us consider the process step by step. For unit displacement from left to right, a molecule within the TPZ must leave the surface of the solid to create a vacancy (0) into which a molecule of L1 can move (L1 f 0), leaving a vacancy behind which will be filled from L1, as it would in the one-liquid case. Similarly, a molecule of L2 must leave the TPZ to fill a vacancy in L2 (L2 f 0), as it would also do in the one-liquid case, but leaving behind another vacancy (0). For a displacement from right to left, the sequence will simply be reversed. Thus, there will be a transient population of vacancies, without which the TPZ could not move. The various moves and the creation of the necessary vacancies will occur at different frequencies, but since the TPZ retains its statistical integrity, all moves must be completed for unit displacement. The chaotic motion of the liquid molecules means that the moves into and out of the TPZ will not usually occur within the same time interval or at precisely the same location. Hence, the two steps are unlikely to be a single sequential event. For each type of move, the time required will be given by (h/kBT) exp(ΔG*/NAkBT), where the prefactor (h/kBT) is the time required when the energy barriers are zero and ΔG* is the associated activation free energy. For L1 f 0 and L2 f 0, ΔG* will be about the same as for the individual liquids, i.e., ΔGW1* and ΔGW2*, respectively. However, for the creation of empty sites within the TPZ, since there are on average 50% of each species, we may expect that the activation free energy will be the average (ΔGW1* + ΔGW2*)/2. Effectively, what we are claiming here is that the work of adhesion of the liquid to the solid in the TPZ is the average of the works of adhesion of the two liquids independently. The total time required for all the moves to be completed will be the sum of the individual times τ12 = τL1f0 + τL2f0 + 2τ0. As a consequence, the time required for unit displacement of the TPZ by one adsorption site will be τ12 ¼

¼





h kB T

       ΔGW1  ΔGW2  ΔGW1  þ ΔGW2  exp þ exp þ 2 exp 2NA kB T N A kB T N A kB T

2 32  sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi h 4 ΔGW1  ΔGW2  5 exp þ exp N A kB T kB T N A kB T

ð8Þ

It follows, therefore, that 2sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    ffi
2   k T ΔG ΔG B W1 W2 4 exp 5 þ exp k012 ¼ h NA kB T NA kB T ð9Þ and since we can write the displacement frequency for each liquid on its own as     kB T
ΔGWi  k0i ¼ exp for i ¼ 1 or 2 ð10Þ h NA kB T we finally obtain 1 ¼ k012

1 1 pffiffiffiffiffi þ pffiffiffiffiffi 0 k1 k02

!2 ð11Þ

More significantly for our present purpose, since ζi µ 1/k0i , then provided λ1 ≈ λ2 pffiffiffiffiffi pffiffiffiffiffi2 pffiffiffiffiffiffiffiffiffi ζ12 ¼ ζ1 þ ζ2 or ζ12 ¼ ζ1 þ ζ2 þ 2 ζ1 ζ2 ð12Þ The last equation states that the contact-line friction for the twoliquid case is simply the sum of the frictions for each liquid alone on the same solid plus a geometric-mean mixing term that accounts for the additional friction due to the comparative lack of vacancies within the TPZ. The equation retains the intuitive appeal of the simple sum, but highlights the subtle difference between the one- and two-liquid cases. The constraint that λ1 ≈ λ2 could be significant, since λ appears in eq 3 to the third power. Nevertheless, provided λ is determined largely by the properties of the substrate, we can expect eq 12 to hold reasonably well. The form of eqs 11 and 12 is such that we recover the standard equations for the one-liquid case when either L1 or L2 is removed. However, it is possible that, near miscibility or for strong net repulsion between the liquids, the mixing term might be different. For example, if there is strong repulsion, then the mixing term might tend to zero, the two liquids acting independently, separated by a low-density gap. Conversely, near miscibility, the two liquids would begin to act as one and the mixing term would dominate. This suggests a generalization of eq 12: pffiffiffiffiffiffiffiffiffi ð13Þ ζ12 ¼ aðζ1 þ ζ2 Þ þ b ζ1 ζ2 where the constant a takes values from 0 (at miscibility) to 1 and the constant b increases from 0 (or even negative values) at pure repulsion to a maximum of 2 at normal immiscibility, before falling back to 1 at miscibility. In the following sections we test eq 12 against the results of large-scale molecular dynamics simulations and the experimental data of Goossens et al.67

’ SIMULATION MODEL Here, we summarize the main elements of the simulations: the geometry, the potentials, and the specific parameters, most of which have already been used to study the wetting of single liquids on flat surfaces,58,71 fibers,76,79 and porous substrates.75,80 The key parameters are the potentials between the solid (S) and 14961

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

Figure 3. Density profiles through the interfaces from two-liquid simulations for CL1
L2 = 0.8 (squares), 0.7 (circles), and 0.6 (inverted triangles) and from a one-liquid simulation (tilted squares).

liquid (L1 and L2) atoms, which are represented by standard pairwise Lennard-Jones 12
6 interactions: 2 !12 !6 3 σ σ ij ij 5 Uðrij Þ ¼ 4CA
B εij 4 ð14Þ
rij rij where rij is the distance between any pair of atoms i and j. The coupling parameter CA
B enables us to control the relative affinities between the atoms. The subscript A
B stands for the various possible interaction pairs: the diagonal interactions of the matrix L1
L1, L2
L2, and S
S and the off-diagonal interactions L1
L2, L1
S, and L2
S. 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 Å. The pair potential is set to zero for rij g 2.5σij (the cutoff length). CA
B is set to 1.0 for L1
L1, L2
L2, and S
S interactions. The L2
S coupling is chosen as 0.6 according to two opposing considerations. It has to be weak enough to allow for the drop to spread but strong enough to modify the spreading rate. The L1
L2 coupling is set to 0.6, 0.7, or 0.8. These levels ensure that the liquids are immiscible, with only slight depletion of liquid at the L1
L2 interface. Figure 3 shows the resulting density profiles together with that from the one-liquid simulations. For the weakest coupling, 0.6, the depletion is no more than 14% of the bulk density, which is consistent with our assumptions about the TPZ made in developing our model of two-liquid wetting. Changing the L1
L2 interaction changes the interfacial tension between the liquids. For a given CL1
L2, we vary CL1
S (from 0.2 to 1.0 in increments of 0.1). This modifies θ01,2 and therefore enables us to simulate a range of spreading behavior from nonspreading to complete wetting. To limit evaporation and increase viscosity, the liquid molecules are constructed as eight-atom chains with a confining potential between nearest neighbors i and j, Uconf(rij) = Arij6. The constant A is set to εij/σij6. The solid substrate is a square planar lattice comprising five atomic layers. The distance between the solid atoms is set to 21/6σij (3.93 Å), i.e., the equilibrium distance given by the Lennard-Jones potential. The atoms can vibrate around their initial equilibrium position according to a harmonic 0 2 potential defined by Uh(r Bi) = D|r Bi
Br0 i | with Br i the instantaneous position of a solid atom i and Br i its initial position. The constant D is given the value 2.5(εijσij2).

ARTICLE

The mass of all three types of atoms is set to 12 g/mol (the molar mass of carbon) and the time step to 5 fs. A reduced temperature of T* = 1.0 is used, i.e., T* = kBT/εij with T = 33.33 K. During spreading, the temperature of only the solid is kept constant via a velocity scaling algorithm. This allows us to mimic a real isothermal solid. Although our model is quite simplistic, it contains all the basic ingredients required to investigate the details of the spreading process. To analyze our data, we first need to compute the interfacial tensions between L1 and L2 and between L1 and its vapor. To do this, we construct the interfaces with planar liquid films and evaluate the interfacial tensions using81 * ! + 3rijz 2 ∂Uðrij Þ 1 N N γ¼ 1
2 rij ð15Þ 2S i ¼ 1 j > i rij ∂rij

∑∑

where rijz is the z component of the separation vectorBr i
Br j and S the area of the interface. The sum over all the liquid atoms is such that only the atoms at the interface contribute. For L1, we obtained γ1 = 2.49 mN/m, in agreement with Bertrand et al.,71 and for the two-liquid simulations, we found γ12 = 2.04, 2.93, and 3.75 mN/m for CL1
L2 = 0.8, 0.7, and 0.6, respectively, with uncertainties that average (0.5 mN/m. Thus, as expected, the interfacial tension increases as the L1
L2 coupling is reduced. These values can be correlated with the thickness of the interface.71 For CL1
L2 = 0.8, the interface is thicker than the liquid
vapor interface (9.68 ( 0.03 Å compared with 8.67 ( 0.02 Å). However, for CL1
L2 = 0.7 and 0.6, the thickness is smaller, equal to 8.00 ( 0.02 and 6.74 ( 0.02 Å, respectively. This means that the interface becomes thinner as the coupling is reduced and the tension increases. For the spreading experiments, the system is constructed as follows: First, an L2 slab (308 Å  308 Å  177 Å) is formed using 35 874 eight-atom molecules and equilibrated against the solid surface. Periodic boundary conditions (PBCs) are used in the x and y directions to extend the effective size of the simulation. Liquid atoms are then selected to create a spherical L1 drop (3449 molecules) in the center of the film. This two-liquid system is then equilibrated at a very low L1
S coupling (0.2) to ensure that the drop does not spread. The L1
S coupling is then increased to the chosen value to initiate spreading. For the oneliquid simulations, the procedure is exactly the same, except that the surrounding liquid is removed after the drop has been created (Figure 4). The same procedure and parameters were used by Bertrand et al.,71 but the drops were larger (5000
27040 molecules for the lowest L1
S couplings). Spreading dynamics are studied in two stages. First, one-liquid simulations are run in a vacuum; then the simulations are run in the presence of the second liquid. So that we can compare the results, the same set of parameters is used in each case.

’ RESULTS AND DISCUSSION To describe the dynamics of wetting, we need to determine the evolution of the contact angle and contact-line position with time. At each stage, we locate the edge of the drop by a density calculation71 and record its position. We then approximate the drop shape by a spherical cap and fit it to a circle. At the beginning of the simulations, when the drop first comes into contact with the solid, the shape is nonspherical, so these configurations are omitted. Depending on the liquid
solid coupling, we have to wait 0.5
1.5 ns for the one-liquid simulations and 2
3 ns for the 14962

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

Figure 4. Equilibrium drop profiles from a two-liquid simulation (left, θ01,2 = 50.7°) and its parent one-liquid simulations (middle, θ01 = 74.7°; right, θ02 = 107.0°) with CL1
S = 0.8, CL2
S = 0.6, and CL2
L1 = 0.8. The dynamics of the contact angle leading to these profiles are shown in Figure 6.

Table 1. Static Contact Angles θ01, θ02, and θ012 a θ01,2 (deg) coupling CL1
S

a

θ01 (deg)

θ01 (deg), Bertrand et al.71

CL1
L2 = 0.8

CL1
L2 = 0.7

CL1
L2 = 0.6

0.4

130.6 ( 2.4

134.8 ( 1.1

123.2 ( 1.9

109.9 ( 1.0

108.4 ( 1.1

0.5

118.1 ( 1.8

121.3 ( 0.8

107.1 ( 0.8

100.5 ( 1.3

100.5 ( 1.5

0.6

107.0 ( 2.5

105.8 ( 2.4

92.2 ( 1.0

89.6 ( 1.4

89.6 ( 1.2

0.7

91.1 ( 1.8

91.2 ( 1.9

72.4 ( 0.9

77.7 ( 1.8

78.0 ( 0.9

0.8

74.7 ( 1.6

74.6 ( 1.7

50.7 ( 1.6

62.4 ( 1.0

69.5 ( 1.2

0.9 1.0

55.7 ( 1.5 36.5 ( 0.9

55.4 ( 1.6 31.8 ( 0.8

* *

44.2 ( 0.6 25.4 ( 0.7

54.6 ( 0.7 38.9 ( 0.7

An asterisk indicates that the contact line reached the boundaries of the simulation box before equilibrium was attained.

two-liquid simulations before the spherical approximation is valid. In addition, the profile of the first few layers of liquid at the L1
S interface is perturbed for energetic and entropic reasons.82 To avoid this latter problem, we divided the drop into horizontal layers and investigated the profile as a function of the number of layers used, from top to bottom. Usually, no more than the first two layers need to be omitted from the fit. Further details can be found in the paper by Bertrand at al.71 The tangent to the circular profile at the intersection with the substrate defines the contact angle θ and the base radius R. Using this procedure, we can measure both parameters as a function of the number of time steps during our simulations and also θ0 when the drop reaches equilibrium. Static Contact Angle. Bartell and Osterhof83 have established a simple thermodynamic relation between the static contact angles θ01, θ02, and θ01,2. This predicts that ð16Þ

Figure 5. Comparison between the static contact angles measured in the simulations and those predicted by the Bartell
Osterhoff relation for CL1
L2 = 0.8 (squares), 0.7 (circles), and 0.6 (inverted triangles). The line represents perfect agreement.

and is derived from Young’s equation on the basis of negligible adsorption of L1 at the L2
S interface or L2 at the L1
S interface. The relation has been confirmed experimentally, although problems with contact angle hysteresis can make this difficult.84 In our simulations, the solid surface is both physically and chemically perfect and there is no hysteresis. Therefore, eq 16 should provide a useful way to test the quality of our data and therefore validate our methods. Table 1 lists equilibrium contact angles obtained from our simulations together with those reported by Bertrand et al.71 The agreement between the values of θ01 from both studies is very good and usually well within the standard deviations. The results for θ01,2 are compared with the Bartell
Osterhof prediction in Figure 5. The agreement here is excellent for contact angles below 100°, with a slight discrepancy at higher

angles that is attributed to the relatively high uncertainty on θ01 at large angles. Another way to check our simulations is to measure θ01,2 for a given CL1
L2, use eq 16 to predict γ12, and compare that with the value obtained independently using eq 15. For CL1
L2 = 0.4 and CL1
S = 0.8, θ01,2 = 72.6° ( 1.0°, which leads to γ12 = 4.74 ( 0.20 mN/m. The calculated value is 5.10 ( 0.50 mN/m. The agreement is within the error range and provides further evidence of the reliability of our methods. Dynamic Contact Angle and Contact-Line Friction. Figure 6 shows (left) the evolution of the dynamic contact angle with time for two-liquid simulations with (CL1
L2 = 0.8) and (right) a comparison with typical dynamics obtained for single liquids. For clarity, the error bars (∼2°) are not shown. It is clear that the time required to reach equilibrium is longer in the two-liquid case.

cos θ01, 2 ¼

γ1 cos θ01
γ2 cos θ02 γ12

14963

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

Figure 6. Left: evolution of the dynamic contact angle for two-liquid simulations with CL1
L2 = 0.8, CL2
S = 0.6, and CL1
S = 0.4
1.0 (from top to bottom). Right: comparison between the dynamic contact angles obtained for one-liquid simulations with CL1
S = 0.6 (top) and 0.8 (middle) and the two-liquid simulation with CL2
S = 0.6 and CL1
S = 0.8 (bottom).

Figure 7. θ1,2 (left) and cos θ1,2 (right) versus V for two-liquid simulations with CL1
L2 = 0.8 and CL1
S = 0.4
0.8 (squares). The black lines are the best fits to eq 2 obtained with G-Dyna.

This is due to two factors: the equilibrium contact angle is smaller (hence, the drop has further to spread), and the presence of the second liquid slows the dynamics. Data fitting to obtain the contact-line frictions is carried out with the subroutines implemented in the G-Dyna software.57 To determine the speed of the contact line, we use a Levenberg
Marquardt85 procedure to fit the dynamic base radius R(t) to a ratio of polynomials on the order of e10, which accurately reproduces the observed evolution of the radius. The Levenberg
Marquardt procedure was used because it is a robust curve-fitting algorithm that is able to find a solution even if the initial guess of the free parameters is far from the final solution. By simple differentiation of the fitted R(t), the speed of the contact line V is calculated and associated with the value of the corresponding contact angle. This is then fitted to the linear version of the MKT to give ζ1 or ζ12, as appropriate. The bootstrap method is used to estimate the error on the free parameters. Figure 7 (left) shows the contact angle data from the two-liquid simulations plotted against the contact-line speed together with the fitted curves for coupling parameters ranging from 0.4 to 0.8. The relation between

the cosine of the contact angle and the velocity (Figure 7, right) demonstrates that the data lie in the linear range of the MKT and confirms that the friction is the relevant parameter to describe the dynamics. The same procedure is used to obtain ζ1 and ζ2 from the oneliquid simulations (ζ2 is obtained from the simulation with CL1
S = 0.6). The spreading velocities with two liquids are, on average, 4 times slower than in the one-liquid case. The corresponding contact-line frictions are therefore expected to be greater. Table 2 lists the values ζ1 and ζ12 obtained in the present work together with those for ζ1 reported by Bertrand et al.71 The values of ζ1 from the two sources are compatible and agree within the uncertainty limits. This indicates that the results are independent of the size of the simulations and that the methods are consistent. Therefore, to improve the statistics for subsequent calculations, we take the average values. Figure 8 (left) compares the contact-line frictions for the different coupling parameters used. The data are plotted against the L1
S coupling. As expected, ζ12 is systematically larger than ζ1. Furthermore, the general trend is the same in both cases, with 14964

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

Table 2. Contact-Line Frictions (mPa 3 s) Obtained for Each Couplinga ζ1,2 (mPa 3 s) coupling CL1
S

ζ1 (mPa 3 s)

ζ1 (mPa 3 s), Bertrand et al.71

CL1
L2 = 0.8

CL1
L2 = 0.7

CL1
L2 = 0.6

average

0.4

**

0.16 ( 0.04

1.03 ( 0.25

0.92 ( 0.17

0.88 ( 0.12

0.96 ( 0.19

0.5

0.17 ( 0.03

0.21 ( 0.06

1.00 ( 0.25

1.04 ( 0.18

1.12 ( 0.15

1.05 ( 0.19

0.6

0.26 ( 0.05

0.36 ( 0.10

1.42 ( 0.35

1.03 ( 0.18

1.18 ( 0.16

1.21 ( 0.23

0.7

0.39 ( 0.08

0.43 ( 0.11

1.54 ( 0.38

1.24 ( 0.21

1.53 ( 0.20

1.44 ( 0.26

0.8

0.56 ( 0.11

0.53 ( 0.14

1.66 ( 0.41

1.65 ( 0.28

1.56 ( 0.21

1.62 ( 0.30

0.9

0.81 ( 0.15

0.75 ( 0.20

*

1.63 ( 0.28

1.54 ( 0.21

1.59 ( 0.25

1.0

1.06 ( 0.21

0.87 ( 0.23

*

1.72 ( 0.29

1.59 ( 0.21

1.66 ( 0.25

a

A single asterisk indicates that the contact line reached the boundaries of the simulation box before equilibrium was attained. Two asterisks indicate that the dynamics were too noisy to obtain a reliable value for ζ1.

Figure 8. Left: contact-line friction versus L1
S coupling for the one-liquid simulation (up triangles) and two-liquid simulation for CL1
L2 = 0.8 (squares), CL1
L2 = 0.7 (circles), and CL1
L2 = 0.6 (down triangles). Right: same data with ζ1,2 averaged (filled squares). The full lines are linear fits (for ζ12, the two data points at the highest L1
S couplings are omitted from the fits; see the text).

friction increasing with coupling. This is fully consistent with the theoretical predictions of the MKT and, increasingly, with experiment. For the two-liquid simulations, there is no clear trend with L1
L2 coupling over the range investigated. This allows us to average the values of ζ12 at each CL1
S. These averages are plotted in Figure 8 (right). As can be seen, the difference between ζ12 and ζ1 increases only slowly with coupling. Since ζ2 is constant (average value 0.31 mPa 3 s) for all the two-liquid simulations, this result is consistent with eq 12. Note that in this case there is no issue with possible differences between λ1 and λ2. Bertrand et al.71 have shown that λ is independent of L
S coupling in these simulations and has a value of 4.3 ( 0.4 Å, which is slightly larger than the solid lattice spacing (3.93 Å). It is perhaps surprising that the L1
L2 coupling has no detectable effect on ζ12, as it does have a significant influence on the thickness of the L1
L2 interface, which one would expect to be reflected in the width of the TPZ. However, within the accuracy of our simulations, this does not seem to affect the dynamics of wetting. As a definitive check on the utility of our new model, in Figure 9 we have plotted the average values of ζ12, those predicted by eq 12, and the result of simply adding ζ1 and ζ2. As we can see, the agreement with eq 12 is very good and much better than with the simple sum at all L1
S couplings for which we have results except 0.9 and 1.0. From these two data points, it appears

that ζ12 has ceased to increase, even though ζ1 continues its usual upward trend. At present, we cannot account fully for this apparent anomaly. Neither the equilibrium contact angles nor the simulations suggest the existence of a precursor film of L1, which would effectively fix the apparent friction. The effect could be an artifact of the simulations or our analysis of the results. However, on the basis of the evidence so far, it does seem that, at high L1
S couplings, the L1 friction begins to dominate and the extra resistance provided by L2 becomes less important. One possibility is that, as CL1
S increases, τL1f0 eventually becomes sufficiently long compared with τL2f0 for the displacement process in the TPZ to become one sequential event. Only one vacancy is then required, which the departing L2 molecule provides. Hence, eq 12 reduces to ζ12 = ζ1 + ζ2 + (ζ1ζ2)1/2. As a result, ζ12 should fall between the outcome predicted by our original eq 12 and the simple sum rule, which is what we observe. Effectively, the L1 f 0 move becomes the principle ratedetermining step. We will investigate this possibility in future research. Overall, however, the very good agreement at all lower L1
S couplings suggests that eqs 11 and 12 and the model on which they are based represent significant progress toward a better understanding of two-liquid dynamic wetting. Supporting Experimental Evidence. Although there are several notable studies of two-liquid dynamic wetting,14
19 there 14965

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir

ARTICLE

given in Table 3 and plotted in Figure 10. The value of ζ2 = 0.018 Pa 3 s was obtained from a separate analysis of their data for water. While there are only three two-liquid systems, the experiments cover a fairly wide range of viscosity (1
33 mPa 3 s). As in Figure 9, the predictions are close to the measurements, although in this case there seems to be a slight but consistent underestimate, even at low viscosity. For the simple sum ζ1 + ζ2 the agreement is worse. Overall, the results are very encouraging, but there is clearly room for more experiments of this type. Indeed, it has recently come to our attention that another group has been working independently on establishing a link between one- and two-liquid dynamic wetting and will shortly publish new experimental data that should provide further insight.86 Figure 9. ζ12 obtained from the simulations (black squares), predicted by eq 12 (open squares), and obtained by simply adding ζ1 and ζ2 (circles) versus CL1
S.

Table 3. Viscosity and Contact-Line Frictions ζ12 (Pa 3 s) for Various Alkane Drops Spreading on a Hydrophobic Surface in Air and Water (from Goossens et al.67)a ζ12 (Pa 3 s)

a

η (Pa 3 s)

ζ1 (Pa 3 s)

measured

predicted

dodecane

0.001

0.511 ( 0.004

1.395 ( 0.008

0.721 ( 0.008

hexadecane

0.003

0.838 ( 0.04

1.585 ( 0.04

1.102 ( 0.04

squalane

0.033

4.65 ( 0.005

5.87 ( 0.002

5.247 ( 0.002

The contact-line friction for water was ζ2 = 0.018 Pa 3 s.

’ CONCLUSION Large-scale molecular dynamics has been used to investigate the spreading of liquid drops on a solid surface, both in the presence and in the absence of a second liquid. The attraction between the spreading liquid and the solid has been varied systematically to simulate a wide range of wetting conditions. Good agreement has been found between the resulting equilibrium contact angles and the Bartell
Osterhof relation, which tends to validate our general approach. More significantly, we have shown that the spreading dynamics can be modeled by an extension of the molecular-kinetic theory using a revised model of the TPZ. In particular, we have derived an equation linking the contact-line friction of a two-liquid system to the sum of the frictions of the individual liquids on the same substrate plus a geometric mean mixing term that accounts for the additional friction of the TPZ. Comparison of this prediction with the results of our simulations shows encouraging agreement. The limited experimental data currently available lend further support. We believe that this is a valuable first step toward a better understanding of two-liquid spreading dynamics. However, there is a particular need for more experimental work. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

Figure 10. Contact-line frictions ζ12 (Pa 3 s) from Goossens et al.67 (black squares), predicted by eq 12 (open squares), and obtained by simply adding ζ1 and ζ2 (circles) versus the viscosity of the alkanes. From left to right, dodecane, hexadecane, and squalane drops spreading under water. The full line is a guide to the eye.

are very few experimental data that can be used to test eqs 11 and 12 directly. To the best of our knowledge, the only published study that attempts to establish a relation between one-liquid and two-liquid wetting using the same liquids and the same solid for all the experiments is that reported by Goossens et al.67 They measured the dynamic contact angle during the spreading of dodecane, hexadecane, and squalane drops on a hydrophobic grafted silicon substrate both in air and immersed in water. The contactline frictions were extracted from the dynamics of spreading using essentially the same procedures as described here. The results are

’ ACKNOWLEDGMENT This research is partly supported by the Fonds National pour la Recherche Scientifique, the European Community, and the Region Wallonne. ’ REFERENCES (1) Blake, T. D.; Haynes, J. J. Colloid Interface Sci. 1969, 30, 421. (2) Blake, T. D. In Wettability; Berg, J. C., Ed.; Marcel Dekker: New York, 1993; pp 251
309. (3) Voinov, O. V. Fluid Dyn. 1976, 11, 714–721. (4) Dussan, V. E. B. J. Fluid Mech. 1976, 77, 665–684. (5) Cox, R. G. J. Fluid Mech. 1986, 16, 169–194. (6) Shikhmurzaev, Y. J. Multiphase Flow 1993, 19, 589–610. (7) Shikhmurzaev, Y. J. Fluid Mech. 1997, 334, 211–249. (8) Seppecher, P. Int. J. Eng. Sci. 1996, 34, 977. (9) Jacqmin, D. J. Fluid Mech. 2000, 402, 57. (10) Yue, P.; Zhou, C.; Fend, J. J. Fluid Mech. 2010, 645, 279–294. (11) De Gennes, P. G. Rev. Mod. Phys. 1985, 57, 827. (12) Blake, T. D. J. Colloid Interface Sci. 2006, 299, 1–13. 14966

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967

Langmuir (13) Bonn, D.; Eggers, J.; Indekeu, J.; Meunier, J.; Rolley, E. Rev. Mod. Phys. 2009, 81, 739–804. (14) Mumley, T.; Radke, C. J.; Williams, M. J. Colloid Interface Sci. 1986, 109, 398–412. (15) Mumley, T.; Radke, C. J.; Williams, M. J. Colloid Interface Sci. 1986, 109, 413–425. (16) Fermigier, M.; Jenffer, P. Ann. Phys. 1988, 13, 37–42. (17) Fermigier, M.; Jenffer, P. J. Colloid Interface Sci. 1991, 146, 226. (18) Foister, R. J. Colloid Interface Sci. 1990, 136, 266–282. (19) Fetzer, R.; Ramiasa, M.; Ralston, J. Langmuir 2009, 25, 8069–8074. (20) Frish, T.; Thoumine, O. J. Biomech. 2002, 35, 1137–1141. (21) Dobereiner, H.; Dubin-Thaler, B.; Giannone, G.; Xenias, H.; Sheetz, M. Phys. Rev. Lett. 2004, 93, 108105. (22) Marsh, R.; Jones, R.; Sferrazza, M. Colloids Surf., B 2002, 23, 31–42. (23) Velzenberger, E.; Kirat, K. E.; Legeay, G.; Nagel, M.; Pezron, I. Colloids Surf., B 2009, 68, 238–244. (24) Asserson, R. B.; Hoffmann, A. C.; Høiland, S.; Asvik, K. M. J. Pet. Sci. Eng. 2009, 68, 209–217. (25) Freer, E. M.; Svitova, T.; Radke, C. J. J. Pet. Sci. Eng. 2003, 39, 137–158. (26) Guy, D. W.; Crawford, R. J.; Mainwaring, D. E. Fuel 1996, 75, 238–242. (27) dos Santos, R. G.; Mohamed, R. S.; Bannwart, A. C.; Loh, W. J. Pet. Sci. Eng. 2006, 51, 9–16. (28) Basu, S.; Nandakumar, K.; Masliyah, J. H. Colloids Surf., A 1998, 136, 71–80. (29) Basu, S.; Nandakumar, K.; Lawrence, S.; Masliyah, J. Fuel 2004, 83, 17–22. (30) Lim, T.-T.; Huang, X. J. Hazard. Mater. 2006, 137, 820–826. (31) Lee, C. H.; Johnson, N.; Drelich, J.; Yap, Y. K. Carbon 2011, 49, 669–676. (32) Kocherginsky, N. M.; Tan, C. L.; Lu, W. F. J. Membr. Sci. 2003, 220, 117–128. (33) Kong, J.; Li., K. Sep. Purif. Technol. 1999, 16, 83–93. (34) Bansal, S.; von Arnim, V.; Stegmaier, T.; Planck, H. J. Hazard. Mater. 2011, 190, 45–50. (35) Kumar, D.; Biswas, S. K. Colloids Surf., A 2011, 377, 195–204. (36) Le Follotec, A.; Pezron, I.; Noik, C.; Dalmazzone, C.; MetlasKomunjer, L. Colloids Surf., A 2010, 365, 162–170. (37) Sonin, A.; Palermo, T.; Lubek, A. Mater. Chem. Phys. 1998, 56, 74–78. (38) Cambiella, A.; Benito, J. M.; Pazos, C.; Coca, J. Colloids Surf., A 2007, 305, 112–119. (39) Kobayashi, I.; Nakajima, M.; Mukataka, S. Colloids Surf., A 2003, 229, 33–41. (40) Spildo, K.; Buckley, J. S. J. Pet. Sci. Eng. 1999, 24, 145–154. (41) Tostado, C. P.; Xu, J.; Luo, G. Chem. Eng. J. 2011, 171, 1340–1347. (42) Bashir, S.; Rees, J. M.; Zimmerman, W. B. Chem. Eng. Sci. 2011, 66, 4733–4741. (43) Perwuelz, A.; Olivera, T. N. D.; Caze, C. Colloids Surf., A 1999, 147, 317–329. (44) Kumar, D.; Biswas, S. K. Colloids Surf., A 2010, 356, 112–119. (45) Jung, Y.; Bhushan, B. Langmuir 2009, 25, 14165–14173. (46) Rosa, M.; de Pinho, M. J. Membr. Sci. 1997, 131, 167–180. (47) Seevaratnam, G.; Ding, H.; Michel, O.; Heng, J.; Matar, O. Chem. Eng. Sci. 2010, 65, 4523–4534. (48) Chatterjee, J. Colloids Surf., A 2001, 178, 249. (49) Thoreau, V.; Malki, B.; Berthome, G.; Boulange-Petermann, L.; Joud, J. J. Adhes. Sci. Technol. 2006, 20, 1819–1831. (50) Ding, H.; Spelt, P. J. Fluid Mech. 2008, 599, 341–362. (51) Ding, H.; Gilani, M.; Spelt, P. J. Fluid Mech. 2010, 644, 217–244. (52) Brochard-Wyart, F.; Gennes, P. G. D. Adv. Colloid Interface Sci. 1992, 39, 1–11. (53) De Ruijter, M. J.; De Coninck, J.; Oshanin, G. Langmuir 1999, 15, 2209–2216.

ARTICLE

(54) Petrov, P. G.; Petrov, J. Langmuir 1992, 8, 1762–1767. (55) Schneemilch, M.; Hayes, R. A.; Petrov, J. G.; Ralston, J. Langmuir 1998, 14, 7047–7051. (56) Ranabothu, S. R.; Kaenezis, S.; Dai, L. L. J. Colloid Interface Sci. 2005, 288, 213–221. (57) Seveno, D.; Vaillant, A.; Rioboo, R.; Adao, H.; Conti, J.; De Coninck, J. Langmuir 2009, 25 (22), 13034–13044. (58) De Ruijter, M. J.; Blake, T. D.; De Coninck, J. Langmuir 1999, 15, 7836–7847. (59) Zhou, G.; Chen, Z.; Ge, W.; Li, J. Chem. Eng. Sci. 2010, 65, 3363–3371. (60) Blake, T.; De Coninck, J.; D’Ortona, U. Langmuir 1995, 11, 4588–4592. (61) Kang, Q.; Zhang, D.; Chen, S. Adv. Water Resour. 2004, 27, 13–22. (62) Koplik, J.; Banavar, J.; Willemsen, J. Phys. Rev. Lett. 1988, 60, 1282–1285. (63) Thompson, P. A.; Robbins, M. O. Phys. Rev. Lett. 1989, 63, 766–769. (64) Qian, T.; Wang, X.-P.; Sheng, P. J. Fluid Mech. 2006, 564, 333–360. (65) Wu, C.; Qian, T.; Sheng, P. J. Phys.: Condens. Matter 2010, 22, 325101. (66) De Coninck, J.; Blake, T. Annu. Rev. Mater. Res. 2008, 38, 15.1–15.22. (67) Goossens, S.; Seveno, D.; Rioboo, R.; Vaillant, A.; Conti, J.; De Coninck, J. Langmuir 2011, 27, 9866–9872. (68) Vega, M. V.; Gouttiere, C.; Seveno, D.; Blake, T. D.; Voue, M.; De Coninck, J. Langmuir 2007, 23, 10628–10634. (69) Li, H.; Sedev, R.; Ralston, J. Phys. Chem. Chem. Phys. 2011, 13, 3952–9. (70) Duvivier, D.; Blake, T.; Seveno, D.; Rioboo, R.; De Coninck, J. Langmuir 2011, 27, 13015–13021. (71) Bertrand, E.; Blake, T.; De Coninck, J. J. Phys.: Condens. Matter 2009, 21, 464124. (72) Puah, L.; Sedev, R.; Fornasiero, D.; Ralston, J.; Blake, T. Langmuir 2010, 26, 17218–17224. (73) Saiz, E.; Tomsia, A. Curr. Opin. Solid State Mater. Sci. 2005, 9, 167–173. (74) Heine, D.; Grest, G.; Webb, E. Phys. Rev. E 2003, 68, 061603. (75) Martic, G.; De Coninck, J.; Blake, T. D. J. Colloid Interface Sci. 2003, 263, 213–216. (76) Seveno, D.; Ogonowski, G.; De Coninck, J. Langmuir 2004, 20, 8385–8390. (77) Benhassine, M.; Saiz, E.; Tomsia, A.; Coninck, J. D. Langmuir 2009, 25, 11450–11458. (78) Blake, T. D.; De Coninck, J. Adv. Colloid Interface Sci. 2002, 96, 21–36. (79) Seveno, D.; De Coninck, J. Langmuir 2004, 20, 737–742. (80) Martic, G.; Gentner, F.; Seveno, D.; De Coninck, J.; Blake, T. D. J. Colloid Interface Sci. 2004, 270, 171–179. (81) Salomons, E.; Mareschal, M. J. Phys.: Condens. Matter 1991, 3, 3645. (82) De Coninck, J.; Dunlop, F.; Menu, F. Phys. Rev. E 1993, 47, 1820. (83) Bartell, F.; Osterhof, H. Ind. Eng. Chem. 1927, 19 (11), 1277–1280. (84) Johnson, R. E. In Wettability; Berg, J. C., Ed.; Marcel Dekker: New York, 1993; pp 1
74. (85) Marquardt, D. W. SIAM J. Appl. Math. 1963, 11, 431–441. (86) Ramiasa, M.; Ralston, J.; Fetzer, R.; Sedev, R. J. Phys. Chem. C 2011; Just Accepted.

14967

dx.doi.org/10.1021/la2034998 |Langmuir 2011, 27, 14958–14967