Matteo
Amabili
,
Emanuele
Lisi
,
Alberto
Giacomello
* and
Carlo Massimo
Casciola
Dipartimento di Ingegneria Meccanica e Aerospaziale, Università di Roma “La Sapienza”, 00184 Rome, Italy. E-mail: alberto.giacomello@uniroma1.it; Fax: +39 064458 5200; Tel: +39 06484854
First published on 9th February 2016
In this contribution we study the wetting and nucleation of vapor bubbles on nanodecorated surfaces via free energy molecular dynamics simulations. The results shed light on the stability of superhydrophobicity in submerged surfaces with nanoscale corrugations. The re-entrant geometry of the cavities under investigation is capable of sustaining a confined vapor phase within the surface roughness (Cassie state) both for hydrophobic and hydrophilic combinations of liquid and solid. The atomistic system is of nanometric size; on this scale thermally activated events can play an important role ultimately determining the lifetime of the Cassie state. Such a superhydrophobic state can break down by full wetting of the texture at large pressures (Cassie–Wenzel transition) or by nucleating a vapor bubble at negative pressures (cavitation). Specialized rare event techniques show that several pathways for wetting and cavitation are possible, due to the complex surface geometry. The related free energy barriers are of the order of 100kBT and vary with pressure. The atomistic results are found to be in semi-quantitative accord with macroscopic capillarity theory. However, the latter is not capable of capturing the density fluctuations, which determine the destabilization of the confined liquid phase at negative pressures (liquid spinodal).
However, the confined gaseous bubbles, which constitute the so-called Cassie state, may break down via two main mechanisms5 illustrated in Fig. 1: intrusion of the liquid into the surface roughness, giving rise to the so-called Wenzel state (Cassie–Wenzel transition), or formation of a critical vapor cavity (cavitation) and its subsequent dispersion in the liquid bulk. It is therefore of paramount importance for technological applications to determine quantitatively the stability of the superhydrophobic Cassie state and the kinetics of vapor loss.
A common strategy to stabilize the superhydrophobic state consists of decreasing the characteristic size of surface roughness in order to increase the capillary forces which can counteract the liquid pressure. Present-day technologies6 are capable of producing regular nanopatterns on arbitrarily large areas, which have been exploited in the study of superhydrophobicity in nanoconfinement.7 This “shrinking” strategy is also exploited in hierarchical surfaces, which combine a nanostructure for enhancing the bubble stability and a microstructure for increasing the entrapped bubble volume.8 In principle, it would be convenient to push the miniaturization down to scales at which the macroscopic models of capillarity show their limits.5 This technological challenge however calls for a more fundamental understanding of superhydrophobicity at the nanoscale. In order to address these open issues, here we use full-atom simulations and compare the results with macroscopic capillarity theory.
The Cassie–Wenzel transition and cavitation through which superhydrophobicity breaks down are, for most of the phase diagram illustrated in Fig. 1, first order transitions, in which the stable and metastable states are separated by free-energy barriers. On the nanoscale, such barriers can be overcome by thermal fluctuations, with the typical exponential dependence of the average time between transitions on the free energy barrier9 (see eqn (1) below). For barriers larger than the thermal energy kBT, this time becomes very long, giving rise to rapid but very infrequent transitions, the so-called rare events. Rare events are challenging both from the experimental and from the simulative point of view, because they require sampling of two very different timescales. In the present theoretical work, dedicated atomistic simulation techniques – the restrained molecular dynamics (RMD)10 – and continuum methods – the continuum rare events method (CREaM)11 – are adopted in order to tackle rare events.
The aim of the present work is therefore to investigate the Cassie–Wenzel transition and cavitation on a nanodecorated surface. The re-entrant structure with a typical size of 5 nm and the two surface chemistries here considered – hydrophilic and hydrophobic – are similar to those already explored in a recent communication.5 Such re-entrant geometries are exploited in order to realize the Cassie state even with liquids having low surface tension (omniphobicity).12–15 Microscopic insights into the phenomena are given by free-energy molecular dynamics simulations, which resolve the atomistic structure of the fluids and of the solid. The validity of the macroscopic capillarity theory – in which the liquid–vapor interface is sharp – is assessed against them. We find that the two descriptions yield surprisingly similar results; however, qualitative differences emerge close to the Wenzel state due to the inherent compressibility of the atomistic model. Multiple transition pathways are possible for the breakdown of superhydrophobicity each characterized by different kinetics; the number of possible pathways increases for hydrophilic surfaces. Furthermore quantitative differences in the free energy barriers are observed far from two-phase coexistence.
This paper is organized as follows. In the first section the continuum rare event method and restrained molecular dynamics are briefly reviewed. In the second section the results are discussed while the last section is left for conclusions.
τ ∼ exp(βΔΩ†), | (1) |
The free energy of a system can be expressed as a functional of some relevant quantity; e.g., in density functional theory it is the density ρ(x) of the system in the ordinary space x for which Ω ≡ Ω[ρ(x)]. In such a complex and high dimensional free-energy landscape computing the transition pathways and the related free-energy barriers is a daunting task. To alleviate this issue a reduced set of variables is often employed in order to characterize the transition. The number and the particular expression of these variables are dictated by the physics of the problem under investigation. In this work we assume that the free energy depends on a single variable, related to the filling of the cavity, see Fig. 2. For the atomistic case, this variable is usually referred to as a collective variable (CV), reflecting the fact that the atomistic degrees of freedom are mapped into a single macroscopic observable. In Section 2.4 we present a model potential to illustrate the approximations introduced by describing the system via a single variable.
(μ,V,T) = −PlVl − PvVv + γlvAlv + γsvAsv + γslAsl, | (2) |
Ω(μ,V,T) ≡ − Ωref = ΔPVv + γlv(Alv + AsvcosθY), | (3) |
By minimizing the grand potential functional in eqn (3) w.r.t. Vv, Alv, and Asv, one finds the (meta)stable configurations of the system. The conditions for stationarity are:
(4a) |
(4b) |
The solutions of the Laplace equation with Young or Gibbs boundary conditions correspond to the stationary points of the functional in eqn (3), which can be minima (stable or metastable states, e.g., Cassie or Wenzel), maxima, or saddle points (transition states, e.g., a critical vapor bubble nucleating from the surface). However, this procedure gives no information about the actual transition pathway(s) between two metastable states. In turn, the transition pathway selects the free-energy barrier which is relevant for determining the kinetics of the transition viaeqn (1).
The Continuum Rare Events Method (CREaM) has been introduced to investigate the states beside the metastable ones and to construct the transition pathways for the wetting or cavitation on arbitrary geometries. The details of this formulation are discussed in ref. 11 and 18 and here only briefly recapitulated. The basic idea consists of finding the constrained stationary points of eqn (3) under the condition of a fixed volume of vapor Vv, which is the progress variable used for the description of the Cassie–Wenzel transition and cavitation. This procedure is carried out using the method of the Lagrange multipliers. The constrained grand potential is defined as I = Ω − λ(Vv − v), where λ is the Lagrange multiplier and v is the chosen value for the volume of vapor. Imposing the stationarity of the constrained functional I yields:11
(5a) |
(5b) |
Vv = v, | (5c) |
The atoms of the lower solid wall are kept fixed at their initial positions forming an fcc lattice. The upper solid wall, instead, is used as a piston to control the pressure Pl of the liquid phase. This is achieved by applying to each particle belonging to the upper wall a constant force in the y direction (for details see ref. 22 and the ESI†). A velocity Verlet scheme is used for the time evolution of the upper solid wall. The temperature of the liquid is kept fixed at T = 0.8 via the Nosé–Hoover chain thermostat.23 The liquid temperature T, together with the liquid pressure Pl, sets the thermodynamic conditions of the system. The cavity dimensions are specified in Fig. 2, with the characteristic length being w ≃ 13. The system extends for 7σ in the z direction. Periodic boundary conditions are applied in the x and z directions. The position of the upper wall fluctuates along the y direction to enforce the required pressure, implying that the height of the computational box changes in time. In the simulation campaign we explore a range of positive and negative liquid pressures, in the interval −0.08 ≤ Pl ≤ 0.16.
As anticipated before, the occurrence of rare events implies that MD trajectories are trapped in the high probability regions of the phase space, the metastable states, and transitions between these regions are infrequent as compared to the timescale accessible to MD simulations.24 A convenient description of the system is given in terms of the collective variables ϕi(r), which depend on the microscopic state of the system r = {r1,…,rN}. For the system in Fig. 2b the simplest choice is a single collective variable counting the number of particles inside the T-shaped nanocavity (yellow rectangle).25 This is a natural choice which can be directly related to the volume Vv of the vapor domain employed in the continuum approach in Section 2.2. Relevant information about the process can be extracted by computing the probability that the observable ϕ(r) assumes a given value Ncav, p(Pl,T;ϕ(r) = Ncav), which depends also on the thermodynamic conditions Pl and T. From p(Pl,T;Ncav) (where the dependence on ϕ(r) is omitted) we define the Landau free energy:
Ω(Pl,T;Ncav) = −kBTlnp(Pl,T;Ncav). | (6) |
The above microscopic expression for the free energy can be compared with the macroscopic grand potential profile Ω(μ,V,T;Ncav) found via the CREaM (see below).
Here Ω(Pl,T;Ncav) is computed via the restrained MD method10 (RMD), which amounts to adding a biasing potential of the form Vbias(r) = k(ϕ(r) − Ncav)2/2 to the physical one. This harmonic-like potential, for suitable values of the spring constant k, restrains the system close to ϕ(r) = Ncav, forcing it to explore also regions of the phase space with low probability which are otherwise unaccessible to brute force simulations. Further details are found in reviews on rare event methods;16,26 in brief, via RMD it is possible to evaluate the gradient of Ω(Pl,T;Ncav) according to:
(7) |
(8) |
The MD engine used in these simulations is the open source code LAMMPS.27 The biasing force is computed using the rare events plugin PLUMED28 which can be interfaced with the LAMMPS.
Fig. 3 (a) Color plot and isolines for the Müller potential (arbitrary units; red corresponds to low values, while violet to large ones). α, β and γ are the minima of the potential. The blue and green points are the constrained minima along the line x = const. The black dashed curve is the actual transition pathway connecting the free-energy minima (taken from ref. 29). (b) Free-energy profile computed along the constrained minima (blue and green lines), and along the minimum energy pathway (dashed black line). (c) Müller potential computed along the line x = 0.7 (grey line in panel a), showing two minima and the free-energy barrier separating them. |
Consider the case in which in order to describe the transition from α to γ only one variable is known (or can be observed), say the x variable in Fig. 3; this is by construction an approximation since the Müller potential is two dimensional. In this section we will try to clarify which are the approximations introduced with this reduced set of variables. Applying the CREaM in one variable to Müller potential amounts to perform a constrained minimization on the line x = const. This procedure yields the set of minima at a fixed x, see, e.g., the blue and green points in Fig. 3a. This set of points is then used to construct a candidate transition pathway connecting the two minima, which is formed by piecewise smooth branches roughly corresponding to the bottom of the valleys of the potential. Finally, the value of the potential computed along the pathways yields the free-energy profile reported in Fig. 3b with solid lines.
A similar result would be obtained by applying an RMD procedure to the same problem, provided that the MD trajectory remains confined to a “valley”. Indeed, the main difference between the CREaM and RMD, even when the same variable is used to describe the system, is the fact that atomistic trajectories are affected by thermal motion. Thanks to thermal fluctuations of the order of kBT MD can escape from shallow minima. For the Müller potential, depending on temperature, thermal fluctuations could overcome the orthogonal barriers shown in Fig. 3c, leading to a more effective sampling of the phase space.
Fig. 3a and b show that the cusps in the free-energy profiles are a symptom of the presence of two neighboring valleys. Around the cusp, the reduced description of the phase space via the variable x is insufficient. In RMD, due to the integration (8), even smooth free-energy profiles can hide jumps between valleys (see, e.g., Fig. 6) which become evident only by considering additional observables (e.g., y for the Müller potential).
In this two dimensional example the exact transition pathway connecting the metastable states α and γ can be computed29 (black dashed line in Fig. 3a). In Fig. 3b the free energy along the actual transition pathway is computed and projected on the x axis for comparison with the CREaM approximation. As shown in Fig. 3a, the CREaM or RMD solutions are close to the exact transition pathway when the “valleys” are deep. The two descriptions differ appreciably only near the transition state, where the reduced description in terms of x apparently breaks down. However, the location of the metastable states and the barriers separating them are similar.
Although the Cassie–Wenzel transition and cavitation are intrinsically high-dimensional problems, a recent study32 has shown that the scenario of a free-energy landscape with deep valleys illustrated in the example above also applies to capillary problems similar to that in Fig. 3. In other words the approximation in terms of a single collective variable is generally viable, except in the vicinity of jumps between neighboring valleys.
Having fixed the thermodynamic and material parameters, we need to find a relationship connecting the atomistic collective variable Ncav and the corresponding order parameter v used in the CREaM. These two variables can be related using the so-called sharp kink approximation according to which the bulk properties of the liquid and vapor phases are extended up to the interface. With this approximation, the total number of particles Ncav in the control volume is simply given by Ncav = ρlVl + ρvv, where ρl and ρv are the bulk density of the liquid and of the vapor phases, respectively. Considering that the total control volume V is fixed, V = v + Vl, we can write Ncav = (ρv − ρl)v + const. Here the constant is chosen such that v coincides for the atomistic and continuum cases at ΔP = 0 (see Fig. 4). Finally, in all the figures v is normalized with the volume Vref of the T-shaped microstructure, which serves as a reference.
Fig. 4 Free-energy profiles for the hydrophobic (a) and hydrophilic (b) chemistries. The atomistic free energy is plotted with a red dash-dotted line. The free energy computed via the CREaM (eqn (5)) is the solid line, with each color representing a different configuration of the liquid–vapor interface. Such configurations are illustrated in the top and bottom strips for the hydrophobic and hydrophilic cases, respectively. The color of the rectangle enclosing each configuration corresponds to the branch of the same color in the CREaM free-energy profile. The insets show the atomistic average density field computed in RMD simulations, with the arrows indicating the corresponding v. |
In CREaM profiles, the macroscopic Wenzel state is attained by construction at v = 0 while the MD profiles have slightly different values. This discrepancy is due to the liquid density depletion near the solid wall in the nanostructure.33 This depletion layer is not taken into account in the sharp interface approximation used to relate Ncav and v. Its thickness depends on the chemistry of the surface and on the liquid pressure, with the larger values corresponding to the hydrophobic solid and to low pressures; this explains why the Wenzel state for the hydrophilic chemistry is closer to v = 0.
In the following discussion, in order to make a quantitative comparison between the microscopic and macroscopic description of the Cassie–Wenzel transition and of cavitation, we divide the free-energy profile in three regions
• v ≥ 0.7Vref which is relevant to the cavitation regime;
• 0.2Vref ≤ v ≤ 0.7Vref which includes the configurations explored in the Cassie–Wenzel transition;
• v ≤ 0.2Vref which corresponds to the Wenzel basin.
In the first region, the free-energy barriers and the critical volumes are compared for atomistic and macroscopic models. Only the hydrophilic chemistry is considered, since the critical bubbles for the hydrophobic case are too large (in the x-direction, see, e.g., the inset of Fig. 4a for Vv > 1) for the definition of the collective variable given in Fig. 2. In the second region, atomistic and macroscopic Cassie–Wenzel transition pathways are compared. Finally, in the third region, the behavior of the two models near the Wenzel state is analyzed.
The volume Vcrv of the critical bubble is reported in Fig. 5b as a function of ΔP. For the CREaM, the critical volume (solid lines) is dictated by the Laplace law, eqn (4a), which, in two-dimensions, gives Rc = γlv/ΔP with Rc the radius of curvature of the critical bubble. The critical volume Vcrv can be easily computed from Rc. There are two possible configurations for the critical bubble. For largely negative pressures the critical bubble is pinned at the outer corner of the T-shaped structure (blue line in Fig. 5b). For moderately negative pressures, the critical bubble meets the solid wall with the Young contact angle (green line in Fig. 5b). The atomistic critical volume (red symbols in Fig. 5b) is in fair agreement with the macroscopic ones at moderately negative pressures for which the critical bubble is not pinned. In conclusion, the Laplace equation predicts rather accurately Vcrv at the nanoscale and over a broad range of pressures.
Both RMD and the CREaM free-energy profiles show cusps, which are signatures of the transition between orthogonal shallow valleys as discussed in the model potential of Section 2.4. The first cusp, near the Wenzel state, is discussed in detail in the next subsection. The second cusp, at v ≃ 0.65Vref in the red RMD profile, corresponds to the transition between an interface pinned at the inner corner (right branch) and one with the liquid touching the cavity bottom (left branch, see the lower strip in Fig. 6).
As the liquid touches the bottom, the atomistic transition pathway switches between two different configurations: a symmetric one, with two vapor bubbles in the arms of the T-structure and an asymmetric one, with a vapor bubble in only one arm, see Fig. 6. In the macroscopic profile these two states correspond to different solutions of eqn (5). The occurrence of this morphological transition in RMD can be explained as a thermally activated jump between two shallow valleys, as illustrated in Fig. 3c. In other words, there is a free-energy barrier in the subspace orthogonal to that spanned by the collective variable Ncav which in RMD is overcome by thermal fluctuations. Clearly, the atomistic description cannot follow each single branch separately because the orthogonal free-energy barriers are of the order of kBT. In contrast, the macroscopic model is capable of following more branches but is unable to jump between any two of them and to evaluate the orthogonal free-energy barriers. The jump occurring in the atomistic simulations, identified by a star in Fig. 6, cannot be detected simply by looking at the free-energy profile, which is smooth (no cusps). However, the jump is easily revealed by examining the corresponding density fields along the transition pathways. The observed abrupt change of configuration corresponds to a free-energy jump between the two valleys identified by the CREaM, shown by the orange dashed arrow between the dark and light green profiles in Fig. 6. These two branches share, at least in the CREaM, the same free-energy gradient. Based on the observed agreement between the atomistic and macroscopic results, the same mean force is expected in eqn (7) irrespective of which of the two branches is visited by the dynamics. In this way the morphological transition is concealed by the thermodynamic integration used to reconstruct the RMD profile (eqn (8)). The macroscopic expression for the free energy, instead, is capable of distinguishing the (absolute, with no undetermined integration constant) free energy of the two branches revealing that the single bubble configuration is energetically favored.
A second atomistic pathway is reported in Fig. 6, which is generated by choosing the initial condition for the RMD simulations in the Wenzel basin (blue dash-dotted line). It is seen that the pathway selected by these initial conditions is different from the one started from Cassie. In particular, the system dynamics cannot explore the valley corresponding to two “symmetric” menisci because this is at higher free energy by ca. 20kBT (see the orange arrow for CREaM calculation). Thus the system is stuck in asymmetric configurations of the meniscus, which first pins at the lower corner of the re-entrant mouth and eventually detaches from it (top strip). This result confirms that the free-energy landscape is extremely complex and suggests that for the Cassie–Wenzel transition and for cavitation two different pathways can be followed. In order to confirm this insight, more sophisticated techniques should be used, such as the string method in collective variables.37
Similar arguments apply also to the system with hydrophilic chemistry. However in this case the number of alternative valleys is very large (see Fig. 4b and ESI†), making a detailed analysis like the one reported for the hydrophobic case very cumbersome.
A direct consequence of the upward concavity of the Wenzel basin is the existence of a liquid (or Wenzel) spinodal ΔPWsp which has no counterpart in the (incompressible) macroscopic model, see the insets in Fig. 7b and in Fig. 7a, respectively. This reflects the physical fact that superhydrophobicity can be restored at sufficiently low pressures. In contrast, the classical nucleation theory fails to capture this feature, predicting finite free energy barriers for all pressures.18
In order to make these observation more quantitative, we report in Fig. 7 the free-energy profiles Ω(v/Vref) (insets) and their derivative dΩ/d(v/Vref) for pressures greater (red) and less (blue) than the liquid spinodal. For the macroscopic model, it is seen that the Wenzel state is attained by construction at v = 0, while the point at which the free-energy derivative jumps from positive to negative values corresponds to the free-energy maximum (Fig. 7a). Decreasing the pressure amounts to shifting the free-energy derivative dΩ/d(v/Vref) ≡ λ by a constant, see eqn (5). Since, however, dΩ/d(v/Vref) has a vertical asymptote for v → 0, a maximum always exists in the continuum model. This maximum is a regular point for large negative pressures, corresponding to a critical bubble nucleating in the corner, while it is a cusp for moderately negative pressures. This cusp is generated by the presence of two valleys having different slopes, cf.Fig. 3.
In the atomistic case, instead, the free-energy derivative changes the sign twice for ΔP ≥ ΔPWsp. These two stationary points are the Wenzel and transition states, respectively, as shown in Fig. 7b. Upon decreasing pressures, the Wenzel state gradually shifts to larger v, see also ref. 42. When ΔP ≤ ΔPWsp, no stationary point exists and the Wenzel state becomes unstable. Because of the peculiar shape of the Wenzel basin, in the atomistic case the maximum is always attained at the cusp where the free-energy derivative has a discontinuity.
Comparison of the present results with previous work on simpler textures11,18 shows that the number of possible pathways dramatically increases with the complexity of the surface texture. For instance, rectangular grooves admit ca. 7 different single-bubble families,11 which constitute chunks of the overall pathway; since two rectangular ends are present in the re-entrant geometry considered here, the number of pathways doubles just because of the trivial replication of the bubble shapes at the left and at the right ends. The re-entrant mouth further increases the number of possible bubble configurations to at least 20 (cf.Fig. 4 in which the symmetric configurations are not shown). It is expected that three-dimensional geometries also induce an analogous increase in the number of pathways. This explosion of complexity may seem to preclude in practice the applicability of the CREaM/Surface Evolver framework to more sophisticated textures. However, we note that only the few pathways with the lowest free energy need to be considered, while the vast majority of them is extremely improbable; symmetry arguments are able to identify bubble configurations that are equiprobable and thus need to be computed only once (e.g., the left/right symmetry mentioned above). Furthermore, simple arguments can be devised a priori on the number and combination of bubbles which are energetically favored.11
The present results allowed for a detailed quantitative comparison of the atomistic and continuum models at the nanoscale. The major qualitative difference concerns the curvature of the free energy profile close to the Wenzel state. For the atomistic model this curvature is positive, accounting for density fluctuations of the confined liquid.38 The macroscopic model, instead, due to the assumption of liquid incompressibility, does not capture density fluctuations and features a negative curvature. This discrepancy is reflected in the different (non-classical) pathways to wetting and nucleation, which in turn lead to different estimates for the kinetics.41 Strictly related to the compressibility is the presence of a liquid spinodal – shifted by confinement as compared to the bulk one – which is only captured by the atomistic model. Quantitative differences emerge in the free energy barriers connected with cavitation, with the macroscopic model (classical nucleation theory) overestimating them. As expected, the largest discrepancies are found away from two-phase coexistence, where the size of the critical bubble becomes nanometric.
The complex free energy landscape connected with wetting and cavitation on structured surfaces has required particular attention in interpreting the results of the rare event atomistic simulations and continuum calculations. The approximation introduced by a reduced description of the transition in terms of a single variable (the volume of the vapor bubble) has been discussed in detail. This convenient choice, which is normally used, e.g., in classical nucleation theory, is capable of identifying most of the pathways, but fails when there is a morphological transition.32 Sufficient but not necessary symptoms of the failure of a reduced description are the presence of cusps in the free energy profile. Furthermore, thermodynamic integration, which is often used in atomistic free-energy methods, may fail to capture significant free-energy jumps. In order to overcome these limitations, a full description of the fluid density field would be required in order to capture the details of the phenomenon and the exact free energy barriers. The latter results call for a more profound understanding of the coarse grained description of a liquid.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5sm02794b |
This journal is © The Royal Society of Chemistry 2016 |