|Home | About | Journals | Submit | Contact Us | Français|
In this paper, we report molecular kinetic analyses of water spreading on hydrophobic surfaces via molecular dynamics simulation. The hydrophobic surfaces are composed of amorphous polytetrafluoroethylene (PTFE) with a static contact angle of ~112.4° for water. On the basis of the molecular kinetic theory (MKT), the influences of both viscous damping and solid-liquid retarding were analyzed in evaluating contact line friction, which characterizes the frictional force on the contact line. The unit displacement length on PTFE was estimated to be ~0.621 nm and is ~4 times as long as the bond length of C-C backbone. The static friction coefficient was found to be ~10−3 Pa·s, which is on the same order of magnitude as the dynamic viscosity of water, and increases with the droplet size. A nondimensional number defined by the ratio of the standard deviation of wetting velocity to the characteristic wetting velocity was put forward to signify the strength of the inherent contact line fluctuation and unveil the mechanism of enhanced energy dissipation in nanoscale, whereas such effect would become insignificant in macroscale. Moreover, regarding a liquid droplet on hydrophobic or superhydrophobic surfaces, an approximate solution to the base radius development was derived by an asymptotic expansion approach.
The dynamic wetting of liquids on solid surfaces is a ubiquitous phenomenon in nature and has many applications in industry. Familiar examples associated with natural processes are capillary suction in plants and self-cleaning of lotus leaves, i.e., lotus effect; and a wide range of practically important applications including surface coating1, phase change heat transfer enhancement2, 3 and surface treatment-assisted flotation4. The early attempts on dynamic wetting research can be traced back to the study of capillary tubes about two hundred years ago. During the last two decades, the long-lasting interest in dynamic wetting has been further stimulated due to the substantial progress in micro/nano-fluidics1 and microelectromechanical systems (MEMS) processing5, 6 that enables advanced surface treatments and micro\nano-structure characterization7, 8. According to the continuum theory of fluid mechanics, dynamic wetting is a thermodynamic process governed by the capillary force, viscosity and gravity9. However, the hydrodynamic analysis of the dynamic wetting yields unphysical predictions regarding the motion of three-phase (solid, liquid and vapor) contact line. When a contact line is moving at a spreading velocity uc, the classical non-slip boundary condition renders the velocity gradient at the three-phase contact line unbounded and hence leads to a shear stress singularity10, 11. In a manner analogous to the boundary layer theory, some researchers truncated the velocity profile to a slip length Lswhere the continuum concept breaks down and the Navier-Stokes equations can then be solved by relaxing the no-slip condition in the vicinity of contact line12–14. This approach is essentially a mathematical treatment to the microscopic phenomenon within the continuum domain and doesn’t provide any physical interpretation to the origin of interfacial slip. Nonetheless, there is little consensus on the order of magnitude of L s 10, 15 and the influence of intrinsic properties of the solid surface on wetting is still not well understood. Therefore, the dynamic wetting needs to be further addressed in the sub-continuum domain, i.e., at the molecular level.
In addition to the viscous dissipation and adhesion energy that are manifest in the macroscopic processes16, the excess surface free energy can be dissipated in the form of contact line friction, which becomes more prominent in the molecular level or in nanoscale, as a result of the oscillations of contact line molecules from their equilibrium positions17. The contact line friction15, 18, 19 is an important concept in the molecular kinetic theory (MKT), a conceptualized model that is based on a series of phenomenological parameters such as the equilibrium frequency K 0 and the unit displacement length λ20. The basic idea underlying MKT is that the contact line motion is essentially a rate process controlled by the corresponding energy barriers20. The continuous and macroscopic displacement of a contact line results from the collective manner of discretized forward or backward jumps of fluid molecules within the contact zone on a solid surface. In this scheme, the spontaneous spreading of contact line stems from the tilted energy barrier that favors the unidirectional displacement of contact line molecules. Hence MKT provides an efficient approach to tackling the interfacial slip paradox in dynamic wetting.
In spite of the tremendous research efforts in studying the dynamic wetting process at the sub-continuum level, there continues to be a lack of understanding on the physical meanings of molecular kinetic parameters and the dominant factors governing the dynamic wetting process still remain elusive. de Ruijter et al.21 applied the MKT to the spontaneous spreading of liquid droplets on various solid materials and found the unit displacement length λ was around 1 nm. The MKT was also applied in processes involving dewetting22, wetting on chemically heterogeneous surfaces23, forced reactive wetting24, and even to solid-liquid-liquid systems25, 26 containing ionic liquids27. In these studies, MKT-based predictions of some macroscopic parameters yielded satisfactory agreement with the experimental values. Besides, the MKT was coupled with molecular dynamics (MD) simulations and was reported to be able to depict the dynamic behaviors of nano-droplets17, 21, 28–30. However, a recent review31 pointed out that the predictive power of the MKT is limited due to the insufficient knowledge on the physical interpretation of the molecular kinetic parameters, which are hard to be related to experimentally accessible quantities or material properties. Therefore, a meticulous analysis of the molecular kinetics in dynamic wetting process is needed to understand the physical meanings of MKT parameters and the fundamental mechanisms governing wetting dynamics.
Another topic of vast scientific interest is the contact zone formation, a multiscale phenomenon spanning from the molecular scale to the macroscopic scale. It is widely accepted that the macroscopic (apparent) contact angle is accompanied by distinct and velocity-dependent microscopic contact angles14, 32. de Gennes et al.33 proposed that the contact zone could be divided into four subdomains, i.e., the molecular region, the proximal region, the central region and the distal region. Each region falls within different length scales and exhibits distinct curvatures controlled by in situ forces. In a recent experimental study, Chen and coworkers34, 35 reported a convex nanobending of contact line in the proximal region and concluded that the microscopic contact angle in MKT actually evolves beyond the molecular region of only several molecules thick. A follow-up study by Lukyanov and Likhtman32 demonstrated that the microscale contact angle forms as a result of the nonlinear distribution of the frictional force at the solid surface and hence the Young’s equation needs to be modified in order to be applied at nanoscale.
Via MD simulation in this study, we aim to provide an exhaustive analysis on the molecular kinetics of the dynamic wetting of water especially on hydrophobic surfaces, on which only limited wetting studies have been carried out. For this purpose, the solid material was chosen to be the amorphous polytetrafluorethylene (PTFE) because PTFE is the main component of Teflon®, which is one of the most widely used hydrophobic coating materials in industry. In particular, adoption of the real chemical structure of PTFE can effectively eliminate the possible artifacts that would otherwise be introduced by artificially-assembled solid structures. From an ab initio perspective, the dynamic wetting of water on hydrophobic surfaces is scrutinized at the molecular level in this work. Besides, the size effect of nano-droplets on the molecular kinetic parameters were investigated and the inherent contact line fluctuations were examined. Finally, an asymptotic solution to the dynamic wetting on hydrophobic or superhydrophobic surfaces was put forward. Our results will help advance the understanding of dynamic wetting on hydrophobic surfaces at the molecular level and may assist active control of wetting-related applications.
Most of previous molecular dynamics studies used the artificially-assembled lattice structures as the solid phase. In general, the wettability of such virtual lattice structure was tuned by varying the depth of the Lennard-Jones potential with other parameters fixed. This approach would incur inconsistency between the lattice constant and the Lennard-Jones potential and the resultant uncertainties would inextricably interfere the recognition of the physical meanings of MKT parameters. To eliminate such effects in wetting analyses, a real solid-liquid system must be simulated and the confined layer method was used in this work to construct the smooth PTFE surface (the details are given in the Supplementary Information). However, the correctness of the force field parameters for as-formed PTFE must be carefully examined. To this end, the glass transition temperature and the variation of the specific volume of PTFE were compared with the experimental data respectively. The detailed procedures to determine the glass transition temperature T g and the cubic thermal expansion coefficient α v for rubbery and glassy states can be found in the Supplementary Information. As shown in Table 1, our MD results are in good agreement with the experimental values of PTFE. Therefore, thus-formed PTFE, which owns comparable properties with the real PTFE material (Teflon®), was adopted as the substrate surface in our study. The least square fit of the MD results is presented in Fig. 1(a). The glass transition temperature, defined as the intersection of the two straight lines representing the rubbery state and the glassy state, was found to be 116.5 °C, which is in excellent agreement with the reported glass transition temperature of 115 °C for Teflon® 36.
The isotropic properties of a surface would also tremendously affect the dynamic spreading of a liquid on it. When a droplet spreads on an anisotropic surface, such as crystalline materials, its leading edge may be inclined to advance along certain directions, giving rise to irregular contact on the surface. The anisotropic effects of crystalline polymer surfaces on dynamic wetting can be found in a previous study37. The isotropy of the amorphous PTFE surfaces in our study was scrutinized in two ways. Firstly, the lateral distributions of C and F atoms in lateral (x and y) directions were checked through calculating their number density distributions38. For an amorphous polymer, the orientations of molecular chains are homogenized and the component atoms exhibit an isotropic distribution. From Fig. 1(b), it is apparent that the distributions of C and F atoms in both x and y directions share the identical profile, evidencing the exceptional isotropy of the as-formed amorphous PTFE layer. Secondly, the droplet shape and profile during wetting were also inspected. The front view snapshots of the morphing water droplet on the amorphous PTFE surface are shown in Fig. 1(c) where the droplet exhibits an approximately symmetric meniscus shape. Obviously, thus-formed PTFE is rather smooth and homogeneous, reflecting the effectiveness of the confined layer method in smoothing out the surface irregularities. Note that there are some water molecules dispersed into the vacuum environment due to evaporation. Figure 1(d) shows the snapshots of the contact area with some irregularities at the perimeter edge due to the inherent fluctuations of the contact line, which will be discussed in detail in the following section. As shown in Fig. 1(c,d), the ideally meniscus-shaped droplet profile and the almost uniform wetting perimeter demonstrate the exceptional isotropy of the amorphous PTFE surface generated by confined layer method.
Water droplets with initial diameters d=15 nm, 20nm, 25nm and 30nm were simulated to investigate the possible effects of droplet size. Each case was simulated 5 times with randomly generated Maxwell velocity distribution at 300K and the dynamic quantities including contact angle and droplet base radius were evaluated as the averaged values. To analyze the dynamic wetting behavior of water on hydrophobic PTFE surfaces, MKT is adopted in this work. The basic idea underlying the MKT is that contact line motion is actually a rate process controlled by the corresponding energy barriers among the adsorption sites on a solid surface11, 39. As illustrated in Fig. 2(a), the macroscopic three-phase contact line motion can be viewed as the statistical results of water particles’ adsorption to and detachment from the adsorption sites20, 21, 28, 39. The adsorption and detachment of molecules are controlled by the outward/forward frequency K + and inward/backward frequency K -, respectively:
Note that λ is the unit displacement length defined in MKT and reflects the statistically averaged distance that liquid molecules can move in a single step. At equilibrium state where the dynamic contact angle θ reaches the static contact angle θ 0, both the adsorption frequency and detachment frequency are equal to a constant value K 0:
The equilibrium state frequency K 0 is primarily controlled by the total energy barrier ΔG. Intrinsically, the MKT assumes that K 0 is significantly smaller than the frequency associated with the bulk thermal motion. This assumption holds true for most cases since the attractive solid-liquid interactions may restrict the random motion of contact line molecules. Herein, we took the ΔG as the total resistance against the contact line and divided it into two main components11 as illustrated in Fig. 2(b). One part, ΔG w, is ascribed to the solid-liquid retarding as a result of the work of adhesion between the solid surface and the liquid:
The other part, ΔGv, arises from the viscous damping between liquid molecules within the contact zone and those in the bulk. Therefore the wetting equation of MKT can be written in its final form as
As reflected in equation 7, the wetting velocity uc pertains to the dynamic contact angle θ via parameters λ, ΔG v and θ 0. To determine these parameters, nonlinear regression with 95% confidence interval was used to fit our MD results to the MKT. It was once suggested by some researchers that the relatively high speed range (i.e., the early stage of wetting in this study) should be excluded in applying MKT. However, recent studies31, 40 found these suggestions are actually misguiding and the contact line friction is actually dominant Figure 3(a) confirmed these findings since the wetting data obtained from MD simulation shows excellent agreement with the MKT predictions throughout the whole range of dynamic contact angle θ. The fitted values of λ, ΔG v and θ 0 for water droplets with different diameters are shown in Fig. 3(b–d).
Figure 3(b) unveils a weak dependence of static contact angle θ 0 on the droplet size as θ 0 of water on PTFE surfaces exhibit a constant value of ~112.4°, which is in line with the experimental value. The landscape of the interatomic potentials on a solid surface can be depicted as an array of potential energy wells (or barriers) as shown in Fig. 2(a). The adsorption sites are located inside each well, where water molecules are inclined to jump into in an advance movement or jump out in a retreat movement. The unit displacement length λ is the statistical result of the distance that liquid molecules within contact zone can move in a single step under the influence of both solid-liquid retarding and viscous damping. Therefore, given the prescribed solid-liquid configuration, λ should remain constant regardless of the droplet size. The independence of λ on water droplet size is confirmed in Fig. 3(b) and λ was found to be ~0.621±0.014 nm. In this study, we suggest that λ be correlated to the backbone configuration of the PTFE chains. Hence we compared the unit displacement length λ to the backbone bond (C-C) length, which is 0.1529 nm, and found that they are on the same order of magnitude and more specifically, λ is ~4 times as long as the backbone bond length.
Figure 3(c) also shows an apparent increase of the molar free energy ΔG with the droplet size. From Fig. 2(b), the motion of contact line molecules is resisted by the viscous damping ΔG v and the solid-liquid retarding ΔG w. The solid-liquid retarding refers to the mechanical resistance of the solid surface against the contact line motion in the form of molecular friction. Indeed, ΔG w exhibits a constant value of ~10.35±0.47 kJ/mol, which can be explicitly calculated from equation 6. The viscous damping ΔG v is essentially related to the viscous dissipation within the droplet. In prior studies11, 17, 31, ΔG v is simply attributed to the viscosity as
Equation 8 infers a constant ΔG v for droplets with different diameters. Nonetheless, ΔG v in Fig. 3(c) shows a systematic positive correlation with the droplet size, indicating ΔG v is not solely related to liquid viscosity η. Here we put forward an alternative analysis to clarify the contradiction. ΔG v is regarded as the cohesion forces exerted on the contact line molecules by all other water molecules, either at the vicinity of contact line or in the bulk liquid. In this scenario, ΔG v should be a function of the strength of intermolecular potential and the fraction of contact line molecules
The meaning of f c is illustrated in Fig. 4. Pij is the strength of intermolecular potential and is supposed to be related to the bulk viscosity. Apparently, increasing the intermolecular potential Pij is equivalent to enhancing the viscosity and gives rise to a stronger viscous damping. This trend coincides with equation 8. However, a larger f c implies that the motion of a single water molecule at the contact line is damped by fewer water molecules, thus resulting in a smaller ΔG v, so we have
For the current study, Pij is related to the viscosity η and is a constant for TIP4P water model, thus ΔG v should be mainly controlled by f c. As shown in Fig. 4, the contact zone is considered to span λ inward and upward from the contact line. The number n c of water molecules within this region can be approximated by
The total volume of water droplet is scaled as
Eventually f c can be calculated from its definition as
Apparently a larger droplet has a larger droplet base, giving rise to , therefore
Thus, from the chain rule, we have
Equation 15 provides a qualitative prediction for the magnitude of viscous damping and is capable of clarifying the observed increase of viscous damping with respect to droplet diameter in Fig. 3(c). Note that in equation 5, the equilibrium frequency K 0 is an exponential function of ΔG. Thus a slight increase of ΔG v with droplet size may eventually induce a significant decrease of K 0 as shown in Fig. 3(c).
The singular phenomena in the interfacial flow are usually tackled either in the viscous regime or in the inertial regime41, 42. Recently a pioneering study heralded the existence of a distinct scheme at the onset of drop coalescence43. Similarly, the MKT implies the contact line friction, which stems from the oscillations of contact line molecules from their equilibrium positions, as a nonnegligible mechanism of energy dissipation. From Rayleigh dissipation function44, the friction coefficient ξ in the MKT framework can be calculated as:
and a static friction coefficient ξ0 can be obtained as:
On one hand, ξ 0 in equation 18 measures the dissipation rate of surface free energy in the form of contact line friction; on the other hand, it can be viewed as the static resistance that contact line experiences at the equilibrium state. Figure 3(d) shows that the static friction coefficient of water on a smooth PTFE surface is on the same order of magnitude with the dynamic viscosity η. This parameter is central to many energy analyses in wetting-related phenomena. For example, in the dropwise condensation process, the accurate prediction of droplet growth and motion entails the knowledge of the magnitude of contact line friction16. The self-cleaning of lotus leaves by virtue of roll off of condensate or rain droplets, i.e., the so-called lotus effect45, partially results from the relatively small ξ 0 on the superhydrophobic surfaces.
According to MKT, water molecules in the contact zone are assumed to move inward and outward with an identical frequency K 0 at the equilibrium state. Therefore, it can be viewed as a stochastic process in which each water molecule is forced to move or jump 2K 0 τ times in a period of τ and has equal probability to either advance or retract in each step of movement. From this perspective, the motion of water particles at the equilibrium state under the MKT framework can be deemed as a Bernoulli trial ε, as shown in Fig. 5(a) and the probability of each outcome is equal to 0.5. Hence the probability p(k) in an event where water particles move a distance of (2k−2K 0 τ)λ, k = 0, 1, 2, … 2K 0 τ, can be calculated as:
where represents the binomial coefficient. Apparently, the mathematical expectation E(ε) of this Bernoulli trial is equal to 0, which coincides with the fact that the mean wetting velocity is 0 at the equilibrium state. The standard deviation σ from the mathematical expectation is given by
Here σ characterizes the deviations of water molecules from their equilibrium locations in a period of τ and is subsequently transformed into the standard deviation of wetting velocity, σ u, as:
For a certain system, τ represents the characteristic time in the wetting process and can be calculated as the characteristic length d divided by the characteristic wetting velocity , which is derived from equation 17, so we have
Equation 22 implies a linear relationship between the time constant τ and the droplet size d. To characterize the strength of the inherent fluctuations, σu needs to be nondimensionalized and a nondimensional number Z is obtained by combining equations 21 and 22:
For a certain type of liquid on a given solid surface, the unit displacement length λ has been confirmed to be constant. Thus Z is supposed to be only affected by the droplet size and the fluctuation of wetting velocity becomes less significant for a larger droplet. More specifically, λ is 0.621 nm for the current case, so we have Fig. 5(b) indicates the contact line fluctuations can only be neglected when the characteristic length exceeds micron scale. In a recent study, Wang et al.18 proved that the oscillation of microscopic contact line suggests a mechanism of enhanced energy dissipation. Subsequently, the prominent contact line fluctuation in the submicron scale is indicative of intensified dissipation of surface free energy in the nanoscale. Such effects need to be accounted for in wetting studies and imply an extra term F in modifying the Young’s equation at nanoscale:
This equation was also suggested by a more recent study32 regarding the contact angle at nanoscale, which ascribes the local contact angle variation at nanoscale to the fluctuations of microscopic force. More specifically, the addition of F in calculating contact angle at nanoscale would give rise to the convex nanobending of contact line34 for large-scale droplets.
The base radius R, defined as the radius of a circle with the same area as the (irregular) contact area of a water droplet, is also investigated in this work. In previous studies10, 39, a scaling law of was put forward to describe the development of base radius in spontaneous wetting process. The scaling approach is an effective method especially when rigorous analytical solution is unachievable for complex non-equilibrium phenomena. Nonetheless, the derivation of such a scaling law (Supplementary Information) is based on the assumption of a small contact angle, i.e., on a hydrophilic surface. However, the PTFE surface in this study is hydrophobic and the static contact angle θ 0 is larger than 110°. We thus resorted to an alternative approach to derive a heuristic scaling law of contact line spreading on hydrophobic surfaces. In the dynamic wetting process, the droplet is assumed to follow a meniscus shape, which is given by:
Equations 17 and 25 can be expanded as the Taylor series at θ = π due to the obtuse contact angle throughout the wetting process and the simplicity of the Taylor series of the trigonometric functions at θ = π. Therefore, equations 17 and 25 can be rewritten as
where a and τ are the coefficients obtained by fitting the MD results and τ can be taken as the characteristic time in the wetting process. Figure 6(a) shows the comparison of our MD data with equation 28. The proposed R ~ tanh(t) relationship demonstrates excellent agreement with the MD results. Also as shown in Fig. 6(b), the linear relationship between the time constant τ and droplet size d is in accord with equation 22.
The contact angle in the studied wetting process ranges from 160° to 110° with a divergence from the truncation angle (θ = 180°) in the Taylor expansion, from which the slight deviation of equation 28 from the MD data may origin. However, it can be expected that the R~ tanh(t) relationship would yield a more accurate prediction of the base radius evolution on a superhydrophobic surface with θ > 150°.
In this study, the dynamic wetting of water on hydrophobic surfaces is explored using large-scale molecular dynamics simulations. The dynamic wetting process is elucidated in the molecular kinetic framework and the molecular kinetic parameters are interpreted from an ab initio perspective. Our results indicate that the contact line friction originates from both the solid-liquid retarding and viscous damping. The solid-liquid retarding arises from the work of adhesion, while the viscous damping cannot be simply related to the viscosity. The contact line friction coefficient of water on smooth PTFE surfaces is on the same order of magnitude as the dynamic viscosity of water, indicating the contact line friction cannot be neglected in the analysis of processes involving moving contact line, such as condensate droplet growth and droplet coalescence. Based on these observations, we formulated a nondimensional number to characterize the inherent contact line fluctuations that are responsible for the nanobending of contact line. Finally, an asymptotic solution was derived to describe the moving contact line on a hydrophobic surface. It is believed that these findings can aid wetting-related studies and can be applied to various engineering applications.
All the simulations in this study were carried out on the single-precision molecular dynamics package Gromacs 5.1.246. The total potential energy between two atoms i and j separated by r ij is the sum of Lennard-Jones potential and the Coulombic pairwise potential:
where εij is the depth of the potential well, σ is the zero-crossing distance, ke is the Coulomb constant and q is the charge of each atom. The time step in each simulation was chosen to be 1 femtosecond (fs). The short range forces were cut-off at 1.2 nm with the Verlet scheme. The particle mesh Ewald method was adopted to handle the long range Coulombic interactions. LINear Constraint Solver (LINCS) algorithm46 were incorporated to enforce all the bonding constraints applied on each atom. A constant temperature of 300 K was maintained in both NPT (constant number N, pressure P and temperature T) and NVT (constant number N, volume V and temperature T) ensembles through a velocity-rescaling scheme.
To construct the amorphous PTFE surface, a single chain molecule consisting of 30 monomers was firstly created, which gives the exact chemical formula of CF3-(CF2)28-CF3. The force field parameters are based on the Optimized Potentials for Liquid Simulations-all atom (OPLSAA) force filed47, 48, which has shown excellent capability in simulating dynamic wetting polyethylene and polyvinyl chloride49. The detailed preparation of PTFE surfaces can be found in the Supplementary Information.
The initial structures of water droplets with diameters of d = 15 nm, 20 nm, 25 nm and 30 nm were prepared in the open source software Packmol50. The amount of the water molecules that are encapsulated into the prescribed spherical space should be consistent with the water density of 998.2 kg/m3, for example, 270930 water molecules for d = 25 nm. The TIP4P51 water model was used to account for the topology information of water molecules, since it shows great agreement with the liquid water properties52. The water sphere was then placed 0.2 nm away from the confined PTFE surface to allow interactions between the solid and the liquid atoms to take effect. The dynamic wetting process was then simulated with temperature maintained at 300 K through an NVT ensemble. The measurement of the instantaneous contact angle followed the algorithm reported by previous MD studies53, 54.
This work is financially supported by NSF CBET under grant number 1550299 and Institute for Critical Technology and Applied Science (ICTAS) at Virginia Tech. The authors also acknowledge Advanced Research Computing at Virginia Tech (http://www.arc.vt.edu) for providing computational resources and technical support that have contributed to the results reported within this paper.
J.C. and L.Z. designed the project. L.Z. carried out the simulation and data analysis. J.C. and L.Z. prepared and reviewed the manuscript.
The authors declare that they have no competing interests.
Electronic supplementary material
Supplementary information accompanies this paper at doi:10.1038/s41598-017-11350-6
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.