Molecular Dynamics Simulations of Slip on Curved Surfaces
Simulations Moléculaires Dynamiques du glissement sur des surfaces incurvées
^{1}
Department of Chemical Engineering, Imperial College London, London
SW7 2AZ, – UK
^{2}
Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge
CB2 1EW – UK
email: daniel.ross08@imperial.ac.uk  esb30@cam.ac.uk
^{*} Corresponding author
We present Molecular Dynamics (MD) simulations of liquid water confined within nanoscale geometries, including slitlike and cylindrical graphitic pores. These equilibrium results are used for calculating friction coefficients, which in turn can be used to calculate slip lengths. The slip length is a material property independent of the fluid flow rate. It is therefore a better quantity for study than the fluid velocity at the wall, also known as the slip velocity. Once the slip length has been found as a function of surface curvature, it can be used to parameterise Lattice Boltzmann (LB) simulations. These larger scale simulations are able to tell us about how fluid transport is affected by slip in complex geometries; not just limited to single pores. Applications include flow and transport in nanoporous engine valve deposits and gas shales. The friction coefficient is found to be a function of curvature and is higher for fluid on convex surfaces and lower for concave surfaces. Both concave and convex surfaces approach the same value of the friction coefficient, which is constant above some critical radius of curvature, here found to be 7.4 ± 2.9 nm. The constant value of the friction coefficient is 10,000 ± 600 kg m^{−2} s^{−1}, which is equivalent to a slip length of approximately 67 ± 4 nm.
Résumé
Nous vous présentons ici des simulations Moléculaires Dynamiques (MD) d’eau liquide confinée dans des géométries nanométriques, comprenant des pores graphitiques cylindriques et semblables à une fente. Ces résultats équilibrés sont utilisés pour calculer des coefficients de frottement, qui à leur tour peuvent être utilisés pour calculer des longueurs de glissement. La longueur de glissement est une propriété matérielle indépendante du débit de fluide. Elle représente donc une meilleure quantité à étudier que la vitesse de fluide contre le mur, également connue sous le nom de vitesse de glissement. Une fois que la longueur de glissement a été trouvée, en fonction de la courbure de la surface, elle peut être utilisée pour paramétrer les simulations de Lattice Boltzmann (LB). Ces simulations sur une échelle plus grande nous indiquent comment le transport de fluide est affecté par le glissement dans des géométries complexes, pas seulement limitées à des pores uniques. Les applications comprennent l’écoulement et le transport dans des dépôts de soupape de moteur nanoporeux et des schistes de gaz. Le coefficient de frottement dépend de la courbure et est supérieur pour les fluides passant sur des surfaces convexes et inférieur pour des surfaces concaves. Les surfaces concaves et convexes approchent de la même valeur de coefficient de frottement. Cette valeur est constante audessus d’un certain rayon de courbure critique, lequel a été défini ici comme étant de 7,4 ± 2,9 nm. La valeur constante du coefficient de frottement est de 10 000 ± 600 kg m^{–2} s^{–1}, ce qui équivaut à une longueur de glissement d’approximativement 67 ± 4 nm.
© D.A. Ross and E.S. Boek, published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Introduction
When studying nanoporous media, conventional Computational Fluid Dynamics (CFD) and LatticeBoltzmann (LB) [1–3] methods may not be suitable to predict flow and transport properties in their current conditions. The reason is that these techniques make the assumption of a noslip boundary condition at fluidsolid interfaces, which is not necessarily the case when pore sizes become comparable to molecular diameters. If slip is not taken into account, then the mean velocity in the channel can be underestimated drastically. Quantifying the amount of slip in a channel of given curvature and solidfluid combination is a challenge that requires explicit molecular scale simulation techniques such as Molecular Dynamics (MD). Previously, we have examined imbibition in nanopores [4] and wettability alteration [5] using MD simulations. Here, we extend this to study slip within graphitic slit pores and carbon nanotubes [6].
It is fairly straightforward to determine the slip length in a straight channel or cylindrical pore using MD simulations. However, in reality most porous materials are not as ideal as that, the geometry can be complex and nontrivial. Unfortunately, the slip in a complicated geometry, with many twists and turns and even roughness at the atomic scale, is not the same as that for a straight channel. It turns out that the slip length changes with the curvature of the solidfluid interface in both the direction of flow and the perpendicular direction [6].
The aim of this work is to find the effect of surface curvature on the friction coefficient, and therefore the slip length. Given that we know this relation, we can upscale our systems from the molecular scale to the continuum scale. This is implemented via a slip parameter in a nanoscale LB model for complex geometries and long time simulations. This allows for avoidance of underestimating flow rates due to assumptions of no slip. At larger scales we can easily perform simulations on large complex geometries that are inaccessible to molecular simulations. From these, parameters such as permeability can be extracted for further upscaling. This is of particular use to hydrocarbon recovery from tight pore spaces, such as gas shales and tight rocks.
In addition to this work, there have been extensive works employing MD simulations to study slip; a small sample is given in references [7–14].
1 Model
In this study, we use the open source MD simulation package: GROMACS (GROningen MAchine for Chemical Simulations) [15].
A simple approach for calculating the slip length is to apply a driving force to fluid molecules in a pore that has been filled at a realistic density. This driving force is applied via updating atom velocities each timestep to reflect some acceleration. The result, after a long equilibration period, is some overall flow through the pore. The steadystate velocity profile can then be extracted via binning and ensemble averaging over many timesteps. As an added benefit, it becomes trivial with this approach to find the viscosity by comparison of the velocity profile to Poiseuille flow. Finally, using the wall velocity, u_{wall}, and the shear rate at the wall, du/dx_{wall}, the slip length, b, can be found:(1)
However, this method has its weak points. For instance, the time required for equilibration is large, in that steadystate must be achieved before the velocity profile can be extracted. This time is dependent on the size of the driving force and can be in excess of 10 s of ns. Similarly the time over which to average results is long and, depending on the system, the slip length can vary with shear rate at these unrealistically high driving forces. Typically forces in excess of O(10) times gravity are required for computational feasibility. With lower forces comes longer equilibration times as well as times for result analysis. Ideally, each data point should be made up of 510 repeated simulations in order to ensure sufficient statistical accuracy. Given that roughly 1020 data points are necessary in order to ensure the correct functional relationship between slip and geometry has been found, up to 200 simulations would be required. This makes the above method even more impractical, as even a single simulation is time consuming, with even slightly more realistic driving forces.
Additionally, in the case of fluids which exhibit high levels of slip, it becomes difficult to extract a slip length via this method due to an virtually flat velocity profile. In this case, it is possible to work back from the NavierStokes equation, together with the slip boundary condition shown in Equation (1), to find the slip length from the mean fluid velocity. The resulting expressions are then found for slip in a channel or cylinder respectively:(2) (3)
where u_{m} is the mean fluid velocity in the direction of flow, μ is the shear viscosity, ρ is the fluid density, a is the applied acceleration, R is the cylinder radius, and H is the height of the channel. This has been used to generate the NonEquilibrium MD (NEMD) results for comparison with the equilibrium method described below. Here, a relatively low acceleration of 0.0001 nm ps^{−2} is used in order to avoid high shear rate effects [16]. The mean fluid velocity is extracted by averaging over the final 35 ns (17,500 timesteps) of a 40 ns simulation for each acceleration and geometry.
An alternative approach for quantifying and calculating the amount of slip in a given geometry, used herein, was proposed by Bocquet and Barrat [17]. The method involves calculating the friction coefficient from Equilibrium MD (EMD) simulations and then relating this back to a sliplength. The advantage of this is that equilibration times are significantly reduced due to no need for nonequilibrium simulations. This is also particularly useful as the friction coefficient can be calculated for any flow direction from the same data set. The friction coefficient, λ, is calculated using the following GreenKubo type relation:(4)
where A is the total surface area, k_{B} is the Boltzmann constant, T is temperature and F(t) is the sum of the tangential forces applied by the fluid on the solid surface. 〈F(t)F(0)〉_{equ} is the autocorrelation function averaged over equilibrium configurations. The slip length can be recovered from the friction coefficient using the following relation:(5)
In Equation (5), μ is the shear viscosity, which can be found from a singlephase fluid simulation using a similar GreenKubo type expression as used for the friction coefficient [18]:(6)
where P_{αβ} are the offdiagonal elements of the pressure tensor and V is the system volume.
The systems presented here are fully atomistic, i.e. no coarse graining has been done above the atomic scale. The OPLSAA (Optimized Potentials for Liquid Simulations  All Atom) forcefield [19] is used here.
Nonbonded interactions are handled by a truncated and shifted LennardJones potential, with a cutoff length, r_{c}, of 0.9 nm. For atomatom separations larger than r_{c}, the van der Waals potential is set to zero. To avoid a discontinuity at r_{c}, the potential is shifted up, such that the truncated potential is zero at r_{c}. The LennardJones potential, U_{LJ}, is calculated as follows:(7)
where σ_{ij} is the distance between atoms i and j at which the interparticle potential is zero and ε_{ij} is the depth of the potential well which dictates the strength of the interaction. Nonbonded pair interactions are calculated as the geometric average, as is the normal case for the OPLS forcefield:(8) (9)
Bonded interactions are handled with harmonic stretching and angle potentials. There are no dihedral bonds in this system. See OPLSAA forcefield [19] for further details on force constants and parameters used in the water model.
Initially, we start with CarbonNanoTube (CNT) pores of varying sizes, between 1 and 10 nm diameters. For the base case we also simulate flat graphitic slabs. Similarly, we have varied the separation between adjacent slabs. Each of these geometries is filled with water and can be seen in Figure 1. Future work shall involve extending these systems to include longchain hydrocarbon mixtures in place of water.
Figure 1. Top left: Water filled, single layered, armchair carbon nanotube. Top right: Single layered armchair carbon nanotube with water on the outside. Bottom left: Water confined between two symmetrical, single layered, graphitic plates. Bottom right: Graphitic structure. 
The SPC/E (Extended Simple Point Charge) water model [20] has been used here, however consistent results have been found with the TIP3P model [21]. Parameters for the SPC/E water model are listed in Table 1.
SPC/E water parameters
Solid surfaces are represented by frozen carbon atoms arranged to mimic graphitic structures. To this end, bonds are inherently confined at a length of 0.14 nm and bond angles at 120°. These atoms are uncharged with σ and ε values of 0.35500 nm and 0.29288 kJ mol^{−1} respectively. The lack of flexible bonds in this model could be a limitation of this work and the effects should be investigated in future.
The final results simulation only requires equilibrium MD data of fluid confined to these geometries. However, the equilibration method affects the number of fluid molecules present and thus the calculated friction coefficient.
The equilibration method we employ here involves multiple stages. Firstly, the pore geometries are generated in MATLAB [22] using the bond lengths and angles defined earlier. The dimensions of each geometry are roughly 5 nm in length for each CNT and 5×5 nm for the flat slabs. Empty space is added around the pores, such that the shortest distance between adjacent images is 5 nm. This is to ensure periodicity does not affect the results. To prevent fluid from being on the exterior of the pore during equilibration, barriers are added at either end of each pore. Any fluid molecules allowed outside of the pore during equilibration would interact across the solid pore boundary and influence the fluid inside the pore. The size of the simulation box is increased along one dimension to allow for fluid reservoirs, before filling the empty space with randomly placed water molecules. Fluid is then removed from the pore itself and from areas that would be inaccessible during the simulation, i.e. the region closed off by the barriers.
The system is allowed to equilibrate at 300 K for 20 ps and then a pressure of 1 bar is applied via the use of two pistons for 500 ps. All systems are in a cubic simulation box and are periodic in all directions. Each reservoir is at least 5 nm long plus enough to allow for fluid that would enter the porespace. This allows for a bulk fluid region between the piston and the pore of at least 2 nm, preventing the formation of a pseudopore.
Each piston is a single layered square lattice structure of carbon like atoms with unrealistic bond lengths, comparable to σ. Atoms are held in alignment by high constraining force constants in the harmonic potential calculations.
Once the system has equilibrated at constant pressure, system snapshots are saved every 10 ps for 1 ns. The reservoirs, barriers and pistons are removed from the resulting set of geometries and the mean fluid density is calculated. Geometries with fluid densities about the mean are then chosen for the results run. The systems are reequilibrated at 300 K and at constant volume for 200 ps before a final 1 ns results run. An example final system, before reservoirs and barriers are removed, is shown in Figure 2. The fluid filled systems after reservoir and barrier removal can be seen in Figure 1.
Figure 2. Cross section of equilibrated system of water inside a carbon nanotube, before removal of the reservoirs. Barriers are included during equilibration only. 
The results simulation is in the canonical ensemble, NVT, such that the total number of molecules (N), total volume (V) and total temperature (T) are kept constant. The system temperature is maintained by the velocityrescaling thermostat [23] although an alternative to this is the NoséHoover thermostat [24, 25]. Both of these result in the correct canonical ensemble, although velocityrescaling was found to provide slightly improved performance.
During the results run, coordinates are output every 5 timesteps (0.01 ps) and the forces are recalculated post simulation, excluding all but the interactions between fluid atoms and the solid carbon surface. The sum of these forces are autocorrelated for use in Equation (4).
2 Results
Figure 3 shows the results of friction coefficient vs radius for fluid inside and outside of a CNT. The central points are the base case for fluid confined between two slabs of Graphene, where the radius plotted is half the separation. The base case is the case of zero curvature, effectively infinite radius. As expected, the separation has no effect on the amount of slip between flat plates.
Figure 3. Friction coefficient calculated using Equation (4) for water inside and outside of a carbon nanotube. The central dataset is for water between two flat surfaces. The lines of fit shown are used to find the constants in Equation (11). 
However, for the cylindrical cases the radius has a significant impact on the amount of friction. This can be attributed to the change in curvature with radius, R. A small radius amounts to a large curvature:(10)
Not only this, but the type of curvature also impacts the amount of friction, whether concave or convex. Fluid on the outside of a CNT experiences a larger amount of friction, reducing the amount of slip. Contrary to this, the fluid on the inside of a CNT experiences a smaller amount of friction.
Another key result here is that there is a clear logarithmic approach to the base case for both fluid inside and outside of a CNT. The relationship between friction coefficient and radius takes the following form:(11)
The constants A and B can be found from Figure 3 and have the following values:(12)
The friction coefficient in all cases becomes effectively constant, at a value of 10,000 ± 600 kg m^{−2} s^{−1}, above some critical radius of curvature, approximately 7.4 ± 2.9 nm from these results. This is equivalent to a curvature of 0.14 nm^{−1}. These results are consistent with the work of Falk et al. [6].
In Figure 4, friction coefficients have been converted to slip lengths, using Equation (5), for comparison with slip lengths extracted from NEMD simulations using Equations (2) and (3). The viscosity used here was found from a 1 ns equilibrium simulation by solving Equation (6) and was found to be 0.69 ± 0.10 cp. A density of 999.12 ± 0.23 kg m^{−3} was found from the same simulation. As can be seen, there is good agreement between the NEMD results and the EMD results, confirming that this EMD method is an appropriate method for atomistically smooth and symmetrical surfaces.
Figure 4. Comparison between EMD and NEMD slip lengths calculated using Equation (5) and Equations (2) and (3) respectively. The lines are the fits from the EMD method and are equivalent to those in Figure 3. NEMD results only exist for the inside and flat cases as there is no simple analytical solution to the NavierStokes equation for the outside case. 
The same curvature dependence of the friction coefficient has been found for systems involving decane, ethanol and octamethylcyclotetrasiloxane [26]. One exception is in the case of decane, which shows a significant reduction in the friction coefficient for small slab separations. This is due to the strong layering of decane at the solid surface. These layers influence each other and form stable domains which impact the decane structure at the surface and thus the decane/graphene commensurability. In the cylindrical cases, the stiff chain structure of the decane molecules results in an orientation that is preferentially parallel to the cylindrical walls which becomes less preferential as the radius of curvature increases. This is also in agreement with works from Supple and Quirke [27]. The works cited above suggest that the simple functional relationship between friction coefficient and radius of curvature shown here is not applicable for more complicated molecules. Although, it may be good enough for simulating systems where most surfaces have some curvature, or for relatively large surface separations (>6 nm). When scaling up simulations from MD to LB for long chain molecules, such as decane, alterations to the functional fit may be required.
Chen et al. [28] calculated the slip dependence on curvature directly by measuring the velocity profile of fluid surrounding rotating cylinders, fitting it to the NavierStokes equation, and extrapolating to get a slip length. They use a simple, single atom, LennardJones fluid. They found that the slip length at a flat wall, given the same atomic interactions, is the same as for all other cases, including curved and rotating boundaries. This is not necessarily contrary to this work as the systems studied are very different. In their work, they used a solid surface with an atomic number density equal to that of the fluid, this leads to a surface which is approximately 30 times less dense than in the present work. This difference in surface could have a significant impact on the results, as the likely reason for the deviation in the amount of slip seen between flat and curved surfaces is the smoothing of the surface energy landscape experienced by the fluid in proximity to the surface with increasing concave curvature and the roughening exhibited with increasing convex curvature. This is discussed in detail by Falk et al. [6]. In addition to this, the smallest radius of curvature tested is approximately 3.5 nm for flow external to a cylinder and 16 nm for internal flow. For internal flow, concave curvature, the radius is already much larger than the critical curvature for TIP3P water in the graphitic systems studied in this work. However, for flow externally, 3.5 nm is within the range at which curvature should have an impact on the slip length. The lack of a change seen in their work is likely due to the differences discussed above between the surface density and the fluid model. It is unknown for sure where the critical curvature would be for less dense surfaces, however it is our suspicion that it would be smaller or nonexistent due to the inability to reach such small cylinders with such a low solid density. The critical radius of curvature would also shift depending on the wettability of the surface and structure of the fluid, as seen in the work of Falk et al. [26]. Another difference between this study and theirs is that they use a solid surface which can hold a temperature, unlike the unmoving surface in this work, which may result in a change in the mechanism of slip.
Petravic and Harrowell [29] pointed out issues with the Bocquet and Barrat method of calculating the friction coefficient and presented an alternative equilibrium theory of the friction associated with the confined fluid and showed how this is related to the intrinsic interfacial friction. They showed that for an amorphous liquidlike wall the friction coefficient, calculated using Bocquet and Barrat method, was no longer constant with separation in the slab case even though NEMD calculations show otherwise. Amorphous surfaces have not been testing in this current work and as such no comment can be made on the applicability of this method in those cases. However, here it has been shown that for graphitic surfaces and CNT the resulting friction coefficient does correspond with that of NEMD simulations, see Figure 4.
The evolution of fluid density within each pore with time can be seen in Figure 5. As expected, the larger systems take longer to fill, however all systems shown here were at their mean fluid densities before the results run. The two data sets that exhibit a lower mean equilibrium density belong to the systems of fluid inside the 1 and 1.5 nm CNT. The reason for this is explained later.
Figure 5. Evolution of the density within each pore during the equilibration process. Only data after 500 ps is used for calculation of the friction coefficient. 
It should be noted that the equilibration process did not involve any barostatting other than the use of pistons. These final mean densities in the pores were compared to simulations that used a Berendsen et al. [30] barostat. The resulting densities were close to those found here, with no bias in either direction. However, the results presented here are for liquid water only. The choice of method may have some impact on gas phase systems.
Figure 6 shows the water density profile for all simulations. For fluid on the outside of a CNT the density profiles all align perfectly and so only one line is shown. For both the flat surfaces and the inside simulations the larger cylinders resulted in overlapping density profiles, however deviation is seen for lower radii CNT, starting at 1.5 nm. This is due to the disappearance of a bulk region of fluid in the centre of the channel. Below this radius, the effects of the wall is seen by fluid at the other side of the pore. This causes a change in the density profile and may have an effect on the calculation of the friction coefficient.
Figure 6. Density profiles of water within each system. 
The statistical accuracy is also questionable down at this radius, as the number of molecules is significantly reduced. There are approximately 330 molecules in the 1.5 nm CNT. In future calculations, the CNT will be made longer to accommodate more molecules at these small radii.
Conclusion
In this study, we have performed MD simulations to find the friction coefficient as a function of the concave and convex surface curvatures.
Given the clear functional form of this friction coefficient with curvature it should be fairly straightforward to use these results in upscaled simulations. For example, this could include Lattice Boltzmann simulations of more complex geometries consisting of multiple pores of varying shapes and sizes. This requires the implementation of a nonzero velocity boundary condition at solidfluid interfaces. The local surface curvature also needs to be calculated, however this can easily be done during preprocessing. The local shear rate can be extracted trivially from the local distribution functions. Combined with the slip lengths found in this study, the wall velocity can be found explicitly from Equation (1) and applied in the simulation. These larger scale simulations can then further our understanding of nanofluidics with respect to complex geometries and can be used to calculate macroscale properties such as permeability.
Acknowledgments
This research is funded by EPSRC and Shell Research Limited, UK.
References
 Yang J., Crawshaw J., Boek E.S. (2013) Quantitative determination of molecular propagator distributions for solute transport in homogeneous and heterogeneous porous media using lattice Boltzmann simulations, Water Resources Research 49, 12, 8531–8538. [CrossRef] [Google Scholar]
 Yang J., Boek E.S. (2013) A comparison study of multicomponent Lattice Boltzmann models for flow in porous media applications, Computers & Mathematics with Applications 65, 6, 882–890. [CrossRef] [MathSciNet] [Google Scholar]
 Boek E.S., Venturoli M. (2010) LatticeBoltzmann studies of fluid flow in porous media with realistic rock geometries, Computers and Mathematics with Applications 59, 7, 2305–2314. [CrossRef] [MathSciNet] [Google Scholar]
 Stukan M.R., Ligneul P., Crawshaw J.P., Boek E.S. (2010) Spontaneous imbibition in nanopores of different roughness and wettability, Langmuir 26, 16, 13342–13352. [CrossRef] [PubMed] [Google Scholar]
 Stukan M.R., Ligneul P., Boek E.S. (2012) Molecular Dynamics Simulation of Spontaneous Imbibition in Nanopores and Recovery of Asphaltenic Crude Oils Using Surfactants for EOR Applications, Oil & Gas Science and Technology – Revue d’IFP Energies nouvelles 67, 5, 737–742. [Google Scholar]
 Falk K., Sedlmeier F., Joly L., Netz R.R., Bocquet L. (2010) Molecular Origin of Fast Water Transport in Carbon Nanotube Membranes: Superlubricity versus Curvature Dependent Friction, Nano Letters 10, 10, 4067–4073. [CrossRef] [PubMed] [Google Scholar]
 Koplik J., Banavar J.R., Willemsen J.F. (1989) Molecular dynamics of fluid flow at solid surfaces, Physics of Fluids A: Fluid Dynamics 1, 781. [CrossRef] [Google Scholar]
 Barrat J.L., Bocquet L. (1999) Large Slip Effect at a Nonwetting FluidSolid Interface, Physical Review Letters 82, 23, 4671–4674. [Google Scholar]
 Cieplak M., Koplik J., Banavar J. (2001) Boundary Conditions at a FluidSolid Interface, Physical Review Letters 86, 5, 803–806. [CrossRef] [PubMed] [Google Scholar]
 Galea T.M., Attard P. (2004) Molecular dynamics study of the effect of atomic roughness on the slip length at the fluidsolid boundary during shear flow, Langmuir 20, 8, 3477–3482. [CrossRef] [PubMed] [Google Scholar]
 Sokhan V.P., Nicholson D., Quirke N. (2001) Fluid flow in nanopores: An examination of hydrodynamic boundary conditions, The Journal of Chemical Physics 115, 8, 3878. [CrossRef] [Google Scholar]
 Sokhan V.P., Nicholson D., Quirke N. (2002) Fluid flow in nanopores: Accurate boundary conditions for carbon nanotubes, The Journal of Chemical Physics 117, 18, 8531. [CrossRef] [Google Scholar]
 Kannam S.K., Todd B.D., Hansen J.S., Daivis Peter J. (2013) How fast does water flow in carbon nanotubes? The Journal of Chemical Physics 138, 9, 094701. [CrossRef] [PubMed] [Google Scholar]
 Hansen J.S., Todd B.D., Daivis P.J. (2011) Prediction of fluid velocity slip at solid surfaces, Physical Review E 84, 1, 016313. [CrossRef] [Google Scholar]
 Pronk S., Páll S., Schulz R., Larsson P., Bjelkmar P., Apostolov R., Shirts M.R., Smith J.C., Kasson P.M., van der Spoel D., Hess B., Lindahl E. (2013) GROMACS 4.5: a highthroughput and highly parallel open source molecular simulation toolkit, Bioinformatics (Oxford, England) 29, 7, 845–854. [CrossRef] [PubMed] [Google Scholar]
 Thompson P.A., Troian S.M. (1997) A general boundary condition for liquid flow at solid surfaces, Nature 389, 6649, 360–362. [CrossRef] [Google Scholar]
 Bocquet L., Barrat J.L. (2007) Flow boundary conditions from nano to microscales, Soft Matter 3, 6, 685–693. [CrossRef] [Google Scholar]
 Zwanzig R. (1965) TimeCorrelation Functions and Transport Coefficients in Statistical Mechanics, Annu. Rev. Phys. Chem. 16, 67–102. [CrossRef] [Google Scholar]
 Jorgensen W.L., TiradoRives J. (1988) The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin, Journal of the American Chemical Society 110, 6, 1657–1666. [CrossRef] [PubMed] [Google Scholar]
 Berendsen H.J.C., Grigera J.R., Straatsma T.P. (1987) The missing term in effective pair potentials, The Journal of Physical Chemistry 91, 24, 6269–6271. [Google Scholar]
 Jorgensen W.L., Chandrasekhar J., Madura J.D., Impey R.W., Klein M.L. (1983) Comparison of simple potential functions for simulating liquid water, The Journal of Chemical Physics 79, 2, 926. [Google Scholar]
 MATLAB (2014) version 8.3 (R2014a). The MathWorks Inc., Natick, Massachusetts. [Google Scholar]
 Bussi G., Donadio D., Parrinello M. (2007) Canonical sampling through velocity rescaling, The Journal of Chemical Physics 126, 1, 014101. [Google Scholar]
 Nosé S. (1984) A molecular dynamics method for simulations in the canonical ensemble, Molecular Physics 52, 2, 255–268. [Google Scholar]
 Hoover W. (1985) Canonical dynamics: Equilibrium phasespace distributions, Physical Review A 31, 3, 1695–1697. [Google Scholar]
 Falk K., Sedlmeier F., Joly L., Netz R.R., Bocquet L. (2012) Ultralow Liquid/Solid Friction in Carbon Nanotubes: Comprehensive Theory for Alcohols, Alkanes, OMCTS, and Water, Langmuir 28, 40, 14261–14272. [CrossRef] [PubMed] [Google Scholar]
 Supple S., Quirke N. (2005) Molecular dynamics of transient oil flows in nanopores. II. Density profiles and molecular structure for decane in carbon nanotubes, The Journal of Chemical Physics 122, 10, 104706. [CrossRef] [PubMed] [Google Scholar]
 Chen W., Zhang R., Koplik J. (2014) Velocity slip on curved surfaces, Physical Review E 89, 2, 023005. [CrossRef] [Google Scholar]
 Petravic J., Harrowell P. (2007) On the equilibrium calculation of the friction coefficient for liquid slip against a wall, The Journal of Chemical Physics 127, 17, 174706. [CrossRef] [PubMed] [Google Scholar]
 Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., DiNola A., Haak J.R. (1984) Molecular dynamics with coupling to an external bath, The Journal of Chemical Physics 81, 8, 3684. [Google Scholar]
Cite this article as: D.A. Ross and E.S. Boek (2016). Molecular Dynamics Simulations of Slip on Curved Surfaces, Oil Gas Sci. Technol 71, 46.
All Tables
All Figures
Figure 1. Top left: Water filled, single layered, armchair carbon nanotube. Top right: Single layered armchair carbon nanotube with water on the outside. Bottom left: Water confined between two symmetrical, single layered, graphitic plates. Bottom right: Graphitic structure. 

In the text 
Figure 2. Cross section of equilibrated system of water inside a carbon nanotube, before removal of the reservoirs. Barriers are included during equilibration only. 

In the text 
Figure 3. Friction coefficient calculated using Equation (4) for water inside and outside of a carbon nanotube. The central dataset is for water between two flat surfaces. The lines of fit shown are used to find the constants in Equation (11). 

In the text 
Figure 4. Comparison between EMD and NEMD slip lengths calculated using Equation (5) and Equations (2) and (3) respectively. The lines are the fits from the EMD method and are equivalent to those in Figure 3. NEMD results only exist for the inside and flat cases as there is no simple analytical solution to the NavierStokes equation for the outside case. 

In the text 
Figure 5. Evolution of the density within each pore during the equilibration process. Only data after 500 ps is used for calculation of the friction coefficient. 

In the text 
Figure 6. Density profiles of water within each system. 

In the text 