Article pubs.acs.org/Langmuir
Pinning of the Contact Line during Evaporation on Heterogeneous Surfaces: Slowdown or Temporary Immobilization? Insights from a Nanoscale Study Jianguo Zhang,* Florian Müller-Plathe, and Frédéric Leroy* Eduard-Zintl-Institut für Anorganische und Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Straße 4, D-64287 Darmstadt, Germany S Supporting Information *
ABSTRACT: The question of the effect of surface heterogeneities on the evaporation of liquid droplets from solid surfaces is addressed through nonequilibrium molecular dynamics simulations. The mechanism behind contact line pinning which is still unclear is discussed in detail on the nanoscale. Model systems with the Lennard-Jones interaction potential were employed to study the evaporation of nanometer-sized cylindrical droplets from a flat surface. The heterogeneity of the surface was modeled through alternating stripes of equal width but two chemical types. The first type leads to a contact angle of 67°, and the other leads to a contact angle of 115°. The stripe width was varied between 2 and 20 liquid-particle diameters. On the surface with the narrowest stripes, evaporation occurred at constant contact angle as if the surface was homogeneous, with a value of the contact angle as predicted by the regular Cassie−Baxter equation. When the width was increased, the contact angle oscillated during evaporation between two boundaries whose values depend on the stripe width. The evaporation behavior was thus found to be a direct signature of the typical size of the surface heterogeneity domains. The contact angle both at equilibrium and during evaporation could be predicted from a local Cassie−Baxter equation in which the surface composition within a distance of seven fluid-particle diameters around the contact line was considered, confirming the local nature of the interactions that drive the wetting behavior of droplets. More importantly, we propose a nanoscale explanation of pinning during evaporation. Pinning should be interpreted as a drastic slowdown of the contact line dynamics rather than a complete immobilization of it during a transition between two contact angle boundaries. while shrinking in size. In the mixed mode,12,16,20 both the contact line length and contact angle change during evaporation. The most interesting pinning phenomenon is related to the constant contact line (CCL) mode,6,12,20 where the contact angle decreases while the contact line is immobilized and its length remains constant, i.e., the drop keeps its contact area or footprint, while becoming flatter. The CCL mode or pinning phenomenon is usually considered to arise from surface heterogeneities, e.g., topographical roughness and/or chemical heterogeneities.21 As a matter of fact, most surfaces in nature are heterogeneous, and some of the surfaces used in experiments and claimed to be homogeneous are only less heterogeneous on the molecular level than other surfaces. When a droplet is placed on a heterogeneous surface,22−32 many phenomena such as contact line pinning25 and contact angle hysteresis22,27,31 may be observed. These wetting properties influence the
1. INTRODUCTION Droplet evaporation on heterogeneous or patterned surfaces has attracted a lot of interest owing to new potential applications such as DNA stretching and DNA mapping,1,2 coating,3 inkjet printing,4 and nanopatterning.5 Understanding the fundamental mechanism of droplet evaporation on heterogeneous surfaces is a key step to controlling deposit patterns and producing smart interfaces or devices. However, the relationship between the macroscopic evaporation behavior of droplets and the nanoscale structure of surfaces is still elusive. In particular, the pinning mechanism behind the immobilization of the three-phase contact line is still unclear. Pinning is considered to be the main cause of the coffee-ring effect,6,7 and such a phenomenon should be suppressed to improve some applications such as inkjet printing.8 Droplet evaporation has been investigated by both experiments and theories.6,7,9−20 Different evaporation modes have been reported to describe the dynamics of the contact line and contact angle. In the constant contact angle (CCA) mode,10,12,17,20 the length of the contact line decreases while the contact angle remains constant, i.e., the drop keeps its shape © XXXX American Chemical Society
Received: March 24, 2015 Revised: June 4, 2015
A
DOI: 10.1021/acs.langmuir.5b01097 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir evaporation behavior. Kusumaatmaja and Yeomans32 carried out an extensive study through Lattice Boltzmann calculations of the behavior of droplets whose volume was quasistatically increased or decreased on heterogeneous surfaces. This work illustrated the richness of scenarios that may be met when evaporation occurs in a slow hydrodynamic mode. In their work, the contact line is a strictly one-dimensional object that may be distorted depending on the surface heterogeneities. Other works have taken a different perspective to consider the contact line and especially how it drives the droplet’s shape. It has been reported that when a droplet sits on a heterogeneous surface, the contact angle is determined by the wetting area adjacent to the three-phase contact line. This view is in line with the idea that the contact angle is a local wetting property determined by the intermolecular interactions in the vicinity of the contact line.27,33 For heterogeneous surfaces, the Cassie34 or Wenzel35 equation is usually employed instead of Young’s equation to calculate the static contact angle. However, these equations may not be valid when the distribution of heterogeneities on the surface is not uniform or when the droplet size is not large enough compared to the typical size of surface36−39 heterogeneities. These considerations indicate that the contact angle has a very close relationship with the nanoscale structure around the contact line. Thus, more complex evaporation behavior is expected when a droplet wets a nanoscale-structured surface. To date, a few studies of droplet evaporation on both structurally rough and chemically heterogeneous surfaces have been reported,20,21,26,39−41 and a pinning phenomenon and CCL−CCA transition have been found. The mechanisms behind these phenomena have been discussed in terms of pinning and depinning forces20,21 or the existence of an energy barrier between two stable contact-angle states.19,42,43 For example, Chen et al.20 calculated the pinning force Fp and depinning force Fd for a droplet on a microlevel-structured surface and claimed that pinning is caused by Fp > Fd. Orejon et al.43 proposed an intrinsic energy barrier, which comes from the surface heterogeneity and needs to be overcome before depinning of the contact line occurs, to explain the pinning. Note that these two explanations originate from a macroscale experiment. When considering the contact angle to be a local wetting property, a microscale picture of pinning is still missing. As we will discuss the pinning from the point of view of atomicscale interactions in the following text, we find that during pinning, the contact line moves slowly rather than being immobilized between two surface domains with different attraction for the fluid. In the present work, the evaporation of nanodroplets from a chemically striped surface is studied by molecular dynamics simulation. We proceed in a way similar to the experimental work of Extrand,27 where the volume of droplets was adjusted such that the contact line met the boundary between lyophobic and lyophilic domains. For simulation convenience, a cylindrical droplet is used in our simulations and the heterogeneity of the surface is modeled through alternating stripes of equal width of two types. The evaporation behavior was investigated by monitoring the dynamics of the contact angle and the contact line. The study focuses on characterizing when, where, and why pinning happens.
Figure 1. Cylindrical droplet on a substrate with alternating A and B stripes.
As in our previous work,17 the particles interact through the 12−6 Lennard-Jones potential. The total potential energy Upot of the system reads ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij 4 ε ∑ ∑ ij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ i=1 j>i ⎣⎝ ij ⎠ N
Upot =
(1)
Interaction parameters σl = 0.3405 nm and εl = 0.992 kJ/mol were used for the liquid−liquid interactions, and σs = 0.4085 nm and εs = 9.92 kJ/mol were used for the solid−solid interactions. All quantities are given in reduced units with respect to the Lennard-Jones parameters for the liquid particles. The mass of the liquid atoms was m = 6.63 × 10−26 kg. The solid atoms were 10 times heavier than the liquid ones. The value T* = kBT/ε = 1 of the reduced temperature corresponds to 119.8 K, i.e., the fluid roughly corresponds to Ar, and the solid atoms are slightly larger (σs = 1.2σl). The time unit is τ = (mσ2/ε)1/2. The solid is made up of two kinds of stripes A and B with equal width (Figure 1). To generate a wetting contrast between the neighboring stripes, different energy parameters εsl* = 0.4 and 0.7 were chosen for the fluid−solid A and fluid− solid B interactions, respectively. The solid−liquid distance parameter σsl* was 1.1. As reported in our previous work,17 the equilibrium contact angles on pure A and pure B surfaces are 115.0 ± 6.9 and 66.7 ± 4.2°, respectively, at T* = 0.83. The liquid consists of 10 800 atoms. These atoms were initially distributed on a cubic lattice at near-liquid density. The solid substrate was modeled as an fcc lattice containing six layers of 1000 atoms each. It was placed underneath the lattice of the liquid atoms. To study the stripe width effect on evaporation, the width was set to 1σs, 2σs, 6σs, 10σs, and 20σs. The corresponding systems are named d1, d2, d6, d10, and d20, respectively, in the following text. The rectangular periodic simulation box size was set to 100σs, 11.2σs, and 1468.8σs in the x, y, and z directions, respectively, where z is the surface normal. The Berendsen thermostat44 with a coupling time of 0.25τ* was used to control the temperature of the substrate atoms only to mimic the heating process. The simulations were carried out using LAMMPS.45 The time step and the cutoff distance were set to 0.0025τ* (5 fs) and 4.4σl (1.5 nm), respectively. Owing to the short extension of the system in the y direction, the droplet is periodic in this direction and has a cylindrical shape.
2. SIMULATION DETAILS The simulation setup is shown in Figure 1. The whole system consists of a solid substrate and of a liquid droplet in its vapor. B
DOI: 10.1021/acs.langmuir.5b01097 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir The droplet was initially relaxed at T* = 0.67 on a uniform B surface (εsl = 0.7) and was subsequently equilibrated at T* = 0.83. Then, some parts of the substrate were turned into A-type surfaces in order to build the periodic pattern of alternating A and B stripes (Figure 1). An additional 10 000τ run was performed to equilibrate the system on the patterned surface. Evaporation was then triggered by removing the gas molecules at a rate of 1 atom per 400 time steps, equal to 1 atom per 0.1τ. Trajectories were collected until the droplet disappeared. As mentioned above, the droplets are cylindrical. It is thus impossible to define a contact radius, and a cylindrical droplet has in fact two contact lines. We use the pseudocontact radius, which we define as half the distance between the two contact lines, to monitor the dynamics of the solid−liquid contact area. The same cluster analysis method as in our previous simulations17 was applied to differentiate the liquid particles from the vapor ones and to monitor the dynamics of the contact angle, pseudocontact radius, and droplet volume. The spatial mass-density distribution was used to identify the liquid−vapor interface. As a result of the cylindrical geometry, this interface is represented by an isodensity contour line, which was used to calculate the contact angle values.17 As reported in our previous work17 and illustrated for example in Figure 3 of ref 46, the interaction between the solid and the liquid leads to a layer structure with an oscillating mass distribution in the liquid near the liquid−solid interface, especially for surfaces which lead to low contact angles. The density points within the realm of layers were excluded from the fitting procedure for the contact angle determination. Only the points beyond a distance 5σ from the surface were considered. However, the contact angle values that result from this geometric analysis reflect the interfacial structure and interactions at the origin of the thermodynamics that drives the shape of the droplet. The position of the contact line was defined as the respective point where the circle which was fitted to the droplet contour intersected with the uppermost substrate atoms. Equilibrium simulations without evaporation were performed to investigate the influence of the nanoscale structure around the contact line and the associated length scale on the contact angle. In these calculations, the droplet was initially placed on a pure A surface (ε = 0.4, equilibrium contact angle 115°). After the droplet reached equilibrium at T* = 0.83, the A-type atoms in a stripe of width 2σs along the axis of the droplet were turned into B-type ones. The droplet on the new surface was equilibrated again, and a trajectory of about 104τ was used to determine the contact angle. The process of increasing the width of the B-type stripe was repeated in a symmetrical way by adding a row of B atoms on each side of the growing stripe. The distance R between a given AB boundary and the closest contact line was then the same on both sides of the droplet. Each addition of B rows was followed by an equilibration run and a production run of 104τ to calculate the value of the contact angle at a given R. In all cases, the droplet shape remained stable during the equilibration and production calculations. At the end of the growth process, the droplet sat on a pure B surface (ε = 0.7, equilibrium contact angle 66.7°). The results of this study are shown in Figure 2b and will be commented on below. The contact line may cross boundaries between the A and B domains in two different manners during evaporation. We denote A → B as the motion of the contact line when it moves from an A stripe to the neighboring B stripe
Figure 2. (a) Area of width w between the two dashed red lines considered to parametrize the local Cassie−Baxter (LCB) equation. The A (green) and B (yellow) domains denote surfaces with different energy parameters. The solid red line represents the location of the contact line. (b) Contact angle variations depending on the distance R between the contact line and the AB boundary. The black solid squares with error bars refer to the simulation data. The red stars refer to the prediction of the LCB equation with a value of 3.5σl for w. The red dashed line is the fit to a hyperbolic tangent equation. The two green dashed horizontal lines mark the contact angles on pure A and pure B surfaces.
and B → A as the motion of the contact line when it moves in the opposite direction.
3. RESULTS AND DISCUSSION 3.1. Cassie−Baxter Equation for Nanodroplets. The Cassie−Baxter (CB) equation is often used to analyze the contact angle of macroscopic droplets on chemically heterogeneous surfaces based on the wetting properties of the components that constitute the surface. On surfaces with two components A and B, it takes the following form cos θ = fA cos θA + fB cos θB (2) where fA and f B are the area fractions of A and B domains and fA + f B = 1; θ is the equilibrium contact angle on the composite surface, whereas θA and θB are the equilibrium contact angles on pure A and B surfaces, respectively. The values of fA and f B in eq 2 are determined by the respective size of the droplet and by the typical wavelength of the heterogeneity domains.36,37 In particular, if the droplet is large enough that the contact line covers a large number of A and B domains in equal proportions, then we have fA = f B = 0.5. In contrast, if the droplet’s size is comparable to the typical size of surface heterogeneities, then fA and f B may not be equal because the fraction of A and B domains covered by the three-phase contact line may be different than the fraction given by the average (global) chemical composition of the surface.36 We adopt a contact angle equation with the analytical form of eq 2 to understand the behavior of the contact angle θ of the nanometer-sized droplets on the A−B binary surfaces depending on contact angles θA and θB taken by the nanodroplet on the corresponding pure surfaces. By doing so, we assume that Young’s contact angle equation is valid to describe systems on the nanometer scale. Several works support this assumpC
DOI: 10.1021/acs.langmuir.5b01097 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir tion.47,48 The dimensions of the systems at hand are such that both situations where the droplet size is much larger than the stripe width and comparable to it must be considered. To describe these cases, we will show that a local formulation of eq 2 is suitable. We call the corresponding equation a local Cassie−Baxter (LCB) equation. In the LCB equation, only the surface composition in the local area around the three-phase contact line is assumed to influence the contact angle. To define such a local area, a distance parameter w, as shown in Figure 2a, is used. A rectangle of width 2w centered around the contact line (solid red line in Figure 2a) is considered for the calculation of parameters fA and f B. Here, fA and f B are calculated as the fraction of the rectangle defined above occupied by A-type atoms and B-type atoms, respectively. We refer to the equilibrium calculations described at the end of section 2. The variation of the contact angle with R, which is the distance between the contact line and the closest AB boundary, is shown in Figure 2b. Negative values of R indicate that the contact line is on the A stripe. Positive values of R indicate that the contact line is on the B stripe. R = 0 is obtained when the contact line is exactly located at the AB boundary. The simulated contact angle (black squares in Figure 2b) changes from θA = 115° on the pure A surface to θB = 66.7° on the pure B surface (dashed green lines). The LCB equation defined above was used to predict the contact angle variation. w is the adjustable parameter that quantifies the local nature of the interactions that drive the contact angle. We found that the variation of θ with respect to R was best reproduced when fA and f B in the LCB were obtained from w = 3.5σl (Figure 2b). This result means that the intermolecular interactions within a distance of about 7σl around the contact line determine the value of the contact angle. A smaller value of 1.57σ was reported by Ritchie et al.38 in an equilibrium molecular dynamics study of water droplets on heterogeneous surfaces with nanoscale roughness. It can thus be suggested that the local interactions which drive the contact angle are restricted to a few atomic diameters. It has already been shown on the macroscopic scale that the contact angle is a local rather than a global wetting property.27,33,36,39 Our result suggests that this conclusion remains valid on the nanometer scale and also gives an indication of the typical length scale associated with the local character of the interaction that drives the contact angle. Note that a value w = 3.5σl indicates that atomic defects can considerably influence the contact angle value and cause contact angle hysteresis. Following this analysis, we introduce the contact band area which is defined as the band of width 7σl centered on the contact line. From a macroscopic standpoint, the contact band may be interpreted as the thickness of the contact line. If the contact line defined as a thick line sits on a A or B domain, then the contact angle is expected to be equal to θA or θB, respectively. In contrast, if the contact line is located on an A area but its thickness overlaps part of a neighboring B domain, then one expects to observe a value of θ that is different than θA, as predicted by the LCB equation. Besides, a hyperbolic tangent equation (eq 3) was also used to fit the simulation data in Figure 2b θ=
⎡ 2(R − R 0) ⎤ 1 1 (θA + θB) − (θA − θB)tanh⎢ ⎥ 2 2 ds ⎣ ⎦
changes. This parameter is taken as a measure of the distance from a given AB boundary at which the droplet no longer behaves as if it were on a pure-A-type or pure-B-type surface. A value of ds = 5σl was obtained. Therefore, large variations of the contact angle are expected as soon as the contact line is located within a distance of 2.5σl from a given AB boundary. 3.2. Evaporation Pattern and Pinning Mechanism. We discuss the evaporation of system d20 first. The time evolution of the droplet during the evaporation is shown in Figure 3a. We
Figure 3. (a) Snapshots of the droplet during evaporation. The white and blue colors of the substrate denote A-type (high contact angle) and B-type (low contact angle) stripes, respectively. (b) Timedependence of contact angle θ (black line, top panel), pseudocontact line radius RCL (black line, center panel), and distance RCB between the contact lines and selected stripe boundaries (bottom panel). L1 and L2 denote the distances between the left contact line and boundaries 1 and 2 (bottom of Figure 3a), respectively. R2, R3, and R4 refer to the distances between the right contact line and the second, third, and fourth boundaries, respectively. The vertical dashed blue lines separate the four evaporation phases. The dashed red line in the top panel shows the contact angle values as predicted by the local Cassie−Baxter equation with a value of 3.5σl for w. The two horizontal green lines in the center panel indicate the position of the closest stripe boundary from the center of mass of the droplet. The two horizontal green lines in the bottom panel indicate the range of ±3.5σl around the stripe boundary.
also provide a movie of the whole evaporation of system d20 as Supporting Information. The droplet is initially centered on less-wettable stripe A, and both contact lines sit on betterwettable stripe B as shown in the top snapshot in Figure 3a. During evaporation, the contact lines first moved toward each other gradually. At a certain time point at about 2000τ, the droplet suddenly shifted to the left (second snapshot from the top in Figure 3a). During this motion, the left contact line moved toward the neighboring A stripe and the right contact line quickly crossed a B → A boundary. At about 3100τ, the droplet was centered on more-wettable stripe B (third snapshot in Figure 3a), with both contact lines on an A domain. Then, because of further evaporation, the droplet reduced its volume further and both contact lines gradually crossed the closest AB boundaries. Finally, the droplet completely vacated both A stripes (bottom snapshot in Figure 3a). To characterize the whole evaporation process, the contact angle, the pseudocontact line radius (defined as half the distance between the contact lines), and the distance between selected contact lines and AB boundaries are shown in Figure 3b. The time dependence of the contact angle and of the
(3)
where ds is a fit parameter which denotes the contact line− domain boundary distance over which the contact angle D
DOI: 10.1021/acs.langmuir.5b01097 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
to a value close to θB. A considerable slowdown can be observed in the pseudocontact radius decay while the contact lines cross the A → B boundaries after the jump phase. This drastic slowdown may be interpreted as the contact lines being pinned at the respective AB boundaries. Thus, the third phase is denoted stick. In the stick phase, evaporation follows an evaporation mode which approaches the constant contact line (CCL) mode observed in experiments.19−21 In the CCL mode, the contact angle decreases with time while the contact radius remains constant as a result of the pinned contact line. The last phase (beyond 10 000τ) corresponds to the droplet evaporation from a pure B surface. One would in principle expect to observe evaporation through the CCA mode and the slip phase to repeat. However, the droplet has become so small that it is difficult to characterize the evaporation through the wildly fluctuating contact angle. The contact angle and contact line dynamics in the three first phases described above indicate that the evaporation follows a slip−jump−stick sequence. This sequence is similar to what was observed in the Lattice Boltzmann calculations of Kusumaatmaja and Yeomans.32 In that work, a drop having the form of a planar section comparable to a projection of a cylindrical droplet was placed on a patterned surface with two kinds of alternating stripes similar to those in the present work. The droplet’s volume was decreased, and repeated slip−jump− stick sequences were observed. Although the present MD simulations are characterized by a fast change in the contact angle during the jump phase, this change occurs in a continuous way in contrast to the instantaneous jumps found by Kusumaatmaja and Yeomans. As we shall discuss in more detail below, the contact angle during evaporation in the MD simulations is driven by the area around the contact line, i.e., what we call the contact band or contact line thickness. If one assumes that the droplet shape is driven by the position of the contact line described strictly as a one-dimensional object as in the view of Kusumaatmaja and Yeomans, the sequence of evaporation is expected to occur as follows: (i) CCA as long as the CLs are on the B stripes, (ii) discontinuous jump of the contact angle and of the contact line positions as they meet AB boundaries, and (iii) strictly pinned (immobile) contact lines as the contact angle decay toward θB. Our simulations support the idea that the contact line position alone does not suffice to understand the continuous variation of the contact angle with respect to time although it can explain the basic sequence of stick, slip, and jump events. The chemical composition in the area adjacent to the contact line should also be considered. In other words, the concept of a contact line should be complemented with an additional parameter corresponding to its thickness so that fine variations in the surface composition and topography can be captured. We also note that the slip and stick phenomena were reported in experimental works.19,43 It should be noted that the slip−jump−stick sequence does not mean that evaporation must always start from the slip phase; it could start with a stick phase if the contact lines are initially located close to AB boundaries. Such an example will be given below. In the top panel of Figure 3b, we compare the variation of the actual contact angle (black line) with the variation of the contact angle predicted by the LCB equation with w = 3.5σl (red line). Excellent agreement between both curves is observed during the whole evaporation, except during the jump phase. This result is remarkable as it shows that the contact angle in the nonequilibrium system can in fact be
pseudocontact line radius (top and center panels in Figure 3b, respectively) divides the evaporation process into four phases as indicated by the vertical dashed blue lines. The first phase, named slip, follows the constant contact angle (CCA) evaporation mode. Indeed, the contact angle remains constant while the contact lines move toward each other such that the contact area decreases with time. The same behavior was observed in our work on the evaporation of cylindrical droplets17 from smooth and homogeneous surfaces. During this phase, the contact angle (as indicated by the red dashed line) has an average value equal to 66.7°, which is the contact angle θB on the pure B substrate. We conclude that the droplet behaves as if it ignored the presence of the A stripe below its center. This can be explained from the observation that there is enough space on a d20 surface for the contact lines to be located at distances larger than 3.5σl from the closest AB boundaries if the droplet is sufficiently large. Therefore, neither the left nor the right contact band overlaps any of the nearest AB boundary. The bottom panel in Figure 3b shows that the distance between the left and right contact lines and the respective closest AB boundaries is indeed larger than 3.5σl (lines L1, L2, R3, and R4) during this initial slip phase. The droplet thus behaves as if it were evaporating from a homogeneous B surface, hence the CCA evaporation mode. It is worth noting that unlike the evaporation of spherical nanodroplets, the evaporation of cylindrical nanodroplet is not influenced by the curvature of the contact line, as discussed in detail in our previous work.17 On homogeneous surfaces, spherical droplets evaporate following a mixed mode when the droplet size is small and the CCA mode when the droplet size is large. We thus anticipate that the evaporation of spherical droplets on concentrically alternating AB domains would occur following similar behavior. The second phase of the evaporation of d20 is called jump because the pseudocontact line radius decreases quickly while the contact angle increases quickly. During this phase, the left contact line suddenly reaches a distance of
