Dissipative Particle Dynamics Simulation of Contact Angle Hysteresis

Jan 26, 2006 - We have studied two types of topological substratessthe continuous solid substrates (CSS) and the discontinuous solid substrates (DSS)s...
0 downloads 0 Views 630KB Size
Download PDF
Langmuir 2006, 22, 2065-2073

2065

Dissipative Particle Dynamics Simulation of Contact Angle Hysteresis on a Patterned Solid/Air Composite Surface Bin Kong and Xiaozhen Yang* State Key Laboratory of Polymer Physics and Chemistry, Joint Laboratory of Science and Materials, Institute of Chemistry, Chinese Academy of Sciences, Beijing 100080, China ReceiVed July 22, 2005. In Final Form: December 27, 2005 We have studied two types of topological substratessthe continuous solid substrates (CSS) and the discontinuous solid substrates (DSS)sby using the dissipative particle dynamics (DPD) method for a better understanding of the contact angle hysteresis on two such substrates. After the validation of DPD in the system, we found that DSS has a different distribution of the metastable states from that of CSS and that DSS has relatively larger contact angle hysteresis at lower temperature. Obtained results also show that CSS is more suitable for making an ultrahydrophobic or ultralyophobic surface than DSS from the point of view of dynamic wettability.

1. Introduction Contact angle hysteresis is an important phenomenon in the wetting process and has been studied widely.1-19 This is because, for ideal surfaces, a unique equilibrium contact angle (θe) is determined theoretically by thermodynamic equations.9,20-27 However, the contact angle on practical surfaces is seldom singlevalued. The observed contact angle of a liquid advancing across a dry solid surface is always larger than that as it recedes. The former is called the “advancing contact angle” (θa), and the latter the “receding contact angle” (θr). Traditionally, the difference in advancing and receding contact angles has been denoted as contact angle hysteresis (θhyst) (i.e., θhyst ) θa - θr). The magnitude of θhyst is more important than that of θe (or θa) in evaluating the dynamic wettability of a real surface2,8,14,16 because θhyst represents the motion ability of a droplet on an inclined plate. Understanding contact angle hysteresis is essential for designing ultrahydrophobic surface, and a low θhyst is preferred because liquid droplets more easily slide on such inclined surfaces. * To whom all correspondence should be addressed. E-mail: [email protected]. (1) Bico, J.; Marzolin, C.; Quere, D. Europhys. Lett. 1999, 47, 220-226. (2) Chen, W.; Fadeev, A. Y.; Hsieh, M. C.; Oner, D.; Youngblood, J.; McCarthy, T. J. Langmuir 1999, 15, 3395-3399. (3) Decker, E. L.; Garoff, S.Langmuir 1996, 12, 2100-2110. (4) Derjaguin, B. V. C. R. Acad. Sci. USSR. 1946, 51, 361. (5) Drelich, J. Pol. J. Chem. 1997, 71, 525-549. (6) Dussan, E. B. Annu. ReV. Fluid Mech. 1979, 11, 371-400. (7) Eick, J. D.; Good, R. J.; Neumann, A. W. J. Colloid Interface Sci. 1975, 53, 235-248. (8) Furmidge, C. G. SJ. Colloid Sci. 1962, 17, 309-. (9) Johnson, R. E.; Dettre, R. H. J. Phys. Chem. 1964, 68, 1744-1750. (10) Li, D.; Neumann, A. W. Colloid Poly. Sci. 1992, 270, 498-504. (11) Marmur, A. J. Colloid Interface Sci. 1994, 168, 40-46. (12) Marmur, A. AdV. Colloid Interface Sci. 1994, 50, 121-141. (13) Marmur, A. J. Imaging Sci. Technol. 2000, 44, 406-409. (14) Miwa, M.; Nakajima, A.; Fujishima, A.; Hashimoto, K.; Watanabe, T. Langmuir 2000, 16, 5754-5760. (15) Neumann, A. W.; Good, R. J. J. Colloid Interface Sci. 1972, 38, 341-. (16) Oner, D.; McCarthy, T. J.Langmuir 2000, 16, 7777-7782. (17) Schwartz, A. M.; Tejada, S. B. J. Colloid Interface Sci. 1972, 38, 359-. (18) Schwartz, L. W.; Garoff, S. Langmuir 1985, 1, 219-230. (19) Schwartz, L. W.; Garoff, S. J. Colloid Interface Sci. 1985, 106, 422-437. (20) Cassie, A. B. D.; Baxter, S. Trans. Faraday Soc. 1944, 40, 0546-0550. (21) Cassie, A. B. D. Contact Angles. Discuss. Faraday Soc. 1948, 3, 11-16. (22) Drelich, J.; Miller, J. D. Langmuir 1993, 9, 619-621. (23) Good, R. J. A J. Am. Chem. Soc. 1952, 74, 5041-5042. (24) Henderson, J. R. Mol. Phys. 2000, 98, 677-681. (25) Israelachvili, J. N.; Gee, M. L. Langmuir 1989, 5, 288-289. (26) Wenzel, R. N. Ind. Eng. Chem. 1936, 28, 988-994. (27) Young, T. Philos. Trans. R. Soc., Part 1 1805, 65-87.

Derjaguin4 first introduced an important conception of metastable states into the study of the hysteresis, which is believed to be the origin of hysteresis. The droplets pin at different metastable states, and various contact angles are observed. In the study of contact angle and hysteresis, two kinds of models are often used: the sessile drop model and the vertical plate model. Johnson and Dettre9 have discussed the hysteresis of concentric ring models for the sessile drop system, and Neumann and Good15 have investigated the hysteresis of vertical plates with horizontal strips. They showed that heterogeneity may cause many metastable states and hence a number of possible stationary contact angles. Schwartz and Garoff18 investigated vertical plates with patchwise heterogeneous surfaces and formulated the free energy theoretically. The minimization of the free energy is the criterion for the determination of the contact angle. Because the full problem is of such complexity that only numerical solutions are feasible in general, they described two approximate solutions to the energy-minimization problem. They found that different patch structures (squares, circles, etc.) lead to different hysteresis values. However, in Schwartz’s study, they did not classify the heterogeneous surfaces by topology. From the topology point of view, continuous solid substrates (CSS) and discontinuous solid substrates (DSS) are two typical topology substrates. They are very simple with a continuous solid phase and a gas phase isolated on the surface or a continuous gas phase and a solid phase isolated on the surface. They are therefore completely different but appear to be fundamental with very basic patterns that can be extracted from natural complex patterned surfaces. We believe that we must understand the hysteresis behavior of the two first before we get into more complicated patterned substrates. In practice, they have been recognized to be two important patterns used in experiments. For DSS, many kinds of pillared ultrahydrophobic surfaces have been prepared.1,2,16,28,29 For CSS, grid surfaces have also been discussed widely.1,2 McCarthy et al.2,16 paid great attention to the hysteresis of the two types of patterns in their experiments. They obtained many interesting results. However, their effort did not show sufficient evidence for the hysteresis of CSS. We propose to use computer simulation to study such an object by using the dissipative particle dynamics (DPD) method.30-32 DPD is a promising approach to the study of hydrodynamic (28) Feng, L.; Li, S. H.; Li, H. J.; Zhai, J.; Song, Y. L.; Jiang, L.; Zhu, D. B. Angew. Chem., Int. Ed. 2002, 41, 1221-.

10.1021/la051983m CCC: $33.50 © 2006 American Chemical Society Published on Web 01/26/2006

2066 Langmuir, Vol. 22, No. 5, 2006

behavior on the mesoscopic scales. It has been applied in the simulation of droplet deformation in gravitational and shear fields,33-35 which show that it is suitable for studying the capillary phenomenon. The advantage of DPD is that momentum is conserved locally and long-range hydrodynamic interactions are preserved. Many factors make the measurement of the contact angle hysteresis in experiments more complex than one can image (e.g., gravity, roughness, liquid sorption and evaporation, liquid penetration, and surface swelling will have significant effects on the results, which is the cause of the poor reproduction of capillary experiments36,37). In simulations, those unexpected factors that may confuse the measurement could be excluded definitely, as shown in section 3.5. In the present study, we mainly investigated the dynamic hysteresis behaviors of the two kinds of solid/air composite surfaces: DSS and CSS. The homogeneous surface and the symmetrical sector composite surface were simulated first to validate that the DPD method is suitable for the study of contact angles. A series of investigations on the two topologies, CSS and DSS, were then performed in the present study. It was found that the two topologies have remarkably different dynamic behaviors, respectively. The crucial quantitysthe contact angle hysteresis θhyst for each topological substrateswas measured. The influence of temperature on θhyst was determined. For making ultrahydrophobic or ultralyophobic surfaces, which topology is suitable for such a purpose? It was discussed here according to the obtained results.

2. Background There are two aspects of the spreading and wetting of a liquid droplet on surfaces. One is the kinetics of the spreading (e.g., the growth of the lateral droplet radius and the contact angle as a function of time). Experiments and simulations have revealed universal laws such as Tanner’s spreading law or R ≈ t1/2 scaling of the contact radius and so on.38-50 The other is the thermodynamics of the capillary system,9,20-27,51-63 (i.e., the contact angle and contact angle hysteresis of droplets on surfaces). In (29) Feng, L.; Song, Y. L.; Zhai, J.; Liu, B. Q.; Xu, J.; Jiang, L.; Zhu, D. B. Angew. Chem., Int. Ed. 2003, 42, 800-802. (30) Espanol, P.; Warren, P. Europhys. Lett. 1995, 30, 191-196. (31) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155-160. (32) Koelman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, 21), 363-368. (33) Clark, A. T.; Lal, M.; Ruddock, J. N.; Warren, P. B. Langmuir 2000, 16, 6342-6350. (34) Jones, J. L.; Lal, M.; Ruddock, J. N.; Spenley, N. A. Faraday Discuss. 1999, 112, 129-142. (35) Liu, Y.; Kong, B.; Yang, X. Z. Polymer 2005, 46, 2811-2816. (36) Erbil, H. Y.; McHale, G.; Rowan, S. M.; Newton, M. I. Langmuir 1999, 15, 7378-7385. (37) Kwok, D. Y.; Neumann, A. W. AdV. Colloid Interface Sci. 1999, 81, 167-249. (38) Bacri, L.; Debregeas, G.; Brochard Wyart, F. Langmuir 1996, 12, 67086711. (39) Bekink, S.; Karaborni, S.; Verbist, G.; Esselink, K.Phys. ReV. Lett. 1996, 76, 3766-3769. (40) Deconinck, J.; Dortona, U.; Koplik, J.; Banavar, J. R. Phys. ReV. Lett. 1995, 74, 928-931. (41) Dortona, U.; Deconinck, J.; Koplik, J.; Banavar, J. R. Phys. ReV. E 1996, 53, 562-569. (42) Duffy, B. R.; Wilson, S. K. Appl. Math. Lett. 1997, 10, 63-68. (43) Garnier, G.; Bertin, M.; Smrckova, M. Langmuir 1999, 15, 7863-7869. (44) Heine, D. R.; Grest, G. S.; Webb, E. B. Phys. ReV. E 2003, 68. (45) Milchev, A.; Binder, K. J. Chem. Phys. 2002, 116, 7691-7694. (46) Nieminen, J. A.; Abraham, D. B.; Karttunen, M.; Kaski, K. MolecularDynamics of A Microscopic Droplet on Solid-Surface. Phys. ReV. Lett. 1992, 69, 124-127. (47) Nieminen, J. A.; Alanissila, T. Phys. ReV. E 1994, 49, 4228-4236. (48) Tanner, L. H. J. Phys. D: Appl. Phys. 1979, 12, 1473-. (49) Voue, M.; De Coninck, J. Acta Mater. 2000, 48, 4405-4417. (50) Wu, X. H.; Phan-Thien, N.; Fan, X. J.; Ng, T. Y. Phys. Fluids 2003, 15, 1357-1362.

Kong and Yang

this article, we focus on the thermodynamics, especially the contact angle hysteresis of liquid droplets on air/solid composite surfaces. The contortion of the three-phase line is an interesting subject for droplet wetting on nonideal surfaces, and it may affect θe and θhyst.5,15,18,19,22,64,65 In sections 4.2.2 and 4.3.2, we can see that the results are in close agreement with the Cassie equation and the contortion of the three-phase line has no remarkable effect. Hence, the effect of the contortion of the three-phase line is negligible in the DPD simulation system used here, and it is not considered in this article. For an ideal solid-liquid-vapor capillary system with a surface that is chemically homogeneous, rigid, flat, on the atomic scale, and free from chemical interaction by vapor or liquid adsorption, Young’s equation27 describes the relationship between the contact angle and interfacial tension of the thermodynamic equilibrium condition.

cos θy )

(γsv - γls) γlv

(1)

where θy is the equilibrium contact angle (θe) of the ideal surface, which can also be called the intrinsic contact angle (θi). γlv, γsv, and γsv are the liquid-vapor, solid-vapor, and solid-liquid interfacial tensions, respectively. For a solid/air composite surface,9,15,20-22,24,25 the Cassie equation20,21 is often used to describe the thermodynamic equilibrium condition

cos θc ) f1 cos θe1 - f2

(2)

where f1 and f2 are the fractional surface areas of the solid and air regions of surfaces and f1 + f2 ) 1. θc may be called the Cassie contact angle, which is the θe of the composite surface. θe1 is the intrinsic contact angle of the solid material. According to the Cassie equation, the contact angle of the heterogeneous surface is determined by f1. Because contact angle hysteresis is caused by metastable states,4 providing particular drop energy may lead to different hysteresis.4,9,15 In 1996, Decker and Garoff3 probed the energy barriers of a random heterogeneous surface in experiments with mechanical vibrational noise. In a DPD system, temperature represents the average kinetic energy of the particles in the system,66 and setting different temperatures gives the droplet (51) Adao, M. H.; de Ruijter, M.; Voue, M.; De Coninck, J. Phys. ReV. E 1999, 59, 746-750. (52) Bresme, F.; Quirke, N. J. Chem. Phys. 1999, 110, 3536-3547. (53) Iwahara, D.; Shinto, H.; Miyahara, M.; Higashitani, K. Langmuir 2003, 19, 9086-9093. (54) Lundgren, M.; Allan, N. L.; Cosgrove, T. Langmuir 2002, 18, 1046210466. (55) Lundgren, M.; Allan, N. L.; Cosgrove, T.; George, N. Langmuir 2003, 19, 7127-7129. (56) Milchev, A.; Milchev, A.; Binder, K. Comput. Phys. Commun. 2002, 146, 38-53. (57) Pal, S.; Weiss, H.; Keller, H.; Muller-Plathe, F. angmuir 2005, 21, 36993709. (58) Pesheva, N.; De Coninck, J. Phys. ReV. E 2004, 70. (59) Samsonov, V. M.; Dronnikov, V. V.; Volnukhina, A. A.; Muravyev, S. D. Surf. Sci. 2003, 532, 560-566. (60) Schneemilch, M.; Quirke, N.; Henderson, J. R. J. Chem. Phys. 2003, 118, 816-829. (61) Tang, J. Z.; Harris, J. G. J. Chem. Phys. 1995, 103, 8201-8208. (62) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Noca, F.; Koumoutsakos, P. Nano Lett. 2001, 1, 697-702. (63) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345-1352. (64) Boruvka, L.; Neumann, A. W. J. Chem. Phys. 1977, 66, 5464-5476. (65) Boruvka, L.; Neumann, A. W. J. Colloid Interface Sci. 1978, 65, 315330. (66) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987.

DissipatiVe Particle Dynamics Simulation

Langmuir, Vol. 22, No. 5, 2006 2067

different kinetic energy. In this article, we probed the metastable states of different topologies by varying the temperature (T).

In the interest of preserving simplicity in the model, ωij values are assumed to have the following functional form:

3. Methodology 3.1. Introduction of the DPD Method. DPD is based on the idea of thermostating a simulation “box” by invoking dissipative and fluctuation forces, which counterbalance each other to ensure constant thermal energy in the system. In addition, a conservative force is included. In the DPD method, the particles move in free space under the three forces acting on them. The dynamics of the system is computed at a large number of successive time steps by solving Newton’s equation of motion numerically for each particle. The total force acting on a particle i at any instant is the sum of three forces, dissipation (FDi ), fluctuation (FRi ), and conservation (FCi ),

Fi ) FDi + FRi + FCi

(3)

which can be written in another form

Fi )

∑j

(FDij + FRij + FCij )

(4)

where the summation ∑j runs over all momentum carriers except i. The dissipative force between a pair of particles (i and j) is assumed to be directly related to their relative velocity

FDij ) -γωDij [eij‚(νi - νj)]eij

(5)

where γ is a scalar coefficient characteristic of the system. ωDij is the weight function of rij, the distance between i and j. eij is a unit vector coincident with rij. FRij , the fluctuation force on i due to j, may be written as

FRij

)

σζijωRij eij

(6)

where ζij is a random number drawn from a Gaussian distribution whose first moment is equal to zero and whose variance is equal to unity. ωRij is the weight function corresponding to the fluctuation forces. The coefficient σ defines the variance of the fluctuation force, which must scale as 1/xδt, where δt is a small time interval during which the force would remain constant. A simple but realistic functional form of FCij , which is universally assumed in mesoscopic models, is

FCij ) RωRij ‚eij

(7)

where ωCij is the weight function and R is a parameter. A DPD model is completely characterized in terms of parameters γ, R, and σ and weight functions ωDij , ωRij , and ωCij . The relationship between σ and γ is furnished by the fluctuationdissipation theorem30

σ2 ) kT 2γ

(8)

with

(ωRij )2 ) ωDij

(9)

Taking kT as the unit of energy, eq 8 gives

σ2 ) 2γ

(10)

ωRij

)

(ωDij )1/2

)

ωCij

) ωij

{[

) 1-

)0

]

rij for rij < rc rc (11) for rij > rc

rc, serving as the unit of length, is the “cutoff” distance above which all three forces vanish. Equation 11 represents a repulsive, springlike, truncated conservation force, which means that the momentum carriers are assumed to be soft spheres. The time evolution of DPD particles is described by Newton’s law:

dr ) Vi dt

(12)

1 dVi ) (FCij dt + FDij dt + FRij xdt) m

(13)

The equation of motion was integrated numerically for each particle in the box using the so-called DPD-VV scheme.67-69 DPD-VV is based on the standard molecular dynamics velocity Verlet algorithm,66,70,71 which is a time-reversible and symplectic second-order integration scheme. The standard velocity Verlet has been shown to be relatively accurate in typical MD simulations, especially for large time steps.66 The DPD-VV scheme is presented succinctly by the following set of steps:

(1) Vi r Vi +

11 C (F ∆t + FDi ∆t + FRi x∆t) 2m i

(2) ri r ri + Vi∆t (3) calculate FCi , FDi , and FRi (4) Vi r Vi +

11 C (F ∆t + FDi ∆t + FRi x∆t) 2m i

(5) calculate FDi (1) DPD-VV differs from the standard velocity Verlet scheme in one important respect. The dissipative forces in DPD depend on the velocities, which in turn are governed by the dissipative forces (eqs 5 and 13). This matter is not accounted for by the standard velocity Verlet scheme. DPD-VV, however, accounts for this complication in an approximate fashion by updating the dissipative forces for a second time at the end of each integration step. This improves its performance considerably while keeping it computationally efficient because the additional update of dissipative forces is not particularly time-consuming. 3.2. Models. The present systems are composed of three phases: two immiscible fluids (1, 2) and a solid (3). The solid phase is “created” by freezing the motion of the particles so that they are locked in a certain configuration. As a feature, DPD is good for fluids.30-35 Because the DPD conservation interaction consists of only the repulsion force without the attraction force, it is hard to immediately describe a system consisting of matter with different densities, for instance, the gas/liquid system. As (67) Besold, G.; Vattulainen, I.; Karttunen, M.; Polson, J. M. Phys. ReV. E 2000, 62, R7611-R7614. (68) Nikunen, P.; Karttunen, M.; Vattulainen, I. Comput. Phys. Commun. 2003, 153, 407-423. (69) Vattulainen, I.; Karttunen, M.; Besold, G.; Polson, J. M. J. Chem. Phys. 2002, 116, 3967-3979. (70) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: New York, 2002. (71) Verlet, L. Phys. ReV. 1967, 159, 98-.

2068 Langmuir, Vol. 22, No. 5, 2006

Figure 1. Part of an initial model. Particles are located regularly at the corners of simple cubic lattices.

Figure 2. Schematic of θa and θr. (a) Initial model for θr, θini < θr. (b) Metastable state of θr. (c) Metastable state of θa. (d) Initial model for θa, θini > θa.

is well known, gas is a fluid. Omitting the viscosity effect on the dynamics, we adopted a liquid to be the gas phase in the present study. Now, the modeled system of fluid1/fluid2/solid describes a real gas (or vapor)/liquid/solid system. According to the derivation of the thermodynamic equations,9,20-27 the equations of the liquid/liquid/solid and gas (or vapor)/liquid/ solid systems do have the same forms. In this article, for convenience, we regard the two immiscible fluids in our systems as vapor (1) and liquid (2) and use traditional nomenclature (γlv, γsv, and γsv) in the following discussions. Simulations are carried out with the usual periodic boundary conditions. The center box size is 50 × 50 × 50. In initial models, each of the particles was located regularly at the corners of simple cubic lattices (Figure 1). The box can be seen as two parts. The lower part is a substrate with a density of 12, and it is composed of 6 layers of particles in which the distance between every set of nearest neighbors is 0.44. The upper part is fluid with a density of 6, and it is composed of 85 layers of particles in which the distance between nearest neighbors is 0.56. The total number of DPD particles in the system is 766 476. Figure 1 shows part of an initial model. In the upper part, the liquid phase is introduced as a drop located on the substrate (Figure 3). Particles in the remaining space of the upper part are set to be vapor. In the substrate, four kinds of patterns are designed for different simulations, including a homogeneous solid, a symmetrical sector, and CSS and DSS patterns, as shown in Figure 3. For each kind of pattern except the homogeneous surface, the fractional surface areas of the solid were adjusted by varying a certain number of particles to be solid or vapor. Each simulation was carried out for more than 60 000 steps. The contact angles value were averaged after 20 000 steps. 3.3. Force Parameters. In principle, 18 parameters would be needed to characterize all of the like and unlike interactions in such a system: they are γ11, γ22, γ33, γ12, γ13, γ23, σ11, σ22, σ33, σ12, σ13, σ23, R11, R22, R33, R12, R13, and R23. However, this number is reduced to 12 by virtue of the fluctuation-dissipation theorem (eq 10). The conservation forces were chosen as follows, where phases 1-3 are vapor, liquid, and solid, respectively. We took R11 ) R22 ) 12.5, which means that the vapor and liquid are

Kong and Yang

assumed to have identical isothermal compressibility. The vapor/liquid interaction R12 was taken to be 80, which ensures that the two fluids are immiscible. Vapor/solid and liquid/solid parameters R13 and R23, respectively, are selected according to section 4.1. Groot72 has carried out an extensive investigation of the effect of the magnitude of δt on the temperature stability of the system. He concluded that the instability problem could be overcome only by reducing the time step to 0.05 DPD units or to even smaller values if practical. In the same computations, Groot established through trial and error that the optimum value of σ that met the requirements of fast temperature equilibration, rapid convergence, and stability is around 3.35 for δt in the range from 0.04 to 0.05. In the present simulations, we have followed Groot’s recommendation34,72 in selecting the values for δt and σ (δt ) 0.05 and σ11 ) σ22 ) σ12 ) σ13 ) σ23 ) 3.35). Hence, γ11 ) γ22 ) γ12 ) γ13 ) γ23 ) 5.61. 3.4. Measurement of Contact Angle. There are three main techniques for measuring the contact angle and contact angle hysteresis in experiments.36 They are the sessile drop (or captive bubble in a liquid),36,73-76 Wilhelmy plate,77 and inclined plane methods.8,78-82 In our study, for the convenience of computer simulation, an alternative method analogous to the sessile drop method is used to obtain θa and θr. In this method, a liquid droplet is placed on a substrate with an initial contact angle (θini) greater than or less than its θe. The drop will spontaneous vary its contact angle toward θe until it is pinned at a metastable state, as shown in Figure 2 where the contact angle is θa or θr. There is no volume variation of the drop in simulation, which is the main difference between this method and the sessile drop method in experiments. A similar method has been used in the molecular dynamics simulation by Adao et al.51 The values of contact angles were calculated in a similar manner in ref 34. Under equilibrium conditions, a drop exists as an approximate spherical cap on a surface in the absence of gravity. We can determine the contact angle by fitting a spherical segment to the simulated shape. The easiest way to achieve this is to calculate in a simulation run the average number of particles in a cross section of the drop as a function of height, z, above the surface. For a sphere, the number of particles per unit height should be equal to πF[(R2 - z02)+2zz0 - z2], where R is the radius of the sphere, z0 is the distance of the center of the sphere from the surface, and F is the particle density. The contact angle is simply equal to 90 + sin-1(z0/R) degrees. 3.5. Simulation Snapshots. Figure 3 shows snapshots of simulations with the four topologies. From the snapshots of the end states, none of the models showed any liquid particle leaving the droplet or intruding into the cavities in these surfaces. Therefore, the inevitable liquid sorption and volatilization in experiments are excluded in the simulations. This property is (72) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423-4435. (73) Drelich, J.; Miller, J. D.; Good, R. J. J. Colloid Interface Sci. 1996, 179, 37-50. (74) Li, D.; Cheng, P.; Neumann, A. W. AdV. Colloid Interface Sci. 1992, 39, 347-382. (75) Lin, S. Y.; Chang, H. C.; Lin, L. W.; Huang, P. Y. ReV. Sci. Instrum. 1996, 67, 2852-2858. (76) Timmons, C. O.; Zisman, W. A. J. Colloid Interface Sci. 1966, 22, 165-. (77) Neumann, A. W.; Tanner, W. J. Colloid Interface Sci. 1970, 34, 1-. (78) Bikerman, J. J. J. Colloid Sci. 1950, 5, 349-359. (79) Extrand, C. W.; Kumagai, Y. J. Colloid Interface Sci. 1995, 170, 515521. (80) Extrand, C. W.; Kumagai, Y. J. Colloid Interface Sci. 1997, 191, 378383. (81) Kawasaki, K. J. Colloid Sci. 1960, 15, 402-407. (82) Rotenberg, Y.; Boruvka, L.; Neumann, A. W. J. Colloid Interface Sci. 1984, 102, 424-434.

DissipatiVe Particle Dynamics Simulation

Langmuir, Vol. 22, No. 5, 2006 2069

Figure 3. Snapshots of simulations. Vapor particles are not shown for clarity. T ) 1.0. (a, d, g, and j) Top views of the four kinds of substrates. (b, e, h, and k) Side views of initial states. (c, f, i, and l) Side views of end states. (a-c) Homogeneous substrate, θini ) 90°, θr ) 119.8°. (d-f) Symmetrical sector composite substrate, f1 ) 0.5, θini ) 120°, θr ) 137.8°. (g-i) CSS patterned substrate, f1 ) 0.67, θini ) 120°, θr ) 132.2°. (j-l) DSS patterned substrate, f1 ) 0.67, θini close to 180°, θa ) 133.5°.

mainly due to the conservation forces parameters between the different phases we used in the simulations. If these parameters are small enough, then remarkably liquid sorption and volatilization will take place. No other force exists in the simulation except the three counterbalance forces between particles, thus there is no gravity effect.

4. Results and Discussion In section 4.1, the relations between the interfacial tension and parameters of conservation force, interfacial tension and temperature, and θe and temperature were obtained by simulation and calculation. In section 4.2, simulations were carried out at

2070 Langmuir, Vol. 22, No. 5, 2006

Kong and Yang

high temperature (T ) 1.0), and in section 4.3, simulations were carried out at low temperature (T ) 0.3). In section 4.4, simulations were carried out at varied temperature and constant f1. In section 4.5, we show some comparisons with other results and conclusions. In sections 4.2 and 4.3, the homogeneous surface and the symmetrical sector heterogeneous surfaces were simulated first. The former corresponds to the ideal smooth homogeneous surface for which the Young equation should be satisfied, and the latter corresponds to the ideal smooth heterogeneous surface for which the Cassie equation should be satisfied. CSS and DSS surfaces were also simulated to find the relationship between topology and contact angle hysteresis. 4.1. Relationships of Interfacial Tension versus the Conservation Force Parameter and Temperature and θe versus Temperature. Simulations were carried out in a manner similar to that described by Jones et al.34 to obtain the relationship between the parameters of the conservation force and interfacial tension. In simulations, the fluid(1)/fluid(2) interfacial tension, γ12, is conveniently computed using the Irving-Kirkwood equation83

γ12 )

∫[pzz - 21(pxx + pyy)] dz

(14)

where pxx, pyy, and pzz are the three diagonal components of the pressure tensor. The interface is parallel to the x-y plane. In terms of the interparticle forces, eq 14 is expressed as

γ12 )

1

[∑

A

i,j

1 Fijzrijz - (Fijxrijx + Fijyrijy) 2

]

i > j (15)

(Fijx, Fijy, Fijz) and (rijx, rijy, rijz) are, respectively, the x, y, and z components of the force and the separation vectors between a pair of particles i and j. The summation is carried out over all particle pairs. A is the interfacial area. In simulations, the fluid/solid interfacial tension, γFS (i.e., γ13 or γ23), is given by the expression

γFS )

1

[∑

A

i,j

]

1 1 Fijzrijz - (Fijxrijx + Fijyrijy) + 2 A

∑i Fisris

i > j (16)

The second sum is taken over all particles. Fis is the force on i due to its interaction with the solid particles, and ris is the distance of i from the solid boundary. We performed simulations on systems of fluid/fluid and fluid/solid with a planar interface and computed the ensemble averages of the summations in eqs 15 and 16 to obtain the respective interfacial tensions. The results are presented in Figure 4a, which is consistent with a recent report.34 Figure 4a shows that, at a particular temperature, the interfacial tensions increase with the increase in the conservation force parameter and the solid/fluid interfacial tension is greater than the liquid/fluid interfacial tension. With the set of conservation force parameters, interfacial tensions can be determined according to Figure 4a. In the present study, thes conservation force parameter are always taken to be R12 ) 80, R13 ) 20, and R23 ) 91. When T ) 1.0, the corresponding interfacial tensions are γlv ) 12.666, γls ) 34.1218, and γsv ) 27.7888, and according to Young’s equation, θy ) 120°. Temperature is an important parameter of the simulation system. In general, when temperature is varied, θi of a homogeneous surface may be also varied,84 thus the effect of (83) Irving, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817-829.

Figure 4. (a) Relationship between interfacial tensions and parameters of the conservation force. (b) Effect of temperature on θe.

temperature on θi must be considered. θi values calculated according to Figure 4a and Young’s equation are shown in Figure 4b. Figure 4b shows that although the interfacial tensions vary with temperature θi undergoes no remarkable change (θi ) 120 at T ) 1.0 and θi ) 119.1 at T ) 0.3). This property of the DPD simulation system used here is very important, and we will not consider the effect of temperature on θi in the following discussions. This property is also common in real systems.84 4.2. High-Temperature Simulation and No Hysteresis Behavior. For the four topologies of substrates, simulations were carried out first at T ) 1.0. 4.2.1. Homogeneous Substrate. Two simulations were done with a homogeneous model to obtain θa and θr. Figure 3a-c shows the model of θr. The results show that θr ) 119.8 and θa ) 120.4, which are consistent with the expected value according to Young’s equation (120°). This indicates that the DPD method is applicable to the study of the capillary system and that the intrinsic contact angle can be controlled by selecting certain parameters of the conservation force. 4.2.2. Symmetrical Sector Composite Substrates. A series of symmetrical sector composite models were simulated from f1 ) 0.11111 to 0.8889 by changing the ratio of solid sectors. Figure 3d-f shows the model of θr for f1 ) 0.5. The results for θa and θr are shown in Figure 5a. The line in Figure 5a is the fitting line according to the Cassie equation, and the fitting value of θe1 is 122.8°. It shows that both the θa and θr of the symmetrical sector composite substrates are consistent with the expected value of the Cassie equation. This means that there is no hysteresis effect for symmetrical sector composite substrates at T ) 1.0. The good consistency shows that the DPD method is an appropriate approach to studying the contact angle of the solid/air composite surface. 4.2.3. CSS and DSS Pattern Substrates. Two series of CSS and DSS pattern models were built by changing f1, the ratio of solid to air, in the substrates. Because of the difficulty in making such a tiny geometrical pattern, the range of f1 that we can make (84) Neumann, A. W. AdV. Colloid Interface Sci. 1974, 4, 105-191.

DissipatiVe Particle Dynamics Simulation

Langmuir, Vol. 22, No. 5, 2006 2071

Figure 5. High-temperature simulation, T ) 1.0. Contact angle as a function of f1. b and O are θa, and [ and ] are θr. (a) Symmetrical sector composite substrates. The solid line represents the fitted Cassie’s law, and the correlation coefficient is 0.99922, with θe1 ) 122.8°. (b) CSS and DSS patterned substrates. Open symbols represent DSS, and solid symbols represent CSS. The solid line represents the fitted Cassie’s law, and the correlation coefficient is 0.96187, with θe1 ) 120.7°.

Figure 6. Low-temperature simulation, T ) 0.3. Contact angle as a function of f1. b and O are θa, and [ and ] are θr. (a) Symmetrical sector composite substrates. The solid line represents the fitted Cassie’s law, and the correlation coefficient is 0.99967, with θe1 ) 121.0°. (b) CSS and DSS patterned substrates. Open symbols represent DSS, and solid symbols represent CSS. A Cassie line with θe1 ) 121.0° is plotted for comparison.

is from 0.32 to 0.93 for the CSS models and from 0.07 to 0.67 for the DSS models. They basically cover the range from 0 to 1. Figure 3g-i shows the CSS model of θr for f1 ) 0.67 and j-l shows the DSS model of θa for f1 ) 0.67. Figure 5b shows the results of θa and θr for the two series. The line is the fitting line according to the Cassie equation, and the fitting value of θe1 is 120.7°. It can be seen that for both kinds of patterns θa and θr are located near the Cassie line. Although the deviations in CSS and DSS are greater than those of symmetrical sector surfaces, they are still approximately consistent with the Cassie line and there is no hysteresis effect. The greater errors may be caused by the asymmetrical nature of the two kinds of patterns. The origin of the deviations will be discussed further in the first paragraph in section 4.5. In summary, θe can be reached, and the Cassie equation is satisfied approximately at T ) 1.0. Hysteresis is not observed in the simulations. 4.3. Low-Temperature Simulation and Hysteresis Behavior. For the four kinds of substrate topology, simulations were also carried out at T) 0.3. 4.3.1. Homogeneous Substrates. For homogeneous substrates, the results are θr ) 118.2° and θa ) 118.4°, which are very close to the expected value of 119.1° according to Young’s equation. These results testify to the conclusion in section 4.1 that for the DPD system we used here temperature has no definite influence on θi. 4.3.2. Symmetrical Sector Composite Substrates. For symmetrical sector composite substrates, the results are shown in Figure 6a. The line is the fitting line according to the Cassie equation, and the fitting value of θe1 is 121.0°. Figure 6a shows that both θa and θr are still consistent with the Cassie equation. This is reasonable because for this kind of model there are no energy barriers when the droplet advances and recedes, as analyzed by Johnson and Dettre.9 There should be no metastable states or hysteresis, even if the drop energy (temperature) is very low.

4.3.3. CSS and DSS Patterned Substrates. For CSS and DSS patterned substrates, the results at T ) 0.3 (Figure 6b) are remarkably different from those at T ) 1.0 (Figure 5b). The line in Figure 6b is the Cassie line with θe1 )121.0° (the fitting value at T ) 0.3 for the symmetrical sector composite substrate, Figure 6a). It can be seen that θa of CSS, θr of CSS, and θa of DSS obviously deviate from the Cassie equation. Compared with the results at T ) 1.0, we can see that at high temperature (Figure 5b) for each simulation the results are located on both sides of Cassie line but at low temperature (Figure 6b) they are all higher or lower than the Cassie line, which means that the Cassie equation is not satisfied for the three kinds of contact angles. θa of DSS is very interesting; it remains around 160° and does not vary with f1. For example, when f1 ) 0.67, according to the Cassie equation θc should be 131.7° when θe1 )121.0°. Obviously there is a metastable state near 160°. At T ) 0.3, the droplet does not have enough energy to overcome the barriers between the metastable state and the thermodynamic equilibrium, thus it remains in the metastable state. The hysteresis for DSS is significant. θr of DSS is close to the Cassie line. In summary, the results show that at low temperature CSS is very different from DSS. The distribution of metastable states is possibly the origin of these differences. Thus we probed the distribution of metastable states by varying the temperature in the next section. 4.4. Metastable States Probed by Varying the Temperature. To understand the relationship between hysteresis and topology more clearly, we probed the metastable states by varying the temperature (T ) 0.3 to 1.0) with constant f1. At different temperatures, the droplets have different energies, and different metastable states may be reached. First, simulations were carried out for homogeneous and symmetrical sector composite substrates (at f1 ) 0.5). The contact angles undergo no obvious change when T is varied from 0.3 to 1.0. The results show that there is no metastable state and no hysteresis for these two substrate topologies, as discussed in section 4.3.

2072 Langmuir, Vol. 22, No. 5, 2006

Figure 7. Variation of temperature to probe the metastable states of CSS and DSS, with f1 ) 0.67. Open symbols represent DSS, and solid symbols represent CSS. (a) Contact angles at different temperatures. There are many steps on the curves, corresponding to metastable states. b and O are θa, and [ and ] are θr. (b) θhyst at different temperatures. θhyst of DSS is greater than that of CSS. θhyst trends to 0 when the temperature increases. (c) Distribution of metastable states. Contact angles (representing the positions of metastable states) vs the lowest temperatures needed to reach the metastable states. b and O are θa, and [ and ] are θr. The arrows between metastable states indicate the wetting direction of the droplets. The solid line represents DSS, and the dotted line represents CSS. For the receding of DSS (]), only one metastable state was found, thus there is no moving direction for it.

Second, simulations were carried out for for CSS and DSS at f1 ) 0.67, which corresponds to θc ) 131.7° according to the Cassie equation when θel ) 120°. The results are shown in Figure 7a. Figure 7a shows that there are many steps when the temperature is varied, which represent metastable states at such positions. Obviously, there should be an energy barrier between every two metastable states. We can see that when the temperature is lower (close to 0.3) the droplet has less energy and may be pinned at a contact angle far away from θe. When the temperature is higher (close to 1), the droplet has greater energy and may go through some barriers to reach a contact angle much closer to θe. It was reasonable to imagine that if the temperature was high enough (i.e., if the droplet had enough energy) then both the advancing and reducing processes might reach θe and hence the difference between θr and θa might vanish. θhyst values of CSS and DSS are plotted versus temperature in Figure 7b. It can be seen that θhyst of DSS is always greater

Kong and Yang

than that of CSS, hence CSS is more suitable for building an ultrahydrophobic or ultralyophobic surface than DSS from the point of view of dynamic wettability.2,8,14,16 θhyst trends to 0 when the temperature increases. It was reasonable to imagine that θhyst would reach 0 and then the difference for the two topologies would vanish if the temperature was high enough (i.e., if the droplet had enough energy). Contact angles of the metastable states are plotted versus the lowest temperatures needed to achieve them, as shown in Figure 7c. This shows that the distributions of metastable states of DSS and CSS have definite difference. For CSS, metastable states in advancing and receding processes are located beside θe approximately symmetrically and their amounts and distances to θe are similar. For DSS, metastable states in advancing and receding processes are definitely asymmetric, and there are several metastable states that are advancing but only one stable state that is receding. According to Figure 7, it is clear that the apparent contact angle is determined by the distribution of the metastable states and the energy of the droplet. 4.5. Comparisons with Other Results and Conclusions. (1) Adao et al.51 and Samsonov et al.59 have simulated the contact angle of droplets on square patchwise surfaces using molecular dynamics. Their results show that the Cassie equation is valid when the patches are small but invalid when the patches are too large. Obviously, f1 applied to the Cassie equation20 is a constant and does not change in the process of wetting. Therefore, it makes the plot of the f1 and cos θ linear. For Adao’s models, when the patches are small enough, f1 can be treated as approximately constant, and the Cassie equation is valid as expected. Otherwise, in Samsonov’s simulation the invalidation of the Cassie equation is inevitable. In the present study, we use symmetrical sector composite substrates to validate the Cassie equation, which ensures that f1 is constant during the wetting process. The results in sections 4.2.2 and 4.3.2 are in quantitative agreement with the Cassie equation and show that DPD is an appropriate approach for simulations of heterogeneous surfaces. In section 4.2.3, we have found that for CSS and DSS patterned substrates the deviations are much larger than that for symmetrical sector composite substrates but the contact angle is still around the expected value of the Cassie equation. This may be an intermediate case between those simulated by Adao and Samsonov. When the three-phase line moves on the surface, the value of f1 may be not a constant, and the value of f1 in Figure 5b is only an average f1 of the whole substrate.15 Fortunately, the plot is nearly linear. (2) Johnson and Dettre9 have analyzed the concentric ring model theoretically. They found that when drop energy is not high enough θa of such topology is close to θe and θr is much less than θe. In the present study, according to Figure 6b, we can see that the two topologies discussed in the present study (CSS and DSS) follow other rules. For DSS, θr is close to θe and θa is much greater than θe, whereas for CSS, θa and θr deviate from θe similarly. In addition, it is difficult to obtain information on the thermodynamic properties of the system (e.g., θe and interfacial tensions from observed apparent contact angles) unless the substrate is a homogeneous, symmetrical sector or has no hysteresis. (3) McCarthy2,16 has paid attention to the effect of topology on hysteresis. He suggested that for solid/air composite surfaces with a particular composition (f1) but different topologies, different contact angle hysteresis values should be observed. This has been verified in Figures 6b and 7b. He proposed that for CSS (PTFE screen) there should be metastable states in both process

DissipatiVe Particle Dynamics Simulation

Langmuir, Vol. 22, No. 5, 2006 2073

of advancing and receding, which is consistent with our results (Figure 7c). He also proposed that for DSS (PTFE posts) there should be little or no metastable states and the droplet could reach θe easily in both process of advancing and receding, which is not consistent with our results (Figure 7c). Figure 7c shows clearly that for DSS, in the receding process, McCarthy’s proposal is true but in the advancing process many metastable states do exist and they cause very large θa values when the drop does not have enough energy. McCarthy had done many experiments on the DSS pattern and has produced many ultrahydrophobic surfaces with high θa and θr and low θhyst, but no comparative result for the CSS pattern have been reported, probably because of difficulties with the experiments. McCarthy suggested that DSS is preferred for building ultrahydrophobic surfaces with low θhyst, yet from Figure 7b, we can see that θhyst of CSS is less than that of DSS and that CSS is more suitable for building ultrahydrophobic surfaces. (4) Bico and Marzolin1 have also done experiments to study the contact angles of CSS and DSS. In their report, they did not use the same f1 for the two topologies because of difficulties with the experiments. In their experiments, f1 is 0.64 for CSS and 0.05 for DSS. They found that θhyst of CSS is 60°, much greater than that of DSS, 15°. However they proved that the high hysteresis of CSS is caused by the liquid penetrating into the cavities of the surface in the process of receding. In this case, the solid/air composite surface was changed to a solid/liquid composite surface. According to the Cassie equation,20 its θc should be calculated as follows:

the present study, computer simulation was used, and liquid penetration was excluded, thus our results (θhyst of CSS is less than that of DSS) do not contradict Bico’s report.

cos θc ) f1 cos θe1 + f2

Acknowledgment. We thank the National Nature Science Foundation, 863 High Technology Project, and the Special Funds for Major State Basic Research Project (2004CB720606) for continuing financial support.

(17)

The value of θc in eq 17 is much less than that of eq 2. Hence, the high θhyst of CSS that they observed is not due to the topology but rather to the low θr caused by the experimental process. In

5. Conclusions In the present article, DPD simulation was used to investigate contact angle hysteresis of droplet wetting on patterned solid/air composite surfaces. Therefore, DPD has been shown to be an appropriate approach to studying contact angles of homogeneous and solid/air composite surfaces. In the present simulation (where temperature has no influence on θi), we drew the following conclusions. For the homogeneous and symmetrical sector composite substrates, it was validated that there is no metastable state and no contact angle hysteresis and that varying the drop energy (temperature) does not influence the apparent contact angle of the two substrates. For the CSS and DSS patterned substrates, topology has a great influence on the distribution of metastable states. At the same f1, DSS has a different distribution of metastable states from that of CSS. Therefore, the two topologies have different contact angle hysteresis values θhyst: DSS has a greater θhyst than does CSS. We found that at low temperature this difference in θhyst is larger and that when the energy (temperature) increases θr and θa trend to θe and θhyst trends to 0. Hence, CSS is more suitable for building an ultrahydrophobic or ultralyophobic surface than DSS, from the point of view of dynamic wettability.

LA051983M