Slip behavior at ionic solid–fluid interfaces
Tapan G. Desai
Advanced Cooling Technologies, Inc., Lancaster, PA 17601, USA
Received 24 September 2010
In final form 1 November 2010
Available online 3 November 2010
Molecular dynamics simulation results showing that the slip behavior at the interface between an ionic solid and fluid depends on the crystal face exposed to the flowing fluid are presented. The boundary condition of fluid NaCl confined between charged octopolar (1 1 1) surfaces is ‘stick’, but, between neutral (1 0 0) NaCl surfaces is ‘slip’. Also direction dependent change from stick to slip boundary condition is observed for atomically flat, neutral dipolar (and atomically rough, neutral octopolar) (1 1 0) NaCl surfaces. Thus, by changing the surface orientation and/or its roughness, the boundary conditions at ionic solid–fluid interface can be changed from slip to stick.
© 2010 Elsevier B.V. All rights reserved.
In the simplest case of the pressure-driven laminar flow of a Newtonian fluid, the flow field can be described in terms of streamlines obtained from the Navier–Stokes equation. For conventional no-slip or stick boundary condition, the streamline velocity is assumed to be zero at the solid–fluid interface. However, recent experiments [1,2] and simulations [3,4] show that a slip can occur (especially when the fluid does not wet the solid), and the no-slip boundary condition is merely a good approximation at the macroscopic length scale [5,6]. With recent advances in miniaturization of flow devices from the microscale to nano-scale, the ‘slip effect’ becomes more important.
Molecular dynamics (MD) simulations have been widely used to show the correlations between the wetting behavior/surface roughness and the slip velocity of the liquid flow at the solid suface. In this Letter, MD simulations were used to investigate whether the surface orientation and the strength of the charges on the ionic crystal surface (NaCl solid) induce a change in the flow properties of the contacting electrolyte (NaCl fluid). The artificial introduction of charges (for charged surfaces) was avoided by using full atom description of ionic materials. Fluid NaCl confined between charged octopolar (1 1 1) surfaces exhibit stick boundary condition, but between neutral (1 0 0) NaCl surfaces a slip length of 10.9σ was found, where σ is the average ionic radius. However, the same fluid confined between atomically flat, neutral dipolar (1 1 0) NaCl surfaces exhibit stick behavior in the <1 0 0> direction, and a slip length of 4.7σ in the <1 1 0> direction. On comparison of these dipolar (1 1 0) surfaces to atomically rough, neutral octopolar (1 1 0) NaCl surfaces, a similar stick behavior was found in the <1 0 0> direction, parallel to the atomic corrugations, whereas, in the <1 1 0> direction, perpendicular to the atomic corrugations, the slip length is reduced to 1.8σ. Thus, by changing the surface orientation and/or its roughness, the boundary conditions at the ionic solid–fluid interface can be changed from slip to stick.
2. Simulation method
In this work, both solid and fluid NaCl were modeled by the empirical potentials of the Born–Mayer-Huggins (BMH) form developed by Fumi and Tosi . The trajectories of the ions were obtained by a fifth order predictor–corrector algorithm with a MD time step of t = 1.15 10-15 s. The summation method of Wolf et al.  was used to evaluate Coulomb energy, forces and stresses, i.e., direct summation of a damped and ‘charge-neutralized’ 1/r potential, with spherical cut off. This method is much faster than the widely used Ewald  or fast-multipole  methods, yet produces similar and physically more transparent results for bulk systems . In the three-dimensional Ewald summation technique for calculations of long-range Coulombic forces for slab geometry sytems that are periodic in two dimensions and have a finite length in the third dimension, Yeh and Berkowitz  proposed a new method that adds a correction term to the standard Ewald summation formula. It is known that when performing a spherically truntion formula. It is known that when performing a spherically truncated pairwise r-1 sum in a crystal or liquid, wherever one truncates, the system summed over is practically never neutral. In the 1/r summation method, Wolf et al. has shown that the Coulomb potential in an arbitrarily disordered, condensed ionic system is short ranged and, for the slabs, the local charge neutrality is achieved by viewing an ionic crystal as a molecular system consisting of Bravais lattice sites on which complete molecules are placed, with the proviso that molecules may not be broken up so as to preserve charge neutrality.
The spherical truncation cutoff, Rc = 4.28σ and damping parameter, α = 1.0 was used. The simulation model closely follows Ref. . It consists of a NaCl solid slab, which has linear dimensions of 24 24 12σ3, where σ is the average ionic radius of NaCl, σ, σCl + σNa)/2 = 2.76 Å. In this study, length is presented in the units of σ. The free surfaces terminating the crystalline slab are chosen to be either the charged (1 1 1) or the neutral (1 0 0) or (1 1 0) surfaces, and the z-axis is in the direction normal to the free surface. Periodic boundary conditions are applied in all three directions. The box size in the z-direction is fixed such that the spacing between the two free surfaces of the replicating slab is 34σ. The NaCl fluid is confined between these two free surfaces. A schematic representation for the case of flowing NaCl fluid confined between NaCl (1 0 0) surfaces is shown in Figure 1a. The flow properties were analyzed by comparing the velocity profiles of NaCl fluid in contact with fully atomistic NaCl slab terminated by either the charged (1 1 1) or the neutral (1 0 0) and (1 1 0) surfaces. As described in detail elsewhere , the NaCl crystal lattice can be viewed either as a face-centered cubic (FCC) Bravais lattice with a dipolar basis molecule or a simple cubic (SC) Bravais lattice with an octopolar (NaCl)4 basis molecule. Only in the special case of the (1 0 0) surface both the FCC and SC descriptions of the Bravais lattice yield identical arrangements of the ions at the surface. For all other surface orientations different structures are obtained, with the octopolar structure always yielding the lower energy configuration because breaking up of surface octopoles generally increases the surface energy . For example, as shown in the schematic representation in Figure 1b, the atomically flat (1 1 0) surface can only be generated via the FCC Bravais lattice with the dipolar basis. In contrast with this ‘dipolar’ surface, the corresponding ‘octopolar’ (1 1 0) surface generated via the SC Bravais lattice with the octopolar basis is not atomically flat (see Figure 1c); although, it has a significantly lower energy than the flat (i.e., dipolar) surface (see Table 1) [12,13].
The difference between the dipolar and octopolar (1 1 1) surfaces is even more pronounced. In spite of being atomically flat, the dipolar (1 1 1) surface has an infinite energy due to the longrange dipole moment in the direction normal to the surface, consisting of alternating layers of cations and anions . In contrast, the octopolar surface exhibits an interesting structure involving partially filled, charged outer planes as shown in a schematic representation in Figure 1d. While the outermost plane is still charged (positively or negatively), it contains only 1/4 of the ions of a regular, fully occupied flat (1 1 1) plane. The next plane towards the bulk is 3/4 filled and oppositely charged, while the third plane from the surface is fully occupied. Since the octopoles of which the surface is made of are dipole-moment free, this octopolar reconstruction involving the top three surface planes eliminates the long-range dipole moment of the flat, but polar (FCC-based) surface, the result being a finite energy (see Table 1). It also results in formation of alternate pyramid-like structures on the surface, which provide a certain degree of corrugation (or roughness). Figure 1d shows the atomic arrangement in the x viewing direction and Figure 1e in the y viewing direction that represents the alternate pyramid-like structures. Since periodic boundary conditions are applied in all directions, one surface is positively charged while the other is negatively charged. Though various deviations from this simple octopolar ‘ground-state’ structure have been proposed both experimentally [15,16] and by simulations , in this Letter we only consider the octopolar structure in its simplest form described above.
All simulations presented in this work were performed at temperatures well above the melting point of NaCl. To prevent melting, the ions constituting the solid slab are held fixed at the crystalline lattice parameter of the relaxed T = 0 K surface at all temperatures. All the fluid–fluid, solid–solid and solid–fluid interactions between the ions were calculated, the equation of motions for the fixed ions in the slab were not solved. It was determined that fixing the positions of ions in the solid slab does not create any artifacts .
A Poiseuille flow was simulated by applying gravitational force, ƒ, to each fluid ion in the x or y direction (parallel to the interface). The value of ƒ was chosen as ¼th of the average absolute force, acting on an ion due to short-range and coulombic interactions at T = 1600 K and a reduced number density, ρ* = 0.6 in the bulk liquid NaCl. The flow was characterized by calculating the velocity profile, which was computed by binning the velocity of the ions in the z direction (normal to the interface), and averaging the velocity of each bin over 15 x 106 steps. Data from the initial 1 x 106 steps was discarded, because the hydrodynamics requiretime to reach a stationary state in which the average velocity fluctuates around a constant value. The temperature was controlled bya velocity-rescaling thermostat, which was applied only in the directions orthogonal to the flow, in order to avoid viscous heating within the fluid film . The absence of finite size effects was confirmed by the close agreement between the velocity profiles obtained from these simulations and simulations with system size equal to ¼ of the surface area of the slab.
Flow simulations of NaCl fluid at ρ* = 0.6 and T = 1600 K confined between neutral or charged NaCl surfaces, resulted in stick behavior at the interface. Using MD simulations, Koplik et al. have shown that the slip length increases with decreasing density of the fluid . Hence to capture slip in our simulations, the fluid flow was studied at a lower value of number density, ρ* = 0.25 and at high temperature, T = 4600 K. In this Letter, the velocity is described in σ/ps units, where ps stands for picosecond. The velocity profiles from our MD simulations were compared with the analytical solution of Navier–Stokes equation for a flow with slip boundary condition, which can be expressed as  U(z) = Um(1 + dδ/w (z2/w2)), where, Um is the maximum fluid velocity in the absence of slip, δ is the effective slip length and w is half of the total separation distance between the two surfaces. Since no-slip boundary condition could not be directly imposed on this system, the value of Um cannot be estimated. Hence, the velocity profiles, U(z) for all the cases were compared and the smallest value of the maximum in the U(z) was chosen as the value of Um. This value of Um is 5.3σ/ps when the fluid is confined between octopolar (1 1 0) or dipolar (1 1 0) surfaces and the flow is in the <1 0 0> direction.
3. Results and discussion
First the results on the flow behavior of NaCl fluid confined between the charged octopolar (1 1 1) and neutral (1 0 0) surfaces are presented. In Figure 2a, the velocity profile U(z) is plotted as a function of distance between the surfaces, where z = 17 and z = 17 indicates the position of the two surfaces. The open circles reprsent data from the MD simulations, whereas the lines represent the analytical solution. The maximum velocity in the U(z) increases from 5.8σ/ps for (1 1 1) to 8.8σ/ps for (1 0 0) surfaces. From the velocity profile, it is apparent that a ‘microscopic slip’ occurs at the surface in the case where the fluid is confined between neutral (1 0 0) surfaces. Using the above-mentioned analytical solution, the value of δ for the (1 0 0) surface can be estimated as 10.9σ and for the (1 1 1) surface d ≈1.6σ. The small value of slip length of 1.6σ is interpreted as stick behavior. The (1 0 0) surface has identical planar directions on the surface, which lead to identical slip lengths in the x and y directions. The flow between the (1 1 1) surfaces is between a positively charged surface on one side and a negatively charged surface on the other side. However, the strong attraction due to the charged surfaces possibly results in stick boundary condition and a similar flow pattern near the differently charged surfaces.
The change in slip length that characterizes the interactions between the surface and the flowing fluid is influenced by a number of factors including intrinsic affinity of the fluid to the surface and surface roughness. The affinity of the fluid is strongly influenced by the strength of the solid/liquid interaction potential will be apparent from the number density profiles. Hence, in Figure 2b, the number density profiles ρ(z) of the flowing fluid in contact with the (1 0 0) and (1 1 1) surfaces are shown, where z is the coordinate perpendicular to the two surfaces (z = 17 indicates the position of the surface and z = 0 indicates the midpoint between the two surfaces). The number-density profile is simply defined as: ρ(z) = ρ+(z) + ρ (z) . In the case of the (1 0 0) surface, ρ(z) increases monotonically with increasing distances from the solid surface to reach the value of the normalized bulk density. However, in the case of the (1 1 1) surface, the density profile near the surface oscillates with increasing distance from the solid surface due to the strong attraction with the charged surface. The fluid fills the vacant sites in the slab, followed by a pronounced peak at approximately 1σ distance away from the surface. This interesting oscillatory characteristic of the density profile comes from the alternating pyramid-like structure of the charged, octopolar surface.
Joly et al.  showed that the slip length inversely depends on the surface charge and the Bjerrum length. The simulations presented here that compares the flow over vastly dissimilar surfaces, (atomically flat, neutral (1 0 0) and atomically corrugated, charged (1 1 1)), lack the flexibility to change the net surface charge, but suggests that these surface charges can have a significant impact on the flow of fluid. However, in these cases, it is difficult to distinguish between the contributions (if any) from the atomic corrugations and the charge at the surface.
In order to understand the effect on the flowing fluid due to the arrangements of cations and anions on the solid surface (in the absence of atomic corrugations and net charge at surfaces), the velocity profile for the fluid in contact with the neutral, dipolar (1 1 0) surface, which is atomically flat, was determined. In contrast to the (1 0 0) surface, the x, <1 0 0>, and y, <1 1 0> directions on the (1 1 0) surface are quite different. Here the velocity profiles of the fluid are compared when the same gravitational force is applied in the x and y-directions. As illustrated in Figure 3a, the maximum velocity is 6.83σ/ps in the y-direction and 5.3σ/ps in the x-direction. From the analytical solution, the slip length along the y-direction is about 4.7σ, whereas in the x-direction the slip length is zero, and thus the fluid exhibits essentially stick behavior. Similar to Figure 2b, in Figure 3b, the number density profiles of the flowing fluid in these two different directions are compared. In these dynamic conditions, the density of the fluid near the surface is slightly greater when the flow is in the x-direction compared to the y-direction. This indicates that a slightly stronger attraction exists between the fluid and the surface when the fluid is flowing in the x-direction resulting in stick boundary condition. Thus by changing the direction of the flow on the dipolar (1 1 0) surface, the boundary condition changes from stick to slip. This indicates that even in the absence of the net surface charge or atomic corrugations, the flow direction can have a large effect on the slip behavior at the interface.
On charge neutral surfaces, in addition to the flow in dissimilar directions, atomic corrugations can also modify the boundary condition. As mentioned above, atomic corrugations are present only in the case of the octopolar (1 1 0) surface and not in the case of atomically flat, dipolar (1 1 0) surface. The velocity profile is shown in Figure 4 for fluid flowing in the x (parallel to the atomic corrugations) and y-directions (perpendicular to the atomic over the atomic corrugations) on the octopolar (1 1 0) surface. The velocity profile in the x-direction (i.e.<1 0 0> direction) shows stick boundary condition, which matches closely with the velocity profile for dipolar (1 1 0) surface (see Figure 3a). In the y-direction (i.e.<1 1 0> direction), the value of slip length from the analytical solution is reduced to 1.8σ, compared to 4.7atomic corrugations is represented by circles for the dipolar (1 1 0) surface. Also, the maximum velocity in U(z) is reduced from 6.83σ/ps (see Figure 3a) to 5.9σ/ps due to the presence of atomic corrugations. Thus it seems that the flow over these corrugations decrease the apparent slip length. By distorting the streamlinesof the fluid flow near the surface, a similar behavior was seen in MD simulations of hydrophilic LJ liquid flowing in surface-nano-structured channels .
Priezjev and Troian  investigated the effects of periodic roughness by modeling the bounding surface as a sinusoid of amplitude a and wavenumber k = 2π/λ and found that the slip length decreased monotonically with ka. Using this description, the value of ka is zero for the dipolar (1 1 0) in <1 1 0>(value of a is 0 due to no corrugations) and is 1.53 for the octopolar (1 1 0) in <1 1 0>. Thus, the trend, similar to Ref. , shows that the sliplength reduces from 4.7 to 1.8σ with increase in the value of ka from 0 to 1.53.
In conclusion, we have performed MD simulations to study the effect on slip behavior of the NaCl fluid confined between ionic surfaces with different terminating planes. We found that the fluid confined between charged octopolar (1 1 1) NaCl surfaces sticks, but at neutral (1 0 0) surfaces the fluid slips with a slip length of 10.9σ. The slip length in the <1 0 0> direction is reduced from 10.9σ on flat (1 0 0) surfaces to a stick behavior on flat, neutral dipolar (1 1 0) surfaces. Even though these two surfaces are flat, their surface energy is very different which may play a role in changing the boundary condition (see Table 1). On dipolar (1 1 0) surfaces, the slip length is 4.7σ for flow in the <1 1 0> direction. On these surfaces, the number density of the flowing liquid was found to be slightly lower in the <1 1 0> direction compared to the <1 0 0> direction resulting in slip behavior. The flow on the atomically rough, neutral octopolar (1 1 0) NaCl surfaces in the <1 1 0> direction, perpendicular to the corrugations, results in reduction of slip length from 4.7 to 1.8σ. In summation, the slip length at the ionic solid–liquid interface depends on surface energy, direction of the flow, net surface charge and corrugations on the surface. The simulation results presented here suggest that if it is possible to fabricate surfaces with different combinations of nano-sized charge patterns that can mimic the ionic surfaces; these ‘charge-patterned’ surfaces might be used to alter the boundary conditions of the surface.
The author acknowledges Dr. Dieter Wolf (ANL) and Dr. Paul Meakin (INL) for many useful suggestions.
 J. Baudry, E. Charlaix, A. Tonck, D. Mazuyer, Langmuir 17 (2000) 5232.
 V.S.J. Craig, C. Neto, D.R.M. Williams, Phys. Rev. Lett. 87 (2001) 054504.
 J.-L. Barrat, L. Bocquet, Phys. Rev. Lett. 82 (1999) 4671.
 C. Cottin-Bizonne, J.-L. Barrat, L. Bocquet, E. Charlaix, Nat. Mat. 2 (2003) 237.
 B.-Y. Cao, M. Chen, Z.-Y. Guo, Phys. Rev. E 74 (2006) 066311.
 T.M. Squires, S.R. Quake, Rev. Mod. Phys. 77 (2005) 977.
 M.P. Tosi, F.G. Fumi, J. Phys. Chem. Solids 25 (1964) 31.
 D. Wolf, P. Keblinski, S.R. Phillpot, J. Eggebrecht, J. Chem. Phys. 110 (1999) 8254.
 P.P. Ewald, Ann. Phys. 64 (1921) 253.
 L. Greengard, Science 265 (1994) 909. references therein.
 I. Yeh, M.L. Berkowitz, J. Chem. Phys. 111 (1999) 3155.
 T.G. Desai, J. Chem. Phys. 127 (2007) 154707.
 D. Wolf, Phys. Rev. Lett. 68 (1992) 3315.
 P.W. Tasker, T.J. Bullough, Phil. Mag. A 43 (1981) 313.
 F. Finocchi, A. Barbier, J. Jupille, C. Noguera, Phys. Rev. Lett. 92 (2004) 136101.
 C. Noguera, J. Phys.: Condens. Matter 12 (2000) R367. references therein.
 J.-L. Barrat, L. Bocquet, Faraday Discuss 112 (1999) 121.
 J. Koplik, J.R. Banavar, J.F. Willemsen, Phys. Fluids A 1 (1989) 781.
 S.C. Hendy, M. Jasperse, J. Burnell, Phys. Rev. E. 72 (2005) 016303.
 L. Joly, C. Ybert, E. Trizac, L. Bocquet, J. Chem. Phys. 72 (2006) 016303.
 N. Priezjev, S.M. Troian, J. Fluid Mech. 554 (2006) 25