© D.A. Ross and E.S. Boek, published by IFP Energies nouvelles, 2016

Introduction

When studying nano-porous media, conventional Computational Fluid Dynamics (CFD) and Lattice-Boltzmann (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 no-slip boundary condition at fluid-solid 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 solid-fluid combination is a challenge that requires explicit molecular scale simulation techniques such as Molecular Dynamics (MD). Previously, we have examined imbibition in nano-pores [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 non-trivial. 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 solid-fluid 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 time-step to reflect some acceleration. The result, after a long equilibration period, is some overall flow through the pore. The steady-state velocity profile can then be extracted via binning and ensemble averaging over many time-steps. 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:$$b=\frac{u_{\mathrm{wall}}}{{\mathrm{du}/\mathrm{dx}|}_{\mathrm{wall}}}$$(1)

However, this method has its weak points. For instance, the time required for equilibration is large, in that steady-state 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 5-10 repeated simulations in order to ensure sufficient statistical accuracy. Given that roughly 10-20 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 Navier-Stokes 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:$$b=\frac{2\mu u_m}{\rho \mathrm{aR}}-\frac{R}{4}$$(2) $$b=\frac{2\mu u_m}{\rho \mathrm{aH}}-\frac{H}{6}$$(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 Non-Equilibrium MD (NEMD) results for comparison with the equilibrium method described below. Here, a relatively low acceleration of 0.0001 nm ps⁻² 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 time-steps) 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 slip-length. The advantage of this is that equilibration times are significantly reduced due to no need for non-equilibrium 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 Green-Kubo type relation:$$\lambda =\frac{1}{A{k}_{B}T}{\int }_0^{\mathrm{\infty }}\mathrm{}\mathrm{dt}\langle F\left(t\right)F\left(0\right){\rangle }_{\mathrm{equ}}$$(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:$$b=\frac{\mu }{\lambda }$$(5)

In Equation (5), μ is the shear viscosity, which can be found from a single-phase fluid simulation using a similar Green-Kubo type expression as used for the friction coefficient [18]:$$\mu =\frac{V}{k_{B}T}{\int }_0^{\mathrm{\infty }}\mathrm{}\mathrm{dt}\langle P_{\alpha \beta }\left(t_0\right)P_{\alpha \beta }\left(t_0+t\right)\rangle$$(6)

where P_αβ are the off-diagonal 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 OPLS-AA (Optimized Potentials for Liquid Simulations - All Atom) forcefield [19] is used here.

Non-bonded interactions are handled by a truncated and shifted Lennard-Jones potential, with a cut-off length, r_c, of 0.9 nm. For atom-atom 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 Lennard-Jones potential, U_LJ, is calculated as follows:$$U_{\mathrm{LJ}}=4{\epsilon }_{\mathrm{ij}}\left[{\left(\frac{{\sigma }_{\mathrm{ij}}}{r}\right)}^{12}-{\left(\frac{{\sigma }_{\mathrm{ij}}}{r}\right)}^6\right]$$(7)

where σ_ij is the distance between atoms i and j at which the inter-particle potential is zero and ε_ij is the depth of the potential well which dictates the strength of the interaction. Non-bonded pair interactions are calculated as the geometric average, as is the normal case for the OPLS forcefield:$${\sigma }_{\mathrm{ij}}={\left({\sigma }_{\mathrm{ii}}{\sigma }_{\mathrm{jj}}\right)}^{1/2}$$(8) $${\epsilon }_{\mathrm{ij}}={\left({\epsilon }_{\mathrm{ii}}{\epsilon }_{\mathrm{jj}}\right)}^{1/2}$$(9)

Bonded interactions are handled with harmonic stretching and angle potentials. There are no dihedral bonds in this system. See OPLS-AA forcefield [19] for further details on force constants and parameters used in the water model.

Initially, we start with Carbon-NanoTube (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 long-chain hydrocarbon mixtures in place of water.

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.

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⁻¹ 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 pore-space. This allows for a bulk fluid region between the piston and the pore of at least 2 nm, preventing the formation of a pseudo-pore.

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 re-equilibrated 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.

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 velocity-rescaling thermostat [23] although an alternative to this is the Nosé-Hoover thermostat [24, 25]. Both of these result in the correct canonical ensemble, although velocity-rescaling was found to provide slightly improved performance.

During the results run, coordinates are output every 5 time-steps (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.

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:$$\kappa =\frac{1}{R}$$(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:$$\lambda =A{\log }_{10}R+B$$(11)

The constants A and B can be found from Figure 3 and have the following values:$$\begin{array}{c}\begin{array}{c}{A}_{\mathrm{inside}}=7,200\pm 300\\ B_{\mathrm{inside}}=3,800\pm 100\end{array}\\ A_{\mathrm{outside}}=-21,900\pm 6,700\\ B_{\mathrm{outside}}=25,300\pm 2,700\end{array}$$(12)

The friction coefficient in all cases becomes effectively constant, at a value of 10,000 ± 600 kg m⁻² s⁻¹, 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⁻¹. 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⁻³ 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.

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 Navier-Stokes equation, and extrapolating to get a slip length. They use a simple, single atom, Lennard-Jones 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 non-existent 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 liquid-like 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.

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.

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 non-zero velocity boundary condition at solid-fluid 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

All Tables

All Figures