Impact of surface nanostructure on ice nucleation Xiang-Xiong Zhang, Min Chen, and Ming Fu Citation: The Journal of Chemical Physics 141, 124709 (2014); doi: 10.1063/1.4896149 View online: http://dx.doi.org/10.1063/1.4896149 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/12?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Hydrophobic hydration driven self-assembly of curcumin in water: Similarities to nucleation and growth under large metastability, and an analysis of water dynamics at heterogeneous surfaces J. Chem. Phys. 141, 18C501 (2014); 10.1063/1.4895539 Effects of surface interactions on heterogeneous ice nucleation for a monatomic water model J. Chem. Phys. 141, 084501 (2014); 10.1063/1.4892804 Freezing of sessile water droplets on surfaces with various roughness and wettability Appl. Phys. Lett. 104, 161609 (2014); 10.1063/1.4873345 On the characterization of crystallization and ice adhesion on smooth and rough surfaces using molecular dynamics Appl. Phys. Lett. 104, 021603 (2014); 10.1063/1.4862257 Controlled ice nucleation in microsized water droplet Appl. Phys. Lett. 81, 445 (2002); 10.1063/1.1492849
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
THE JOURNAL OF CHEMICAL PHYSICS 141, 124709 (2014)
Impact of surface nanostructure on ice nucleation Xiang-Xiong Zhang,1 Min Chen,1,a) and Ming Fu2 1 2
Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China GE Aviation, Cincinnati, Ohio 45215, USA
(Received 30 May 2014; accepted 9 September 2014; published online 26 September 2014) Nucleation of water on solid surface can be promoted noticeably when the lattice parameter of a surface matches well with the ice structure. However, the characteristic length of the surface lattice reported is generally less than 0.5 nm and is hardly tunable. In this paper, we show that a surface with nanoscale roughness can also remarkably promote ice nucleation if the characteristic length of the surface structure matches well with the ice crystal. A series of surfaces composed of periodic grooves with same depth but different widths are constructed in molecular dynamics simulations. Water cylinders are placed on the constructed surfaces and frozen at constant undercooling. The nucleation rates of the water cylinders are calculated in the simulation using the mean first-passage time method and then used to measure the nucleation promotion ability of the surfaces. Results suggest that the nucleation behavior of the supercooled water is significantly sensitive to the width of the groove. When the width of the groove matches well with the specific lengths of the ice crystal structure, the nucleation can be promoted remarkably. If the width does not match with the ice crystal, this kind of promotion disappears and the nucleation rate is even smaller than that on the smooth surface. Simulations also indicate that even when water molecules are adsorbed onto the surface structure in high-humidity environment, the solid surface can provide promising anti-icing ability as long as the characteristic length of the surface structure is carefully designed to avoid geometric match. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896149] I. INTRODUCTION
The ability to control heterogeneous nucleation with auxiliaries is of fundamental importance and has clear technological applications, such as in the development of icephobic surfaces. For this purpose, both theoretical1–5 and experimental6–8 studies have been conducted to reveal the mechanisms of the nucleation of melt on foreign surfaces. In recent decades, the development of icephobic surfaces has attracted increased attention because of its vast application prospect. Owing to their high water repellence, the anti-icing performance of hydrophobic surfaces with roughness ranging from nanometer to micrometer has been widely examined.9–17 However, given that too many factors (e.g., surface morphology,10, 18 adsorbed impurity ions,19 and local electric field in cracks20 ) can strongly affect ice nucleation, no clear relationship between the icephobicity and hydrophobicity of the substrate has been established yet. Among these factors, investigation of the role of the surface morphology has general and fundamental importance not only for ice nucleation but also for growth of diamond film,21 crystallization of proteins,22 and so forth. The present study aims to provide some new insights on the influence of surface morphology on ice nucleation. For surfaces with atomic-scale roughness, one generally accepted mechanism is that the substrate can remarkably promote ice nucleation when its lattice structure matches a) Author to whom correspondence should be addressed. Electronic mail:
[email protected]
0021-9606/2014/141(12)/124709/7/$30.00
well with the ice crystal.23–26 Free energy barrier calculations reveal that this concept also holds for crystallization of simple Lennard-Jones (LJ) fluid.27 An increase in the dimension of the surface roughness up to nano/micrometer scale results in a more complex nucleation phenomena on solid surfaces. A type of representative structure at this scale is the patterned-grooves on a substrate.8, 28 The dicing technique28 and a number of lithographic methods, such as electron-beam lithography29 and imprint lithography,30 are available to fabricate regular, periodic arrays of grooves on flat substrates, whereas irregular patterned grooves can be created by simply scratching the substrate with nanometer/micrometer-sized diamond powder.8 The grooves obtained are often rectangular in shape. If the width of the groove is large enough and the nucleus growing from one corner of the groove cannot be influenced by the wall on the other side, a theory was developed to predict that the nucleation will be promoted relative to the flat surface.2 The extent of the promotion is determined by the wettability of the substrate and the wedge angle of the groove. This prediction is also supported by molecular simulations.4, 5, 31 Simulation5 and experimental tests6, 32 show that the nucleation promotion is most effective when the wedge angle is consistent with some intrinsic angle of the crystal. If the groove is not wide enough, such as nanoscale cracks on acid surface,20 the nucleus would be confined between two sides of the groove. The behavior of ice nucleation in this kind of surface nanostructures is less understood. In this work, to determine the possibility of manipulating the surface morphology with ordered nanoscale grooves to control ice nucleation, we examine the nucleation ability of
141, 124709-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-2
Zhang, Chen, and Fu
supercooled water on a series of patterned-grooves as a function of groove width. To the best of our knowledge, this paper presents the first quantitative molecular-level analysis of the nucleation ability of supercooled water on groove-patterned surfaces. The results show that the nucleation rates of icing are remarkably sensitive to the width of the grooves. When the width of the groove matches well with the specific lengths of the ice crystal, the nucleation process is remarkably enhanced. By contrast, the grooves provide no such enhancement without the geometric match. II. COMPUTATIONAL METHODS
As shown in Fig. 1, we simulate the nucleation process of water cylinders on structured surfaces using molecular dynamics (MD) simulations. The model surface is made of platinum crystal (lattice constant = 3.922 Å) with (001) face contact with water. The grooves are constructed by periodically removing R rows of atoms along the y-axis direction where R = 0, 2, 4, 5, 6, 7, 8, 9, and 10, with three layers of atoms deleted in depth. The widths of the grooves and lands are identical. The coarse-grained model of water (mW) proposed by Molinero et al.33 is used to describe the interactions among water molecules. The mW model has been proved to be capable of accurately describing the thermodynamics and the structures of liquid water. Most importantly, the high computational efficiency of the mW model enables the observation of the appearance of ice embryo within computational time scale compared with some often-used water models, such as the simple point charged extended (SPC/E) model.34 The interactions between water molecules and the atoms of the model surface are described by simple LJ potential U = 4ε[(σ /r)12 − (σ /r)6 ], where r is the distance between one water molecule and one surface atom, σ = 2.8155 Å, and ε = 0.0572 kJ/mole, which is approximately 20% of the strength of the interaction between platinum atoms and SPC/E water molecules.34, 35 The interactions between surface atoms are also described by the LJ potential, with σ Pt = 2.471 Å and εPt = 66.9607 kJ/mole.35 Each surface atom is connected to its original position through a harmonic spring, U = K[r(t) − r(0)]2 , where r(0) is the initial position of a platinum atom, r(t) is the position at time t and K = 13 eV Å−3 is the spring
J. Chem. Phys. 141, 124709 (2014)
constant. Considering the computational efficiency, a water cylinder with radius about 8.0 nm is used instead of a spherical droplet. The water cylinder is melted and equilibrated at 300 K for about 2.0 ns during which water molecules spread to the surface and penetrate into the surface structure. The system is then quenched to 208 K and allowed to proceed until crystallization of the entire system. To ensure the validity of the estimation of nucleation rate in the simulation, the free energy barrier of the ice nucleation should be sufficiently higher than kT,36 where k is the Boltzmann constant and T is the nucleation temperature. To this end, the solid surface is expected to be highly hydrophobic and the water molecules should not be deeply undercooled. In this work, we set the strength of interaction between water molecules and surface atoms to be approximately 20% of the strength of the interaction between platinum atoms and SPC/E water molecules which results in a contact angle about 88o on smooth surface.34, 35 The nucleation temperature is set to 208 K, which is 6 K above the critical temperature at which the crystallization of mW water transforms from a nucleation-dominant to a growth-dominant process.37 Further increase in temperature leads to failure in observing nucleation phenomenon in all cases within acceptable simulation time. At least 50 simulations with different initial configurations are conducted for each R. The configurations are recorded every 1000 time steps for analysis, with the time step of 10 fs. The classical canonical ensemble (NVT) is used in all of the simulations, where N is the number of molecules in the system, V is the volume of the simulation domain, and T is the nucleation temperature. The temperature is controlled by the Nose-Hoover thermostat method and the motion equations are integrated by the Verlet algorithm. The dimensions of the simulation domain are Lx = 23.7281, 24.3164, 25.1008, 23.532, 23.532, 24.7086, 25.1008, 24.7086, and 23.532 nm for R = 0, 2, 4, 5, 6, 7, 8, 9, and 10, respectively; Ly = 4.1181 nm for R = 0, 2, 4, and 5, Ly = 4.7604 nm for R = 6–10, and Lz = 14.1192 nm for all systems. The sums of water molecule numbers and surface atom numbers are 19 642, 17 878, 18 067, 19 359, 19 359, 19 683, 19 791, 19 683, and 22 221, respectively, for R = 0, 2, 4, 5, 6, 7, 8, 9, and 10. Periodic conditions are applied in x-axis and y-axis directions and nonperiodic condition is employed in z-axis direction. To estimate and compare the enhancement of the ice nucleation on different surfaces, the nucleation rate J is used as the criterion. The nucleation rate is defined as the number of nucleation events happening in unit area and unit time. A successful nucleation event happens when the size of the ice embryo is larger enough than the critical nucleus and is growing steadily. In the framework of classical nucleation theory (CNT),38 the nucleation rate is connected with the free energy barrier G∗ by Eq. (1), J = A exp(−G∗ /kT ),
FIG. 1. Schematic illustration of the simulation system (R = 6). The blue box indicates the simulation domain and the inset illustrates the widths and height of the structure.
(1)
where A is the kinetic prefactor, k is the Boltzmann constant, and T is the nucleation temperature. In nucleation experiments, direct determination of G∗ is difficult. G∗ is usually extracted from the measured nucleation rate J by fitting Eq. (1). In MD simulations of nucleation events, both G∗ and J are able to be calculated
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-3
Zhang, Chen, and Fu
J. Chem. Phys. 141, 124709 (2014)
tion studies.45, 46 In this method, a local orientational parameter qlm (i) for each molecule i is defined as qlm (i) =
q˜lm (i) l m=−l
⎛
Nb (i)
q˜lm (i) = ⎝
1/2 , |q˜lm (i)|
(3)
⎞ Ylm (rij )⎠ /Nb (i),
(4)
j =1
FIG. 2. Plot of the MFPT curve for R = 6. Error bars are the standard deviation of ensemble average. The red solid line is the best fitting of Eq. (2).
with certain techniques. If the nucleation barrier is relatively small and the nucleation events can happen within acceptable simulation time, a number of methods are available to obtain the nucleation rate, including direct observation method,39 Yasukoka-Matsumoto method,40 survival probability method,41 and mean first-passage time (MFPT) method.36 Among these methods, quantities pertinent to nucleation, such as the size of critical nucleus, the induction time, and the nucleation barrier, can be obtained simultaneously by using MFPT method. On the contrary, for large nucleation barrier, special computational methods are needed to sample the phase space. To this end, many sampling methods, such as forward flux sampling (FFS) method42 and biased-umbrella sampling method,43 have been developed. Considering the given advantages, the MFPT method is employed in this work. To employ the MFPT formalism, as shown in Fig. 2, we need first collect a series of independent nucleation trajectories and monitor the time τ (n) in each trajectory at which the size of the maximum ice nucleus in the system appearing or exceeding n for the first time. n is the number of water molecules of the maximum ice nucleus in each recorded configuration. The average τ (n) is then obtained by averaging over all trajectories. Using the average τ (n), the nucleation rate can be derived by fitting Eq. (2),36 τ (n) =
1 {1 + erf [c(n − n∗ )]}, 2J S
where Ylm is the spherical harmonics and rij is the unit vector connecting water molecule i and one of its nearest neighboring molecules j. The neighboring molecules of i are those within a cut off distance of 0.45 nm. Nb (i) is the number of the nearest neighbors of molecule i. Herein, l = 6 is adopted to identify the structures. Then, if product m=6 q6m (i)q ∗ 6m (j ) between molecule i and j is greater than
m=−6
0.5, they are considered to be connected, where q∗ 6m is the complex conjugate of q6m . If one molecule has at least 11 connected neighbors, it is recognized as a solid molecule. Two connected solid molecules belong to the same crystal cluster. Solid clusters with different size exist in the system simultaneously,38 but only the largest cluster is monitored in the MFPT method. III. RESULTS AND DISCUSSION
We simulate ice nucleation on a series of groovepatterned surfaces. The calculated nucleation rate on smooth surface (R = 0) is 1.22 ± 0.14 × 1024 m−2 s−1 and the average induction time is approximately 14.287 ns. Using the nucleation rate on the smooth surface as the benchmark, the nucleation rates of icing as a function of groove width are plotted in Fig. 3. The fluctuation of the nucleation rates reflected from the error bars suggests that 50 nucleation trajectories are sufficient to derive a reasonable nucleation rate in our cases. As shown in Fig. 3, the evolution of the nucleation rate does not exhibit any regularity along with the width of the groove. The
(2)
where J is the nucleation rate, S is the effective contact sur√ face area, c = Z π and the size of the critical nucleus n∗ are fitting parameters with Z the Zeldovich factor, and erf(x) is the error function. The MFPT curve shown in Fig. 2 reaches a clear plateau, which indicates that the model surface is hydrophobic enough to ensure that the nucleation barrier is sufficiently higher than kT. In the MFPT calculation, an order parameter is needed to identify the cluster in solid phase. The order parameter should be able to distinguish molecules in liquid and solid phases effectively and should not be sensitive to the orientation of the crystal at the same time. In this study, we employ the bond order parameter proposed by Steinhardt et al.44 and developed by ten Wolde et al.43 to identify the cluster in solid state which has been applied successfully in numerous crystalliza-
FIG. 3. The nucleation rate and the nucleation barrier as a function of R. The nucleation rate J is normalized by the nucleation rate of icing on smooth surface JR = 0 and the dashed line indicates J/JR = 0 = 1. Error bars are comparable to the size of the squares in the figure. The short-dotted lines connecting the symbols are guides to the eye.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-4
Zhang, Chen, and Fu
nucleation is enhanced in some cases and weakened in the others when compared with that on the smooth surface. By means of a two-dimensional (2D) Ising model, Page and Sear4 calculated the nucleation rates along with the width of the 2D pore. Their results suggested that the nucleation process consists of a nucleation process in the pore and a subsequent nucleation process out of the pore. The nucleation rate in the pore decreased and the nucleation rate out of the pore increased as a function of the pore width and thus resulted in an “inverted-U”-shaped curve. The difference between Page and Sear’s results and ours should be attributed to the inherent structure of the ice crystal which cannot be considered in the Ising model.4, 31 In Eq. (1), the kinetic prefactor A can be expressed as ZfC0 provided the nucleation process is stationary, where Z = 0.12 ± 0.01 is the Zeldovich factor which varies a little in different simulations, C0 is the concentration of nucleation sites which is similar in different cases and f is the attachment rate of the monomers to the critical nucleus.38 The factor f can be expressed as 24Dn∗2/3 /λ2 , where D is the selfdiffusion coefficient of the water molecules, n∗ is the size of the critical nucleus, and λ is the atomic jump distance. It is reasonable to recognize λ as a constant in our liquid-solid phase transformation.45, 47 Surprisingly, the size of the critical nucleus changes a little for different R, with a value of 10 ± 2. According to the simulation results, all of the nucleation events start from the groove for R = 5, 7, and 9 and from the vicinity of the surface for R = 2 and 4. For R = 6, 8, and 10, about 90% of the nucleation events start from the groove and the rest start from the vicinity of the solid surface. Given that the influence of the solid surface on the ordering of the water molecules extends only approximately 1.0 nm from the surface,48 only the mean square displacements of the water molecules within 1.0 nm next to the surface were collected to determine the diffusion coefficient D. D was then calculated by fitting the mean square displacements of 400 ps with a linear function resulting in a value of 1.94 ± 0.18 × 10−9 m2 s−1 . As a consequence, the nucleation rate is mainly determined by the nucleation barrier. Fig. 3 shows that the energy barriers on surfaces with R = 5, 7, and 9 are much lower than the rest and the corresponding nucleation rates are higher. The results imply that the grooves with R = 5, 7, and 9 might have some kind of geometric preference for ice nucleation. To verify this conjecture, the width of the grooves is compared with some specific lengths of the ice crystal. First, the configurations of water molecules in the groove are checked. The snapshots of the water molecules in the first layer in the groove after the whole system has experienced crystallization are shown in Fig. 4. As shown in Fig. 4, the configurations of the water molecules in the grooves with different width differ significantly. At a groove width of R = 2, the water molecules are arranged into an atomic chain because of the confinement. At R = 4, the water molecules remain in amorphous state even after nucleation. At R = 5, 7, and 9, perfect ice structure forms with basal plane parallel to the surface (xy plane), and at R = 6, 8, and 10, the molecules crystallize into ice with many defects. Corresponding to the snapshots shown in Fig. 4, some
J. Chem. Phys. 141, 124709 (2014)
FIG. 4. Snapshots of the water molecules in the first layer in the groove. The water molecules are connected by hydrogen bonds identified by a distance criterion with a cut off of 3.0 Å.
specific widths of the hexagonal ice (Ih)49 parallel to the basal plane are given in Fig. 5. Second, a reliable method is needed to evaluate the accessible width of the grooves L1 for the water molecules. Owing to the repulsive force between the surface atoms and the water molecules, the accessible width L1 is smaller than the width L0 which is the distance between the center of mass of the wall atoms, as shown in Fig. 6. To determine L1, 2000 configurations of the system in liquid state are collected for each R. The distances between the water molecules in the groove and the center of mass of the wall are calculated and the minimum values of 1 and 2 in each configuration are collected. L1 is then determined using Eq. (5), L1 = L0 − (1 + 2) ,
(5)
where angle brackets mean ensemble average. As shown in Fig. 1, the same method is used to derive the accessible depth H of the groove for water molecules. L1 and H are used to
FIG. 5. Some specific widths of the hexagonal ice composed of complete or incomplete six-rings parallel to the basal plane at 208 K. A complete sixring is composed of six water molecules connected via hydrogen bonds. The lattice constants of the hexagonal ice are a = 4.5094 Å, b = 7.8105 Å, and c = 5.2070 Å. The water molecules in the insets are connected by red-colored hydrogen bonds. The numbers below the rings indicate the number of the six.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-5
Zhang, Chen, and Fu
J. Chem. Phys. 141, 124709 (2014)
FIG. 6. Groove widths derived from different methods as a function of R. The comparison between the accessible widths L1 and the specific widths of the ice crystal (L1-d) is represented by red open diamonds. d denotes the specific width of hexagonal ice which is closest to L1. The error bars for L1 and L1-d are comparable with the size of the squares in the figure. The solid lines connecting symbols are guides to the eye.
determine the effective surface area between the water cylinder and the surface, and then the nucleation rate. Kumar et al. suggested that the accessible width can be estimated by Eq. (6),50 where σW O = 2.8155 Å and σ OO = 2.3925 Å are two-body potential distance parameters between water-wall and water-water molecules, respectively. The accessible width is labeled as L1 in Eq. (6) to distinguish from L1. As shown in Fig. 6, L1 overestimates the accessible width, L1 = L0 −
σW O + σOO . 2
(6)
As illustrated in Fig. 6, the deviations between the accessible width and the size of the ice crystal are smaller in R = 5, 7, and 9 systems (within 1.0 Å) than those in R = 2, 4, 6, 8, and 10 systems (∼3.0 Å). As indicated in Figs. 4 and 6, as R = 2 and 4 the deviations are negative and perfect ice structure cannot form in the grooves. At R ≥ 5, the accessible width of the groove is larger than 5.207 Å (N = 1); thus, ice nucleation has the opportunity to start from the bottom of the grooves. Although the nucleation has the chance to originate from the groove as long as R ≥ 5, the rates of ice nucleation suggest that ice nucleation prefers to happen in R = 5, 7, and 9 systems rather than in R = 6, 8, and 10 systems. To achieve a better understanding, as shown in Fig. 7, we calculate the per-molecule potential for the water molecules belonging to the ice nucleus in the groove for R = 5, 7, and 9 systems as a function of (L1-d) which is defined in Fig. 6. The potential of each water molecule is due to its interactions with all other water molecules and surface atoms. As shown in the inset of Fig. 7, the ice nucleus in one groove and the two solid lands next to it are kept, while the rest ice nuclei and solid lands are removed. The solid atoms are connected to their initial positions through harmonic springs. The two remaining lands are moved towards/away from the ice nucleus in the groove with a step of 0.01 Å. After each move, the atoms in the remaining lands are connected to its new positions with harmonic springs and the system is relaxed for 1 ps. The per-molecule potential is then calculated. For comparison, the curves for R = 7 and 9 are moved towards the
FIG. 7. Per-molecule potential of the ice nucleus in the groove of R = 5, 7, and 9 systems. The lines illustrate the variation of the potential as a function of . The open circle, the solid circle, and the open square represent the deviations of L1-d for R = 5, 7, and 9 systems, respectively. The inset illustrates the remaining lands and ice nucleus. The ice nucleus in the groove is colored red and the rest part of the remaining ice nucleus is colored light blue.
curve for R = 5 along the longitudinal axis (by subtracting 0.2958 eV and 0.8837 eV for R = 7 and R = 9 systems, respectively) to let them have the same minimum potential without changing the horizontal values of the curves. As shown in Fig. 7, the per-molecule potential approaches the minimum at min = 0.559, 0.556, and 0.561 Å, respectively, for R = 5, 7, and 9. The values of min for R = 5, 7, and 9 systems are very close to 0.56 Å and almost identical. The ice structure in the grooves tends to be most stable at min. The deviations for R = 5, 7, and 9 systems are 0.7786, 0.6766, and 0.7007 Å, respectively, and the order of the corresponding nucleation rates is R = 5 < 9 < 7; thus, larger deviations exhibit lower nucleation rates. The nucleation rate for R = 7 is approximately one order larger than that in smooth case with R = 0. A typical nucleation process is shown in Fig. 8, in which ice nucleation starts spontaneously from the bottom of different grooves. The ice nuclei grow separately at first and then expand into one cluster afterward. The nucleation rate on R = 7 surfaces at 216 K is approximately two orders lower than that at 208 K, with a ratio of (J216K /J208K )R = 7 = 0.0781. Under homogenous conditions and with the combination of FFS and MD methods, Li et al. reported that nucleation rate of mW water decreased by approximately three orders as the temperature was increased from 208 K to 216 K.51 Ten simulations at 216 K are conducted on smooth
FIG. 8. Five successive snapshots (a)–(e) of a nucleation process in R = 7 system. Only the four largest growing ice nuclei are shown here and different clusters are illustrated with different colors.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-6
Zhang, Chen, and Fu
surfaces; however, nucleation events are not observed even after 450 ns of simulation. This finding further identifies the effectiveness of geometric match in promoting ice nucleation. As the accessible width of the groove deviates from the specific width of Ih six-ring, the effectiveness of the surface in promoting ice nucleation becomes weaker. Fig. 4 shows that in R = 6, 8, and 10 systems, ice structure grown from the bottom of the grooves have many defects. It implies that the increase of the deviation allows water molecules diffusing within a larger space and the ice nucleus with perfect structure might not be the most stable. The reduction of nucleation rates in R = 2, 4, 6, 8, and 10 systems relative to the flat surface suggests that even with the penetration of water molecules into the surface structure under highly humid environment, ice nucleation can still be delayed as long as the characteristic length of the structure is carefully designed. Notably, the geometric match effect weakens as the width increases further. At relatively large widths, the nucleus formed in one corner of the groove is unaffected by the wall in the other side; thus, the nucleation phenomenon is expected to be close to that described by the classical heterogeneous nucleation theory.2 Furthermore, the depth of the groove affects geometric preference. Simulations conducted at 208 K confirm that as the depth of the groove exceeds a certain threshold value, the properties of water molecules in the groove become similar to that of water molecules confined between two infinite planes. In this case, ice crystal can grow with the basal plane parallel to the yz plane even in R = 4 system and the geometric match effect is challenged. Under this situation, the influence of groove width on ice nucleation requires further investigation. Although there are many factors need to be considered, this study provides new insights on the nucleation behavior of ice on groove-patterned surfaces.
IV. CONCLUSION
To investigate the role of nanoscale surface roughness in ice nucleation, a series of grooves with accessible widths ranging from 0.1 nm to 2.0 nm are modeled. The nucleation ability of the surfaces is estimated by determining the nucleation rate of supercooled water on these surfaces. Results suggest that a well-matched accessible width and intrinsic length of the ice crystal remarkably promotes ice nucleation. Permolecule potential in the perfect ice crystal as a function of the deviation shows the presence of a minimum potential with a corresponding deviation of 0.56 Å. The promotion of ice nucleation on well-matched surfaces (R = 5, 7, and 9) is attributed to the reduction of energy cost in forming perfect ice crystal. As deviation increases, energy cost increases and promotion disappears. The reduction of the nucleation rate in R = 2, 4, 6, 8, and 10 systems relative to the flat surface suggests that even when the water molecules are adsorbed onto the surface structure under highly humid environment the solid surface can provide promising anti-icing ability as long as the characteristic length of the surface structure is carefully designed.
J. Chem. Phys. 141, 124709 (2014)
ACKNOWLEDGMENTS
This work is financially supported by GE Aviation and the Science Fund for Creative Research Groups of National Natural Science Foundation of China (Grant No. 51321002). The computations are carried out at the Tsinghua National Laboratory for Information Science and Technology, China. 1 D.
Turnbull, J. Chem. Phys. 18, 198 (1950). Sholl and N. Fletcher, Acta Metall. 18, 1083–1086 (1970). 3 X. Liu, J. Phys. Chem. B 105, 11550–11558 (2001). 4 A. J. Page and R. P. Sear, Phys. Rev. Lett. 97, 065701 (2006). 5 A. J. Page and R. P. Sear, J. Am. Chem. Soc. 131, 17550–17551 (2009). 6 Y. Diao, T. Harada, A. S. Myerson, T. A. Hatton, and B. L. Trout, Nat. Mater. 10, 867–871 (2011). 7 Y. Diao, A. S. Myerson, T. A. Hatton, and B. L. Trout, Langmuir 27, 5324– 5334 (2011). 8 J. Holbrough, J. Campbell, F. Meldrum, and H. Christenson, Cryst. Growth Des. 12, 750–755 (2012). 9 H. Saito, K. Takai, and G. Yamauchi, J. Soc. Mater. Sci., Japan 46, 185–189 (1997). 10 L. L. Cao, A. K. Jones, V. K. Sikka, J. Z. Wu, and D. Gao, Langmuir 25, 12444–12448 (2009). 11 S. A. Kulinich and M. Farzaneh, Appl. Surf. Sci. 255, 8153–8157 (2009). 12 S. A. Kulinich and M. Farzaneh, Langmuir 25, 8854–8856 (2009). 13 M. He, J. X. Wang, H. L. Li, X. L. Jin, J. J. Wang, B. Q. Liu, and Y. L. Song, Soft Matter 6, 2396–2399 (2010). 14 K. K. Varanasi, T. Deng, J. D. Smith, M. Hsu, and N. Bhate, Appl. Phys. Lett. 97, 234102 (2010). 15 S. Farhadi, M. Farzaneh, and S. Kulinich, Appl. Surf. Sci. 257, 6264–6269 (2011). 16 S. Jung, M. Dorrestijn, D. Raps, A. Das, C. M. Megaridis, and D. Poulikakos, Langmuir 27, 3059–3066 (2011). 17 S. A. Kulinich, S. Farhadi, K. Nose, and X. W. Du, Langmuir 27, 25–29 (2011). 18 S. Suzuki, A. Nakajima, N. Yoshida, M. Sakai, A. Hashimoto, Y. Kameshima, and K. Okada, Langmuir 23, 8674–8677 (2007). 19 D. N. Blair, B. L. Davis, and A. S. Dennis, J. Appl. Meteorol. 12, 1012– 1017 (1973). 20 M. Gavish, J. Wang, M. Eisenstein, M. Lahav, and L. Leiserowitz, Science 256, 815–818 (1992). 21 P. A. Dennig and D. A. Stevenson, Appl. Phys. Lett. 59, 1562–1564 (1991). 22 N. E. Chayen, E. Saridakis, and R. P. Sear, Proc. Natl. Acad. Sci. U.S.A. 103, 597–601 (2006). 23 M. Gavish, R. Popovitz-Biro, M. Lahav, and L. Leiserowitz, Science 250, 973–974 (1990). 24 S. Okawa, A. Saito, and T. Matsui, Int. J. Refrig. 29, 134–141 (2006). 25 S. J. Cox, S. M. Kathmann, J. A. Purton, M. J. Gillan, and A. Michaelides, Phys. Chem. Chem. Phys. 14, 7944–7949 (2012). 26 X. X. Zhang, M. Chen, and M. Fu, Appl. Surf. Sci. 313, 771–776 (2014). 27 H. Wang, H. Gould, and W. Klein, Phys. Rev. E 76, 031604 (2007). 28 Z. Yoshimitsu, A. Nakajima, T. Watanabe, and K. Hashimoto, Langmuir 18, 5818–5822 (2002). 29 A. Van Blaaderen, R. Ruel, and P. Wiltzius, Nature (London) 385, 321–324 (1997). 30 J. Zhang, A. Alsayed, K. H. Lin, S. Sanyal, F. Zhang, W. J. Pao, V. S. K. Balagurusamy, P. A. Heiney, and A. G. Yodh, Appl. Phys. Lett. 81, 3176– 3178 (2002). 31 L. O. Hedges and S. Whitelam, Soft Matter 8, 8624–8635 (2012). 32 V. Lopez-Mejias, A. S. Myerson, and B. L. Trout, Cryst. Growth Des. 13, 3835–3841 (2013). 33 V. Molinero and E. B. Moore, J. Phys. Chem. B 113, 4008–4016 (2009). 34 H. Berendsen, J. Grigera, and T. Straatsma, J. Phys. Chem. 91, 6269–6271 (1987). 35 P. M. Agrawal, B. M. Rice, and D. L. Thompson, Surf. Sci. 515, 21–35 (2002). 36 J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. 126, 134103 (2007). 37 E. B. Moore and V. Molinero, Nature (London) 479, 506–509 (2011). 38 D. Kashchiev, Nucleation: Basic Theory with Applications, 1st ed. (Butterworth-Heinemann, Boston, 2000). 39 J. Julin, I. Napari, and H. Vehkamäki, J. Chem. Phys. 126, 224517 (2007). 40 K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109, 8451–8462 (1998). 2 C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-7
Zhang, Chen, and Fu
41 G. Chkonia, J. Wolk, R. Strey, J. Wedekind, and D. Reguera, J. Chem. Phys.
130, 064505 (2009). 42 R. J. Allen, P. B. Warren, and P. R. Ten Wolde, Phys. Rev. Lett. 94, 18104 (2005). 43 P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 104, 9932–9947 (1996). 44 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784–805 (1983). 45 S. Auer and D. Frenkel, J. Chem. Phys. 120, 3015–3029 (2004).
J. Chem. Phys. 141, 124709 (2014) 46 S.
E. M. Lundrigan and I. Saika-Voivod, J. Chem. Phys. 131, 104503 (2009). 47 Y. Lu, X. Zhang, and M. Chen, J. Phys. Chem. B 117, 10241–10249 (2013). 48 X. X. Zhang, Y. J. Lü, and M. Chen, Mol. Phys. 111, 3808–3814 (2013). 49 K. Rottger, A. Endriss, J. Ihringer, S. Doyle, and W. Kuhs, Acta Crystallogr., Sect. B: Struct. Sci. 50, 644–648 (1994). 50 P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H. E. Stanley, Phys. Rev. E 72, 051503 (2005). 51 T. Li, D. Donadio, and G. Galli, Nat. Commun. 4, 1887 (2013).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
THE JOURNAL OF CHEMICAL PHYSICS 141, 124709 (2014)
Impact of surface nanostructure on ice nucleation Xiang-Xiong Zhang,1 Min Chen,1,a) and Ming Fu2 1 2
Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China GE Aviation, Cincinnati, Ohio 45215, USA
(Received 30 May 2014; accepted 9 September 2014; published online 26 September 2014) Nucleation of water on solid surface can be promoted noticeably when the lattice parameter of a surface matches well with the ice structure. However, the characteristic length of the surface lattice reported is generally less than 0.5 nm and is hardly tunable. In this paper, we show that a surface with nanoscale roughness can also remarkably promote ice nucleation if the characteristic length of the surface structure matches well with the ice crystal. A series of surfaces composed of periodic grooves with same depth but different widths are constructed in molecular dynamics simulations. Water cylinders are placed on the constructed surfaces and frozen at constant undercooling. The nucleation rates of the water cylinders are calculated in the simulation using the mean first-passage time method and then used to measure the nucleation promotion ability of the surfaces. Results suggest that the nucleation behavior of the supercooled water is significantly sensitive to the width of the groove. When the width of the groove matches well with the specific lengths of the ice crystal structure, the nucleation can be promoted remarkably. If the width does not match with the ice crystal, this kind of promotion disappears and the nucleation rate is even smaller than that on the smooth surface. Simulations also indicate that even when water molecules are adsorbed onto the surface structure in high-humidity environment, the solid surface can provide promising anti-icing ability as long as the characteristic length of the surface structure is carefully designed to avoid geometric match. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896149] I. INTRODUCTION
The ability to control heterogeneous nucleation with auxiliaries is of fundamental importance and has clear technological applications, such as in the development of icephobic surfaces. For this purpose, both theoretical1–5 and experimental6–8 studies have been conducted to reveal the mechanisms of the nucleation of melt on foreign surfaces. In recent decades, the development of icephobic surfaces has attracted increased attention because of its vast application prospect. Owing to their high water repellence, the anti-icing performance of hydrophobic surfaces with roughness ranging from nanometer to micrometer has been widely examined.9–17 However, given that too many factors (e.g., surface morphology,10, 18 adsorbed impurity ions,19 and local electric field in cracks20 ) can strongly affect ice nucleation, no clear relationship between the icephobicity and hydrophobicity of the substrate has been established yet. Among these factors, investigation of the role of the surface morphology has general and fundamental importance not only for ice nucleation but also for growth of diamond film,21 crystallization of proteins,22 and so forth. The present study aims to provide some new insights on the influence of surface morphology on ice nucleation. For surfaces with atomic-scale roughness, one generally accepted mechanism is that the substrate can remarkably promote ice nucleation when its lattice structure matches a) Author to whom correspondence should be addressed. Electronic mail:
[email protected]
0021-9606/2014/141(12)/124709/7/$30.00
well with the ice crystal.23–26 Free energy barrier calculations reveal that this concept also holds for crystallization of simple Lennard-Jones (LJ) fluid.27 An increase in the dimension of the surface roughness up to nano/micrometer scale results in a more complex nucleation phenomena on solid surfaces. A type of representative structure at this scale is the patterned-grooves on a substrate.8, 28 The dicing technique28 and a number of lithographic methods, such as electron-beam lithography29 and imprint lithography,30 are available to fabricate regular, periodic arrays of grooves on flat substrates, whereas irregular patterned grooves can be created by simply scratching the substrate with nanometer/micrometer-sized diamond powder.8 The grooves obtained are often rectangular in shape. If the width of the groove is large enough and the nucleus growing from one corner of the groove cannot be influenced by the wall on the other side, a theory was developed to predict that the nucleation will be promoted relative to the flat surface.2 The extent of the promotion is determined by the wettability of the substrate and the wedge angle of the groove. This prediction is also supported by molecular simulations.4, 5, 31 Simulation5 and experimental tests6, 32 show that the nucleation promotion is most effective when the wedge angle is consistent with some intrinsic angle of the crystal. If the groove is not wide enough, such as nanoscale cracks on acid surface,20 the nucleus would be confined between two sides of the groove. The behavior of ice nucleation in this kind of surface nanostructures is less understood. In this work, to determine the possibility of manipulating the surface morphology with ordered nanoscale grooves to control ice nucleation, we examine the nucleation ability of
141, 124709-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-2
Zhang, Chen, and Fu
supercooled water on a series of patterned-grooves as a function of groove width. To the best of our knowledge, this paper presents the first quantitative molecular-level analysis of the nucleation ability of supercooled water on groove-patterned surfaces. The results show that the nucleation rates of icing are remarkably sensitive to the width of the grooves. When the width of the groove matches well with the specific lengths of the ice crystal, the nucleation process is remarkably enhanced. By contrast, the grooves provide no such enhancement without the geometric match. II. COMPUTATIONAL METHODS
As shown in Fig. 1, we simulate the nucleation process of water cylinders on structured surfaces using molecular dynamics (MD) simulations. The model surface is made of platinum crystal (lattice constant = 3.922 Å) with (001) face contact with water. The grooves are constructed by periodically removing R rows of atoms along the y-axis direction where R = 0, 2, 4, 5, 6, 7, 8, 9, and 10, with three layers of atoms deleted in depth. The widths of the grooves and lands are identical. The coarse-grained model of water (mW) proposed by Molinero et al.33 is used to describe the interactions among water molecules. The mW model has been proved to be capable of accurately describing the thermodynamics and the structures of liquid water. Most importantly, the high computational efficiency of the mW model enables the observation of the appearance of ice embryo within computational time scale compared with some often-used water models, such as the simple point charged extended (SPC/E) model.34 The interactions between water molecules and the atoms of the model surface are described by simple LJ potential U = 4ε[(σ /r)12 − (σ /r)6 ], where r is the distance between one water molecule and one surface atom, σ = 2.8155 Å, and ε = 0.0572 kJ/mole, which is approximately 20% of the strength of the interaction between platinum atoms and SPC/E water molecules.34, 35 The interactions between surface atoms are also described by the LJ potential, with σ Pt = 2.471 Å and εPt = 66.9607 kJ/mole.35 Each surface atom is connected to its original position through a harmonic spring, U = K[r(t) − r(0)]2 , where r(0) is the initial position of a platinum atom, r(t) is the position at time t and K = 13 eV Å−3 is the spring
J. Chem. Phys. 141, 124709 (2014)
constant. Considering the computational efficiency, a water cylinder with radius about 8.0 nm is used instead of a spherical droplet. The water cylinder is melted and equilibrated at 300 K for about 2.0 ns during which water molecules spread to the surface and penetrate into the surface structure. The system is then quenched to 208 K and allowed to proceed until crystallization of the entire system. To ensure the validity of the estimation of nucleation rate in the simulation, the free energy barrier of the ice nucleation should be sufficiently higher than kT,36 where k is the Boltzmann constant and T is the nucleation temperature. To this end, the solid surface is expected to be highly hydrophobic and the water molecules should not be deeply undercooled. In this work, we set the strength of interaction between water molecules and surface atoms to be approximately 20% of the strength of the interaction between platinum atoms and SPC/E water molecules which results in a contact angle about 88o on smooth surface.34, 35 The nucleation temperature is set to 208 K, which is 6 K above the critical temperature at which the crystallization of mW water transforms from a nucleation-dominant to a growth-dominant process.37 Further increase in temperature leads to failure in observing nucleation phenomenon in all cases within acceptable simulation time. At least 50 simulations with different initial configurations are conducted for each R. The configurations are recorded every 1000 time steps for analysis, with the time step of 10 fs. The classical canonical ensemble (NVT) is used in all of the simulations, where N is the number of molecules in the system, V is the volume of the simulation domain, and T is the nucleation temperature. The temperature is controlled by the Nose-Hoover thermostat method and the motion equations are integrated by the Verlet algorithm. The dimensions of the simulation domain are Lx = 23.7281, 24.3164, 25.1008, 23.532, 23.532, 24.7086, 25.1008, 24.7086, and 23.532 nm for R = 0, 2, 4, 5, 6, 7, 8, 9, and 10, respectively; Ly = 4.1181 nm for R = 0, 2, 4, and 5, Ly = 4.7604 nm for R = 6–10, and Lz = 14.1192 nm for all systems. The sums of water molecule numbers and surface atom numbers are 19 642, 17 878, 18 067, 19 359, 19 359, 19 683, 19 791, 19 683, and 22 221, respectively, for R = 0, 2, 4, 5, 6, 7, 8, 9, and 10. Periodic conditions are applied in x-axis and y-axis directions and nonperiodic condition is employed in z-axis direction. To estimate and compare the enhancement of the ice nucleation on different surfaces, the nucleation rate J is used as the criterion. The nucleation rate is defined as the number of nucleation events happening in unit area and unit time. A successful nucleation event happens when the size of the ice embryo is larger enough than the critical nucleus and is growing steadily. In the framework of classical nucleation theory (CNT),38 the nucleation rate is connected with the free energy barrier G∗ by Eq. (1), J = A exp(−G∗ /kT ),
FIG. 1. Schematic illustration of the simulation system (R = 6). The blue box indicates the simulation domain and the inset illustrates the widths and height of the structure.
(1)
where A is the kinetic prefactor, k is the Boltzmann constant, and T is the nucleation temperature. In nucleation experiments, direct determination of G∗ is difficult. G∗ is usually extracted from the measured nucleation rate J by fitting Eq. (1). In MD simulations of nucleation events, both G∗ and J are able to be calculated
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-3
Zhang, Chen, and Fu
J. Chem. Phys. 141, 124709 (2014)
tion studies.45, 46 In this method, a local orientational parameter qlm (i) for each molecule i is defined as qlm (i) =
q˜lm (i) l m=−l
⎛
Nb (i)
q˜lm (i) = ⎝
1/2 , |q˜lm (i)|
(3)
⎞ Ylm (rij )⎠ /Nb (i),
(4)
j =1
FIG. 2. Plot of the MFPT curve for R = 6. Error bars are the standard deviation of ensemble average. The red solid line is the best fitting of Eq. (2).
with certain techniques. If the nucleation barrier is relatively small and the nucleation events can happen within acceptable simulation time, a number of methods are available to obtain the nucleation rate, including direct observation method,39 Yasukoka-Matsumoto method,40 survival probability method,41 and mean first-passage time (MFPT) method.36 Among these methods, quantities pertinent to nucleation, such as the size of critical nucleus, the induction time, and the nucleation barrier, can be obtained simultaneously by using MFPT method. On the contrary, for large nucleation barrier, special computational methods are needed to sample the phase space. To this end, many sampling methods, such as forward flux sampling (FFS) method42 and biased-umbrella sampling method,43 have been developed. Considering the given advantages, the MFPT method is employed in this work. To employ the MFPT formalism, as shown in Fig. 2, we need first collect a series of independent nucleation trajectories and monitor the time τ (n) in each trajectory at which the size of the maximum ice nucleus in the system appearing or exceeding n for the first time. n is the number of water molecules of the maximum ice nucleus in each recorded configuration. The average τ (n) is then obtained by averaging over all trajectories. Using the average τ (n), the nucleation rate can be derived by fitting Eq. (2),36 τ (n) =
1 {1 + erf [c(n − n∗ )]}, 2J S
where Ylm is the spherical harmonics and rij is the unit vector connecting water molecule i and one of its nearest neighboring molecules j. The neighboring molecules of i are those within a cut off distance of 0.45 nm. Nb (i) is the number of the nearest neighbors of molecule i. Herein, l = 6 is adopted to identify the structures. Then, if product m=6 q6m (i)q ∗ 6m (j ) between molecule i and j is greater than
m=−6
0.5, they are considered to be connected, where q∗ 6m is the complex conjugate of q6m . If one molecule has at least 11 connected neighbors, it is recognized as a solid molecule. Two connected solid molecules belong to the same crystal cluster. Solid clusters with different size exist in the system simultaneously,38 but only the largest cluster is monitored in the MFPT method. III. RESULTS AND DISCUSSION
We simulate ice nucleation on a series of groovepatterned surfaces. The calculated nucleation rate on smooth surface (R = 0) is 1.22 ± 0.14 × 1024 m−2 s−1 and the average induction time is approximately 14.287 ns. Using the nucleation rate on the smooth surface as the benchmark, the nucleation rates of icing as a function of groove width are plotted in Fig. 3. The fluctuation of the nucleation rates reflected from the error bars suggests that 50 nucleation trajectories are sufficient to derive a reasonable nucleation rate in our cases. As shown in Fig. 3, the evolution of the nucleation rate does not exhibit any regularity along with the width of the groove. The
(2)
where J is the nucleation rate, S is the effective contact sur√ face area, c = Z π and the size of the critical nucleus n∗ are fitting parameters with Z the Zeldovich factor, and erf(x) is the error function. The MFPT curve shown in Fig. 2 reaches a clear plateau, which indicates that the model surface is hydrophobic enough to ensure that the nucleation barrier is sufficiently higher than kT. In the MFPT calculation, an order parameter is needed to identify the cluster in solid phase. The order parameter should be able to distinguish molecules in liquid and solid phases effectively and should not be sensitive to the orientation of the crystal at the same time. In this study, we employ the bond order parameter proposed by Steinhardt et al.44 and developed by ten Wolde et al.43 to identify the cluster in solid state which has been applied successfully in numerous crystalliza-
FIG. 3. The nucleation rate and the nucleation barrier as a function of R. The nucleation rate J is normalized by the nucleation rate of icing on smooth surface JR = 0 and the dashed line indicates J/JR = 0 = 1. Error bars are comparable to the size of the squares in the figure. The short-dotted lines connecting the symbols are guides to the eye.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-4
Zhang, Chen, and Fu
nucleation is enhanced in some cases and weakened in the others when compared with that on the smooth surface. By means of a two-dimensional (2D) Ising model, Page and Sear4 calculated the nucleation rates along with the width of the 2D pore. Their results suggested that the nucleation process consists of a nucleation process in the pore and a subsequent nucleation process out of the pore. The nucleation rate in the pore decreased and the nucleation rate out of the pore increased as a function of the pore width and thus resulted in an “inverted-U”-shaped curve. The difference between Page and Sear’s results and ours should be attributed to the inherent structure of the ice crystal which cannot be considered in the Ising model.4, 31 In Eq. (1), the kinetic prefactor A can be expressed as ZfC0 provided the nucleation process is stationary, where Z = 0.12 ± 0.01 is the Zeldovich factor which varies a little in different simulations, C0 is the concentration of nucleation sites which is similar in different cases and f is the attachment rate of the monomers to the critical nucleus.38 The factor f can be expressed as 24Dn∗2/3 /λ2 , where D is the selfdiffusion coefficient of the water molecules, n∗ is the size of the critical nucleus, and λ is the atomic jump distance. It is reasonable to recognize λ as a constant in our liquid-solid phase transformation.45, 47 Surprisingly, the size of the critical nucleus changes a little for different R, with a value of 10 ± 2. According to the simulation results, all of the nucleation events start from the groove for R = 5, 7, and 9 and from the vicinity of the surface for R = 2 and 4. For R = 6, 8, and 10, about 90% of the nucleation events start from the groove and the rest start from the vicinity of the solid surface. Given that the influence of the solid surface on the ordering of the water molecules extends only approximately 1.0 nm from the surface,48 only the mean square displacements of the water molecules within 1.0 nm next to the surface were collected to determine the diffusion coefficient D. D was then calculated by fitting the mean square displacements of 400 ps with a linear function resulting in a value of 1.94 ± 0.18 × 10−9 m2 s−1 . As a consequence, the nucleation rate is mainly determined by the nucleation barrier. Fig. 3 shows that the energy barriers on surfaces with R = 5, 7, and 9 are much lower than the rest and the corresponding nucleation rates are higher. The results imply that the grooves with R = 5, 7, and 9 might have some kind of geometric preference for ice nucleation. To verify this conjecture, the width of the grooves is compared with some specific lengths of the ice crystal. First, the configurations of water molecules in the groove are checked. The snapshots of the water molecules in the first layer in the groove after the whole system has experienced crystallization are shown in Fig. 4. As shown in Fig. 4, the configurations of the water molecules in the grooves with different width differ significantly. At a groove width of R = 2, the water molecules are arranged into an atomic chain because of the confinement. At R = 4, the water molecules remain in amorphous state even after nucleation. At R = 5, 7, and 9, perfect ice structure forms with basal plane parallel to the surface (xy plane), and at R = 6, 8, and 10, the molecules crystallize into ice with many defects. Corresponding to the snapshots shown in Fig. 4, some
J. Chem. Phys. 141, 124709 (2014)
FIG. 4. Snapshots of the water molecules in the first layer in the groove. The water molecules are connected by hydrogen bonds identified by a distance criterion with a cut off of 3.0 Å.
specific widths of the hexagonal ice (Ih)49 parallel to the basal plane are given in Fig. 5. Second, a reliable method is needed to evaluate the accessible width of the grooves L1 for the water molecules. Owing to the repulsive force between the surface atoms and the water molecules, the accessible width L1 is smaller than the width L0 which is the distance between the center of mass of the wall atoms, as shown in Fig. 6. To determine L1, 2000 configurations of the system in liquid state are collected for each R. The distances between the water molecules in the groove and the center of mass of the wall are calculated and the minimum values of 1 and 2 in each configuration are collected. L1 is then determined using Eq. (5), L1 = L0 − (1 + 2) ,
(5)
where angle brackets mean ensemble average. As shown in Fig. 1, the same method is used to derive the accessible depth H of the groove for water molecules. L1 and H are used to
FIG. 5. Some specific widths of the hexagonal ice composed of complete or incomplete six-rings parallel to the basal plane at 208 K. A complete sixring is composed of six water molecules connected via hydrogen bonds. The lattice constants of the hexagonal ice are a = 4.5094 Å, b = 7.8105 Å, and c = 5.2070 Å. The water molecules in the insets are connected by red-colored hydrogen bonds. The numbers below the rings indicate the number of the six.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-5
Zhang, Chen, and Fu
J. Chem. Phys. 141, 124709 (2014)
FIG. 6. Groove widths derived from different methods as a function of R. The comparison between the accessible widths L1 and the specific widths of the ice crystal (L1-d) is represented by red open diamonds. d denotes the specific width of hexagonal ice which is closest to L1. The error bars for L1 and L1-d are comparable with the size of the squares in the figure. The solid lines connecting symbols are guides to the eye.
determine the effective surface area between the water cylinder and the surface, and then the nucleation rate. Kumar et al. suggested that the accessible width can be estimated by Eq. (6),50 where σW O = 2.8155 Å and σ OO = 2.3925 Å are two-body potential distance parameters between water-wall and water-water molecules, respectively. The accessible width is labeled as L1 in Eq. (6) to distinguish from L1. As shown in Fig. 6, L1 overestimates the accessible width, L1 = L0 −
σW O + σOO . 2
(6)
As illustrated in Fig. 6, the deviations between the accessible width and the size of the ice crystal are smaller in R = 5, 7, and 9 systems (within 1.0 Å) than those in R = 2, 4, 6, 8, and 10 systems (∼3.0 Å). As indicated in Figs. 4 and 6, as R = 2 and 4 the deviations are negative and perfect ice structure cannot form in the grooves. At R ≥ 5, the accessible width of the groove is larger than 5.207 Å (N = 1); thus, ice nucleation has the opportunity to start from the bottom of the grooves. Although the nucleation has the chance to originate from the groove as long as R ≥ 5, the rates of ice nucleation suggest that ice nucleation prefers to happen in R = 5, 7, and 9 systems rather than in R = 6, 8, and 10 systems. To achieve a better understanding, as shown in Fig. 7, we calculate the per-molecule potential for the water molecules belonging to the ice nucleus in the groove for R = 5, 7, and 9 systems as a function of (L1-d) which is defined in Fig. 6. The potential of each water molecule is due to its interactions with all other water molecules and surface atoms. As shown in the inset of Fig. 7, the ice nucleus in one groove and the two solid lands next to it are kept, while the rest ice nuclei and solid lands are removed. The solid atoms are connected to their initial positions through harmonic springs. The two remaining lands are moved towards/away from the ice nucleus in the groove with a step of 0.01 Å. After each move, the atoms in the remaining lands are connected to its new positions with harmonic springs and the system is relaxed for 1 ps. The per-molecule potential is then calculated. For comparison, the curves for R = 7 and 9 are moved towards the
FIG. 7. Per-molecule potential of the ice nucleus in the groove of R = 5, 7, and 9 systems. The lines illustrate the variation of the potential as a function of . The open circle, the solid circle, and the open square represent the deviations of L1-d for R = 5, 7, and 9 systems, respectively. The inset illustrates the remaining lands and ice nucleus. The ice nucleus in the groove is colored red and the rest part of the remaining ice nucleus is colored light blue.
curve for R = 5 along the longitudinal axis (by subtracting 0.2958 eV and 0.8837 eV for R = 7 and R = 9 systems, respectively) to let them have the same minimum potential without changing the horizontal values of the curves. As shown in Fig. 7, the per-molecule potential approaches the minimum at min = 0.559, 0.556, and 0.561 Å, respectively, for R = 5, 7, and 9. The values of min for R = 5, 7, and 9 systems are very close to 0.56 Å and almost identical. The ice structure in the grooves tends to be most stable at min. The deviations for R = 5, 7, and 9 systems are 0.7786, 0.6766, and 0.7007 Å, respectively, and the order of the corresponding nucleation rates is R = 5 < 9 < 7; thus, larger deviations exhibit lower nucleation rates. The nucleation rate for R = 7 is approximately one order larger than that in smooth case with R = 0. A typical nucleation process is shown in Fig. 8, in which ice nucleation starts spontaneously from the bottom of different grooves. The ice nuclei grow separately at first and then expand into one cluster afterward. The nucleation rate on R = 7 surfaces at 216 K is approximately two orders lower than that at 208 K, with a ratio of (J216K /J208K )R = 7 = 0.0781. Under homogenous conditions and with the combination of FFS and MD methods, Li et al. reported that nucleation rate of mW water decreased by approximately three orders as the temperature was increased from 208 K to 216 K.51 Ten simulations at 216 K are conducted on smooth
FIG. 8. Five successive snapshots (a)–(e) of a nucleation process in R = 7 system. Only the four largest growing ice nuclei are shown here and different clusters are illustrated with different colors.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-6
Zhang, Chen, and Fu
surfaces; however, nucleation events are not observed even after 450 ns of simulation. This finding further identifies the effectiveness of geometric match in promoting ice nucleation. As the accessible width of the groove deviates from the specific width of Ih six-ring, the effectiveness of the surface in promoting ice nucleation becomes weaker. Fig. 4 shows that in R = 6, 8, and 10 systems, ice structure grown from the bottom of the grooves have many defects. It implies that the increase of the deviation allows water molecules diffusing within a larger space and the ice nucleus with perfect structure might not be the most stable. The reduction of nucleation rates in R = 2, 4, 6, 8, and 10 systems relative to the flat surface suggests that even with the penetration of water molecules into the surface structure under highly humid environment, ice nucleation can still be delayed as long as the characteristic length of the structure is carefully designed. Notably, the geometric match effect weakens as the width increases further. At relatively large widths, the nucleus formed in one corner of the groove is unaffected by the wall in the other side; thus, the nucleation phenomenon is expected to be close to that described by the classical heterogeneous nucleation theory.2 Furthermore, the depth of the groove affects geometric preference. Simulations conducted at 208 K confirm that as the depth of the groove exceeds a certain threshold value, the properties of water molecules in the groove become similar to that of water molecules confined between two infinite planes. In this case, ice crystal can grow with the basal plane parallel to the yz plane even in R = 4 system and the geometric match effect is challenged. Under this situation, the influence of groove width on ice nucleation requires further investigation. Although there are many factors need to be considered, this study provides new insights on the nucleation behavior of ice on groove-patterned surfaces.
IV. CONCLUSION
To investigate the role of nanoscale surface roughness in ice nucleation, a series of grooves with accessible widths ranging from 0.1 nm to 2.0 nm are modeled. The nucleation ability of the surfaces is estimated by determining the nucleation rate of supercooled water on these surfaces. Results suggest that a well-matched accessible width and intrinsic length of the ice crystal remarkably promotes ice nucleation. Permolecule potential in the perfect ice crystal as a function of the deviation shows the presence of a minimum potential with a corresponding deviation of 0.56 Å. The promotion of ice nucleation on well-matched surfaces (R = 5, 7, and 9) is attributed to the reduction of energy cost in forming perfect ice crystal. As deviation increases, energy cost increases and promotion disappears. The reduction of the nucleation rate in R = 2, 4, 6, 8, and 10 systems relative to the flat surface suggests that even when the water molecules are adsorbed onto the surface structure under highly humid environment the solid surface can provide promising anti-icing ability as long as the characteristic length of the surface structure is carefully designed.
J. Chem. Phys. 141, 124709 (2014)
ACKNOWLEDGMENTS
This work is financially supported by GE Aviation and the Science Fund for Creative Research Groups of National Natural Science Foundation of China (Grant No. 51321002). The computations are carried out at the Tsinghua National Laboratory for Information Science and Technology, China. 1 D.
Turnbull, J. Chem. Phys. 18, 198 (1950). Sholl and N. Fletcher, Acta Metall. 18, 1083–1086 (1970). 3 X. Liu, J. Phys. Chem. B 105, 11550–11558 (2001). 4 A. J. Page and R. P. Sear, Phys. Rev. Lett. 97, 065701 (2006). 5 A. J. Page and R. P. Sear, J. Am. Chem. Soc. 131, 17550–17551 (2009). 6 Y. Diao, T. Harada, A. S. Myerson, T. A. Hatton, and B. L. Trout, Nat. Mater. 10, 867–871 (2011). 7 Y. Diao, A. S. Myerson, T. A. Hatton, and B. L. Trout, Langmuir 27, 5324– 5334 (2011). 8 J. Holbrough, J. Campbell, F. Meldrum, and H. Christenson, Cryst. Growth Des. 12, 750–755 (2012). 9 H. Saito, K. Takai, and G. Yamauchi, J. Soc. Mater. Sci., Japan 46, 185–189 (1997). 10 L. L. Cao, A. K. Jones, V. K. Sikka, J. Z. Wu, and D. Gao, Langmuir 25, 12444–12448 (2009). 11 S. A. Kulinich and M. Farzaneh, Appl. Surf. Sci. 255, 8153–8157 (2009). 12 S. A. Kulinich and M. Farzaneh, Langmuir 25, 8854–8856 (2009). 13 M. He, J. X. Wang, H. L. Li, X. L. Jin, J. J. Wang, B. Q. Liu, and Y. L. Song, Soft Matter 6, 2396–2399 (2010). 14 K. K. Varanasi, T. Deng, J. D. Smith, M. Hsu, and N. Bhate, Appl. Phys. Lett. 97, 234102 (2010). 15 S. Farhadi, M. Farzaneh, and S. Kulinich, Appl. Surf. Sci. 257, 6264–6269 (2011). 16 S. Jung, M. Dorrestijn, D. Raps, A. Das, C. M. Megaridis, and D. Poulikakos, Langmuir 27, 3059–3066 (2011). 17 S. A. Kulinich, S. Farhadi, K. Nose, and X. W. Du, Langmuir 27, 25–29 (2011). 18 S. Suzuki, A. Nakajima, N. Yoshida, M. Sakai, A. Hashimoto, Y. Kameshima, and K. Okada, Langmuir 23, 8674–8677 (2007). 19 D. N. Blair, B. L. Davis, and A. S. Dennis, J. Appl. Meteorol. 12, 1012– 1017 (1973). 20 M. Gavish, J. Wang, M. Eisenstein, M. Lahav, and L. Leiserowitz, Science 256, 815–818 (1992). 21 P. A. Dennig and D. A. Stevenson, Appl. Phys. Lett. 59, 1562–1564 (1991). 22 N. E. Chayen, E. Saridakis, and R. P. Sear, Proc. Natl. Acad. Sci. U.S.A. 103, 597–601 (2006). 23 M. Gavish, R. Popovitz-Biro, M. Lahav, and L. Leiserowitz, Science 250, 973–974 (1990). 24 S. Okawa, A. Saito, and T. Matsui, Int. J. Refrig. 29, 134–141 (2006). 25 S. J. Cox, S. M. Kathmann, J. A. Purton, M. J. Gillan, and A. Michaelides, Phys. Chem. Chem. Phys. 14, 7944–7949 (2012). 26 X. X. Zhang, M. Chen, and M. Fu, Appl. Surf. Sci. 313, 771–776 (2014). 27 H. Wang, H. Gould, and W. Klein, Phys. Rev. E 76, 031604 (2007). 28 Z. Yoshimitsu, A. Nakajima, T. Watanabe, and K. Hashimoto, Langmuir 18, 5818–5822 (2002). 29 A. Van Blaaderen, R. Ruel, and P. Wiltzius, Nature (London) 385, 321–324 (1997). 30 J. Zhang, A. Alsayed, K. H. Lin, S. Sanyal, F. Zhang, W. J. Pao, V. S. K. Balagurusamy, P. A. Heiney, and A. G. Yodh, Appl. Phys. Lett. 81, 3176– 3178 (2002). 31 L. O. Hedges and S. Whitelam, Soft Matter 8, 8624–8635 (2012). 32 V. Lopez-Mejias, A. S. Myerson, and B. L. Trout, Cryst. Growth Des. 13, 3835–3841 (2013). 33 V. Molinero and E. B. Moore, J. Phys. Chem. B 113, 4008–4016 (2009). 34 H. Berendsen, J. Grigera, and T. Straatsma, J. Phys. Chem. 91, 6269–6271 (1987). 35 P. M. Agrawal, B. M. Rice, and D. L. Thompson, Surf. Sci. 515, 21–35 (2002). 36 J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. 126, 134103 (2007). 37 E. B. Moore and V. Molinero, Nature (London) 479, 506–509 (2011). 38 D. Kashchiev, Nucleation: Basic Theory with Applications, 1st ed. (Butterworth-Heinemann, Boston, 2000). 39 J. Julin, I. Napari, and H. Vehkamäki, J. Chem. Phys. 126, 224517 (2007). 40 K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109, 8451–8462 (1998). 2 C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
124709-7
Zhang, Chen, and Fu
41 G. Chkonia, J. Wolk, R. Strey, J. Wedekind, and D. Reguera, J. Chem. Phys.
130, 064505 (2009). 42 R. J. Allen, P. B. Warren, and P. R. Ten Wolde, Phys. Rev. Lett. 94, 18104 (2005). 43 P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 104, 9932–9947 (1996). 44 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784–805 (1983). 45 S. Auer and D. Frenkel, J. Chem. Phys. 120, 3015–3029 (2004).
J. Chem. Phys. 141, 124709 (2014) 46 S.
E. M. Lundrigan and I. Saika-Voivod, J. Chem. Phys. 131, 104503 (2009). 47 Y. Lu, X. Zhang, and M. Chen, J. Phys. Chem. B 117, 10241–10249 (2013). 48 X. X. Zhang, Y. J. Lü, and M. Chen, Mol. Phys. 111, 3808–3814 (2013). 49 K. Rottger, A. Endriss, J. Ihringer, S. Doyle, and W. Kuhs, Acta Crystallogr., Sect. B: Struct. Sci. 50, 644–648 (1994). 50 P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H. E. Stanley, Phys. Rev. E 72, 051503 (2005). 51 T. Li, D. Donadio, and G. Galli, Nat. Commun. 4, 1887 (2013).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.135.239.97 On: Thu, 18 Dec 2014 05:22:52
