Molecular dynamics pre-simulations for nanoscale computational fluid dynamics

We present a procedure for using molecular dynamics (MD) simulations to provide essential fluid and interface properties for subsequent use in computational fluid dynamics (CFD) calculations of nanoscale fluid flows. The MD pre-simulations enable us to obtain an equation of state, constitutive relations, and boundary conditions for any given fluid/solid combination, in a form that can be conveniently implemented within an otherwise conventional Navier–Stokes solver. Our results demonstrate that these enhanced CFD simulations are then capable of providing good flow field results in a range of complex geometries at the nanoscale. Comparison for validation is with full-scale MD simulations here, but the computational cost of the enhanced CFD is negligible in comparison with the MD. Importantly, accurate predictions can be obtained in geometries that are more complex than the planar MD pre-simulation geometry that provides the nanoscale fluid properties. The robustness of the enhanced CFD is tested by application to water flow along a (15,15) carbon nanotube, and it is found that useful flow information can be obtained.


Introduction
Nanofluidic technologies are advancing rapidly, and a range of new technical opportunities are emerging, for example, efficient filtration of water using carbon nanotubes (CNT) (Alexiadis and Kassinos 2008;Mattia and Gogotsi 2008), heat removal and control in high heat-flux systems such as nuclear reactors, micro/nano-electro-mechanical systems (MEMS/NEMS), and micro-chemical reactors (Saidur et al. 2011;Wen et al. 2009). However, the prediction of the fluid mass flow rate and heat transfer in nanoscale systems presents a major barrier to their design. The existence of non-continuum effects, such as molecular layering and velocity slip near to liquid-solid interfaces, seemingly precludes efficient continuum computational fluid dynamics (CFD). On the other hand, more accurate molecular dynamics (MD) simulations can be extremely costly in terms of the computational resources they require. For example, to simulate 200 nm 3 of water for 1 ns, with a modern MD code running on 8 CPUs in parallel, can require approximately 2 days of computational time. To simulate the liquid over much larger time and space scales is beyond the reach of current computational capabilities (even with the improved computational efficiency offered by graphical processing units, GPUs). This certainly prevents using such simulations within a practical iterative design process.
In this paper, we demonstrate that CFD simulations can still play a major role in nano-flow prediction for engineering design or, alternatively, for the initiation of steadystate MD simulations. We show that useful predictions can be readily and reliably obtained in complex nanoscale geometrical domains, despite the limits of the continuumfluid assumption, if appropriate, fluid state models, viscosity relationships, and slip models are provided.
Both experiments and molecular simulations show that the continuum-fluid assumptions (e.g. local linear constitutive relationships) are still appropriate for water confined to channels of width *1-2 nm (see Bocquet and Charlaix (2010) and references therein), and MD simulations have been used to show that Lennard-Jones (LJ) fluids confined to geometries of *2-3 nm still show continuum behaviour (Huang et al. 2007;Sofos et al. 2009;Travis et al. 1997). For nano-confined fluids below the continuum limit, however, the strain rate can vary rapidly within several molecular diameters  due to oscillations in the density; therefore, the stress becomes non-local. For quickly varying strain rates in homogeneous fluids, the non-local stress is calculable , but for confined fluids, this is an unsolved problem [see for example Cadusch et al. (2008); Todd (2005)]. In addition, tests of continuum-fluid equation performance at the nanoscale have been restricted to extremely simple geometries, typically Poiseuille flow [e.g. Zhang et al. (2012)] and other canonical cases for which Hagen-Poiseuille equations are solved with slip boundary conditions. There is currently little substantive evidence to suggest that CFD can generally be applied at the nanoscale in arbitrary flow geometries.
In this paper, we make three contributions to nanofluidic modelling: 1. we propose a convenient pre-simulation MD framework that enables us to measure CFD-type properties (e.g. boundary conditions and constitutive relations) for a given solid/liquid combination; 2. we demonstrate that these fluid properties can be used to obtain highly accurate CFD predictions in complex geometries at the nanoscale, provided that a significant portion of the flow exhibits continuum bulk-like behaviour. This is distinct from previous studies [e.g. Zhang et al. (2012)] that focused mainly on 1D flow configurations. Importantly, we demonstrate that accurate CFD predictions can be obtained in cases that are more complex than the MD flow configuration from which the CFD fluid parameters were obtained; 3. we demonstrate that even in highly non-continuum 3D flows, where no significant bulk flow exists (such as flow through some small diameter CNT), qualitatively accurate, and so useful, CFD predictions can be obtained.
The paper is structured as follows: in Sect. 2, we define the MD procedure from which the fluid properties are extracted for a particular liquid/solid combination. In Sect. 3, this is applied to a LJ fluid interacting with solid bounding surface atoms that are effectively frozen. This is followed by nanoscale CFD predictions for three different geometries: a short channel connecting two reservoirs; a long channel connecting two reservoirs; and a channel with a geometrical irregularity connecting two reservoirs. The results in each case are compared with full MD data, which are vastly more computationally expensive and time-consuming to obtain. In Sect. 4, to demonstrate its applicability to highly non-continuum flows, we use our enhanced CFD to model water flow along a CNT. This example is chosen in order to explore the robustness of our approach for cases where the continuum assumption is known to be invalid. Finally, in Sect. 5, we discuss the applications, limitations, and potential developments of the method.
2 Molecular dynamics procedure for generating fluid properties We employ a preliminary MD simulation to obtain fluid properties and boundary conditions that enable the effective use of a Navier-Stokes solver in nanoscale applications it is typically not suitable for. This approach could be classed as a 'sequential molecular-continuum hybrid method' [see Mohamed and Mohamad (2010) for a review of hybrid methods], where 'sequential' refers to the fact that the MD is performed in advance of, and so independent of, the continuum model. In contrast, in 'concurrent' schemes, e.g. HMM (E et al. 2009;Asproulis et al. 2012) or Domain Decomposition (O'Connell and Thompson 1995), the MD simulations are fully coupled to a continuum model. While this is ultimately more likely to produce accurate solutions, it is certainly more computationally expensive. The sequential methodology we adopt instead is far more practical for preliminary and iterative engineering design simulations (at least, based on today's computing capabilities). An example of a similar methodology can be found in Dongari et al. (2011a, b) in which molecular free path distributions within rarefied gases are measured using MD simulations and then used within the Navier-Stokes-Fourier equations. For the isothermal CFD simulations of nano-flows we consider in this paper, we require the following fluid properties and boundary conditions: the viscosity coefficient as a function of density, the pressure as a function of density, the slip length as a function of density and shear rate, and what we here define as the 'surface offset' (d) which defines the position of the surface, as modelled by the CFD, relative to some atomistic reference point (in this paper, the atomic centres). The implications of the choice of the dependencies the fluid properties and boundary conditions on particular variables (here, density) are discussed in Sect. 5.
For efficiency and convenience, we propose a single MD configuration from which all of these fluid properties can be measured and/or controlled. This enables efficient concurrent MD simulations over any range of variables considered by multiple realisations of the same MD geometry/setup. We have chosen to use non-equilibrium MD (NEMD) simulations to gather the data rather than equilibrium MD despite potential, accuracy issues (Kannam et al. 2012), because this approach would be more suitable for simulating a complex fluid-solid interaction, where slip could be strain dependent. A possible more refined approach could be to use a combination of equilibrium MD simulations and NEMD simulations where appropriate. Figure 1 (far left) shows the MD pre-simulation domain; it is symmetrical about its centrelines and uses periodic boundary conditions in the streamwise direction (i.e. in the x-direction) and into the page (i.e. in the z-direction). The domain has bulk, shear and interface zones (as labelled) for measuring state, constitutive and boundary properties, respectively. Pressure and density (and also temperature, if simulating non-isothermal cases) are measured in the bulk zone. In addition to this, in the bulk zone, an artificial streamwise body force (F x ) is applied (see Fig. 1, centre left), which creates a velocity profile in the domain similar to that illustrated (centre right). We assume that the equation of state in the bulk zone is unaffected by the magnitude of strain rate generated. In the shear zone, the fluid is subject to a constant shear stress, s xy , directly resulting from the bulk-zone forcing. A linear flow velocity profile is developed in the shear zone, and this is leastsquares fitted. The measured strain rate and shear stress are then used to obtain a viscosity coefficient, l, through In this paper, we assume that shear viscosity is sufficient to describe the fluid constitutive behaviour, while accepting that the pre-simulation configuration would need to be modified to deal with extensional viscosity. Any significant density oscillations associated with molecular layering are confined to the interface zone. In this zone, we calculate what we term the 'CFD surface displacement' (d), which is the distance that a CFD wall/ surface needs to be displaced from the centres of surface atoms in order to accurately represent the boundary of the fluid (as opposed to the boundary of the solid); see d in Fig.  1. We take this displacement to be the distance from the centre of the surface wall atom to where the fluid density becomes at least 10 % of the bulk, i.e. q ! aq bulk , where a ¼ 0:1. The linear velocity profile obtained in the shear zone is extrapolated into the interface zone to find the apparent slip length, n, as defined from the CFD surface (see Fig. 1, centre right). Across multiple simulations, we obtain the bulk pressure, a viscosity coefficient, the slip length, and the surface displacement, for a range of combinations of bulk density and applied shear stress. In the MD simulations, the applied shear stress and bulk density are varied by modifying the body force and by adding/removing molecules, respectively, using the FADE algorithm . Finally, once all data are collected over the expected range of density and shear stress, 1 functional relationships are constructed for the desired fluid properties (using, for example, fitted polynomials), which are then used in the CFD model. The behaviour of this CFD model ultimately depends on these functional relationships, and this choice requires some experience or needs to form part of an iterative approach (this is discussed in Sect. 5). For the cases considered in this paper, we adopt the following: for pressure, p ¼ pðqÞ; for dynamic viscosity, l ¼ lðqÞ; surface displacement, d ¼ dðqÞ; and slip length n ¼ nðq; _ cÞ, where q is the bulk fluid density and _ c is the strain rate in the shear zone. This dependence on density would normally imply a high-speed high-Mach number flow, but in nanoscale internal flows, it is possible to have substantial fluid compressibility at extremely low Mach numbers due to viscous-related pressure losses [see Gad-el Hak (2010) for a discussion of this]. For this reason, capturing the influence of density on fluid properties is critical to the accurate prediction of nanoscale flows. For all of the examples considered in this paper, the influence of strain rate can be safely ignored, but we consider its effect on slip length for demonstration purposes. The fluids we consider are therefore Newtonian in the bulk; a non-Newtonian fluid, for example, would at least require l ¼ lðq; _ cÞ. Note that for the simulation of well-understood fluids, it would not be necessary to extract all of these properties from MD pre-simulations.

CFD for nanoscale flow of a Lennard-Jones fluid
Owing to the lack of detailed and reliable experimental flow measurements at the nanoscale, in this section, we compare our enhanced CFD predictions with full-scale MD simulation results. This comparison is intended to test whether flow field solutions of comparable accuracy to full MD can be obtained from our enhanced CFD in complex nanoscale geometries, without the need for ad hoc corrections, and at only a fraction of the cost of full MD. A Lennard-Jones (LJ) model of liquid argon is chosen, where the solid wall atoms are fixed/frozen (Thompson and Troian 1997); the exact interatomic potentials used are given in '' Appendix''.

MD pre-simulation results
The MD pre-simulations are performed as described in Sect. 2, using the mdFoam solver (Borg et al. 2010;Macpherson et al. 2007;Macpherson and Reese 2008) that is implemented within the OpenFOAM libraries (Open-FOAM 2013). The MD algorithm has molecules evolving using Newton's equations of motion, and m i are the velocity, total force, position and mass, respectively, of an arbitrary molecule i in a system of N molecules at a time t. The total force per molecule f i is calculated at every time-step from the sum of pair-wise intermolecular forces between molecules, Uðr ij Þ is the potential energy when molecules i and j are separated by The MD pre-simulation domain is constructed as illustrated in Fig. 1 (far left) extending 4.08 nm and 5.44 nm in the x-and z-directions, respectively. These dimensions are principally chosen to be large enough to avoid unwanted 'wrap-around' effects due to the periodic boundaries. In the y-direction, each interface region extends by 0.68 nm, each shear zone extends by 3.4 nm, and the bulk zone extends by 4.08 nm, giving a total height of 12.24 nm. Owing to the symmetry of the problem, properties extracted from the shear and interface zones are mirrored and averaged (thus reducing the overall sampling time). The temperature is maintained at T = 292.8 K by coupling molecules to a velocity-unbiased Berendsen thermostat (Berendsen et al. 1984) with a time constant of s T ¼ 21:6 fs, applied within 36 independent bins placed in the y-direction, with each bin being 0.34 nm thick. It has been shown that using a thermostat on a confined fluid may affect the flow properties (Bernardi et al. 2010); however, a thermostat is used in all the MD simulations in this work so that we can compare with isothermal CFD simulations. The same thermostat is used in both the pre-simulations and the full MD simulations; therefore, the same error is present in both-the verification of the simulation approach is thus not undermined by any physical uncertainty introduced by the thermostat.
The pre-simulation is divided into two steps. In the first part of the simulation, the MD ensemble is set to the target density and allowed to run to a steady state (which takes *1.5 million MD time-steps) with the external artificial force applied for the target strain rate. For any molecule in the bulk region, the external forcing is given by a Gaussian distribution centred around the centreline of the simulation box: where " F is the magnitude of the Gaussian, r s is an estimate of the required width of the curve, andn x is the unit vector in the x-direction. The relationship between the forcing magnitude and the shear stress can be obtained by substituting Eq. (2) into the conservation of momentum and integrating giving: The second part of the pre-simulation is then used to measure the required fluid properties over an interval of around 2 million MD time-steps of 5.4 fs each.
3.1.1 Bulk pressure as a function of bulk density Figure 2a shows MD pre-simulation measurements of pressure, obtained from the standard Irving and Kirkwood expression (1950), varying with the mass density.
The MD pre-simulation results are least-squares-fitted to a second-order polynomial. This then serves as an equation of state within the enhanced CFD solver to connect the mass continuity equation to the momentum equation. In this case, the polynomial is p ¼ 0:001559q 2 À 3:387q þ 2020:6. For reference, data from the NIST database for argon (Linstrom and Mallard 2001) are also plotted in Fig. 2a and is in close agreement with our MD pre-simulation data. Clearly, in this particular case, properties for argon are well known, but we extract the equation of state from our MD presimulation for the purposes of demonstration. The equation of state (and the viscosity equation in the following subsection) could be obtained by performing equilibrium MD simulations of a bulk fluid; however, the data are extracted from the pre-simulations here both for convenience and computational efficiency.

Dynamic viscosity as a function of density
The strain rate is extracted from the MD shear zone by a least-squares linear fit to the relaxed and time-averaged velocity profile. The applied shear stress is measured using the Irving-Kirkwood equation and then compared with the strain rate using Eq. (1) to give a dynamic shear viscosity coefficient for L-J argon at a given bulk density. The viscosity coefficients measured from our MD pre-simulations of LJ argon are shown in Fig. 2b. A least-squares polynomial fit of 2nd order in density is also plotted: l ¼ 7:96 Â 10 À10 q 2 À 1:774 Â 10 À6 q þ 0:001106. This is then used in our enhanced CFD simulations to close the momentum equation. Again, for reference, data from the NIST database for liquid argon are also plotted in Fig. 2b. Note, due to the breakdown of the continuum assumption and the existence of non-local stress, this state-dependent viscosity becomes only approximate when applied to a nano-confined fluid.

CFD surface displacement as a function of density
The surface displacement d defines the location of the CFD boundaries relative to the atomic (actual) walls. If d varies substantially with density (or any other fluid property), the geometry of the enhanced CFD domain becomes dependent on the CFD solution itself. However, for the fluid/ solid combinations considered in this paper, over the density ranges considered, d is effectively constant, see Fig. 3. This allows us to assume that the surface displacement is fixed, which avoids the need for a complicated recursive solution and re-meshing procedure. In previous work, the surface displacement has been set to the liquid-solid interaction length (Joseph and Aluru 2008)  chosen in another arbitrary way, and in some cases neglected altogether. The liquid-solid interaction length is significantly larger than the surface displacement we propose, which may be partly responsible for the large discrepancy between continuum-fluid predictions and MD results found in these previous studies, e.g. Joseph and Aluru (2008).

Slip length as a function of density and strain rate
Liquid slip velocity at surfaces is calculated using the Navier slip condition: where n is the slip length and _ c is the shear rate at the bounding surface. The same least-squares-fitted linear velocity profile from Sect. 3.1.2 is used to calculate the slip length (as defined from the CFD surface). In this work, we only investigate steady isothermal liquid flows; therefore, our slip boundary condition only needs to depend on the strain rate [as in Thompson and Troian (1997)] and the density (Bocquet and Charlaix 2010). In fact, the strain rate dependence is not necessary for the example cases we consider due to the relatively low shear rates, but here it is included for the purposes of illustration. Based on the strain rate/slip length relationship proposed in Thompson and Troian (1997), and assuming a linear dependence on density, a least-squares fit is performed to the following equation: where q is the density, _ c c is the critical shear rate [see Thompson and Troian (1997)], and c 1 ; c 2 and _ c are parameters of the fit to our MD pre-simulations, which are À1:2052 Â 10 À12 kg À1 m 4 ; 3:7468 Â 10 À9 m and 1:5431Â 10 11 s À1 , respectively. We leave the slip length dependence on curvature [as discussed in Einzel et al. (1990)] for future refinement of this model. Figure 4 shows our MD pre-simulation data and the least-squares fit of Eq. (5); results are shown for three different values of density. The slip model approximated by Eqs. (4) and (5) is directly introduced as a Robin boundary condition in the enhanced CFD solver. The slip length for the simple fluids modelled in this paper could be measured using equilibrium MD using similar methods to Bocquet and Barrat (1994);Hansen et al. (2011). However, for ease of application and generality (e.g. for capturing the interaction of solid surfaces with non-Newtonian fluids) we use this non-equilibrium methodology.

The enhanced CFD model
We perform finite-volume CFD simulations using Open-FOAM (2013), an open source set of C?? libraries for solving partial differential equations on unstructured meshes and in parallel. Specifically, we use the laminar, compressible solver sonicLiquidFoam, which we have modified to (a) accommodate a nonlinear equation of state, (b) allow a density-dependent viscosity, and (c) incorporate slip boundary conditions of the form given in Eq. (5). A compressible solver is used despite the very low Mach numbers, since, as discussed above, significant compressibility can occur in microand nano-geometries (Gad-el Hak 2010; Patronis et al. 2013).

Simulation results
To test the reliability of our predictions using CFD enhanced with MD pre-simulation input, we compare them to results from full-domain MD calculations. We also compare results with predictions from compressible CFD with no-slip at the wall, and without modelling the CFD surface displacement  Fig. 4 Slip length n varying with strain rate _ c for three density values q 1 ¼ 1;276 kg/m 3 ; q 2 ¼ 1;447 kg/m 3 and q 3 ¼ 1;668 kg/m 3 . MD presimulation data points (symbols) and fit to Eq. (5) (dashed lines) (referred to as 'no-slip CFD'). We also compare with incompressible CFD with the same slip model but no surface displacement (referred to below as 'incomp. slip CFD'). As test cases, we choose flows that all exhibit non-continuum behaviour (e.g. slip at surfaces), but also contain a significant bulk-flow region, even within the smallest features of the geometry. In Sect. 4 we consider the quality of CFD predictions in cases where such a bulk region does not exist.
The two-dimensional cases we consider in this section involve connected reservoirs that are held at different pressures, an example of a filtration configuration, say. The first case has the reservoirs connected by a straight channel 108.8 nm long (Case 1), the second by a straight channel 231.2 nm long (Case 2), and the third by a 231.2 nm long straight channel with a cylindrical geometrical irregularity/ defect with radius 0.68 nm (Case 3). The channel width (measured as the distance between the centre of opposite solid surface atoms) is 4.08 nm for all three cases. Figure 5 shows the CFD meshes for Case 1 alongside the corresponding full MD domain. Figure 6 shows the CFD and MD domain for Case 3 (L ¼ 231:2 nm, with channel width of 4.08 nm; the width at the defect is 1.7 nm). For all three cases, the pressure at the inlet and outlet reservoir is 650 MPa and 300 MPa, respectively. The full MD simulations are used to evaluate the accuracy 2 of the enhanced CFD predictions by comparison. As is standard CFD practice, the mesh resolution has been tested in each case to ensure mesh independence of the results.
Figures 7a-c and 8a-c show results of pressure and density, respectively, along the centreline of each domain for each case. Both the centreline and the sampling region (where averaging is performed) are indicated on Fig. 5. In Figs. 7 and 8, differences between the CFD and MD results can be seen near the outermost boundaries of the reservoirs. This is because in the full-domain MD, for convenience, the reservoirs are connected by periodic boundary conditions, with a local body force imposing the pressure drop [see Docherty et al. (2014) for details of this approach]; in the CFD, however, boundary pressures can be specified directly, and so periodicity need not be enforced.
Velocity profiles cross-channel are presented in Fig. 9a-d at cross sections A and B (as indicated on Fig. 5) for Cases 1 and 2. The streamwise velocity along the centreline of the channel for Case 3 is presented in Fig. 9e. Finally, in Table 1, predictions for the mass flow rate through each channel are given.
In all three cases, the agreement between our enhanced CFD model (the dashed lines in Figs. 7, 8 and 9) and the MD results (solid lines) is extremely good for all of the flow variables considered. Also, CFD predictions of mass flow rate (arguably the most important bulk property in nano-channel flow cases, and one that no-slip CFD underpredicts very substantially) are all within 4 % of the   values obtained from full MD simulations. This very positive result is reassuring given the non-trivial nature of the geometry considered in Case 3, with the small nonplanar irregularity in the channel. The MD simulation results show this small defect reduces the mass flow rate by more than 10 % (compared with the otherwise identical Case 2). Again, our enhanced CFD technique captures this effect accurately: the flow rate is reduced by 12 %. Table 2 provides an indication of the computational cost for the three full-domain MD simulations. The longest simulations presented in this paper ran in parallel (on 24 CPUs) for 18 days. The laminar-flow CFD itself has a negligible cost by comparison, although the MD pre-simulations also require the computational resources indicated in the last row of Table 2. However, these pre-simulations need only to be performed once for a particular fluid/solid  The percentage difference (error) between the mass flow rates predicted by the CFD models and the full MD results are presented in parentheses Microfluid Nanofluid (2015) 18:461-474 469 combination and then used for any number of flow geometries thereafter.

Water flow through carbon nanotube (CNT) membranes
In this section, we test the robustness of our enhanced CFD technique for cases where, in some region of the flow field, the continuum-fluid assumption is far from being valid, and where there also exist large regions of bulk fluid for which MD is prohibitively expensive. The three-dimensional flow configuration we consider is essentially the same as depicted in Fig. 5, except that the two reservoirs of water, held at different pressures, are now separated by a (15,15) CNT of length 50 nm and diameter approximately 2 nm; since the domain is periodic in the y-and z-directions, this setup represents a regularly repeated array of CNTs. The flow of water through CNTs has recently been the focus of substantial research effort (Alexiadis and Kassinos 2008) mainly due to the extremely high flow rates that have been both predicted (Nicholls et al. 2012) and measured (Mattia and Gogotsi 2008;Whitby and Quirke 2007). These flow rates are often expressed as an enhancement factor, which is the ratio of water flow rates along the CNT to those predicted by classical fluid dynamics (i.e. the Hagen-Poiseuille equation). The low friction associated with this water transport, and the high selectivity of CNTs, makes CNTs (and other nanotubes) excellent candidates for highefficiency desalination and other filtration applications. The high flow rates, often reported as being orders of magnitude greater than classical flow theory predicts (Whitby and Quirke 2007), are typically attributed to both weak surface-fluid interactions and molecular ordering/layering that enables water molecules to pass efficiently along the CNT in a semi-ordered or structured manner. Clearly, this kind of flow is difficult to describe accurately with a continuum-fluid model, but we demonstrate below that reasonable results can still be obtained for some spatially and temporally averaged properties. The likely reason that our enhanced CFD estimates are reasonable is that the flow in the CNT is dictated by the liquid interaction with the smooth graphitic surface (which is adequately modelled), despite the non-continuum conditions within the fluid.

MD pre-simulation results
The MD pre-simulations for this case are constructed identically to these of Sect. 3, including the dimensions of the geometry. Here, though, the TIP4P/2005 molecular water model is used to describe the condensed phase of water, while the solid boundary walls consist of atom-thick graphene sheets that are modelled using 663 frozen carbon atoms. The exact interatomic potentials used are given in '' Appendix''. As water is a well-known fluid, we use data from NIST (Linstrom and Mallard 2001) for the pressuredensity and density-viscosity relationships, both of which are fitted to quadratic polynomials. For the pressure-density relationship, the equation used is p ¼ 0:00684q 2 À 11:49q þ 4655, and for the density-viscosity relationship we use: l ¼ 1:413 Â 10 À8 q À 2:879 Â 10 À5 q þ 0:01555.
MD pre-simulations provide values for the surface displacement (i.e. from the carbon atoms to the fluid) and the slip length over the range of densities within the channel. As the strain rate is low, we only use three data points to model a linear dependency of the slip length on density, i.e.  where c 1 and c 2 are parameters of the fit and are À2:8248 Â 10 À10 kg À1 m 4 and 3:3117 Â 10 À7 nm, respectively. In this case, for the enhanced CFD simulations we take d ¼ 0:266 nm and the slip length relationship is shown in Fig. 10.

Simulation results
We again compare our enhanced CFD predictions against full MD results and against the standard CFD models outlined in Sect. 3. The CFD mesh is chosen to be fine enough to safely give mesh-independent results for mass flow rate; given that the cost of the CFD simulations is extremely small, achieving this poses no particular problem. The pressure difference between the reservoirs is set to be 200 MPa because it is very challenging to obtain useful information from MD using only low pressure differences due to the extended sampling times required to filter low-velocity signals from the thermal noise (Nicholls et al. 2012). These high pressure (and consequently density) differences make the CFD predictions even more challenging. The full MD simulation was performed in parallel on 48 CPUs, with the majority of the computational effort attributable to the two reservoir regions. These have dimensions 4:4 Â 10:6 Â 10:3 nm and are chosen to be large enough to avoid any effects on the CNT flow due to reservoir boundaries. The intermolecular potentials used are the same as those in the pre-simulation, as given in '' Appendix''. Figure 11 shows pressure and density plots along the centreline of the CNT; in the MD, this is done within a cylinder of radius 0.1577 nm about the centreline. Due to the substantial density fluctuations within the MD simulations (see Fig. 11b), a bulk density effectively does not exist, and the choice of the size of this sampling region can substantially affect the bulk density measured. The no-slip CFD model does not exhibit large pressure drops at the inlet and the outlet due to the much lower velocity in the tube (and therefore, there are lower accelerations at the inlet and outlet) than in the slip cases. Cross-sectional velocity profiles in the centre of the CNT are plotted in Fig.  12a. The mass flow rate in the full MD simulation is measured to be 4:3 Â 10 À14 kg/s, which is 23 % greater than that predicted by our enhanced CFD. That this is a significant improvement on conventional CFD model predictions is indicated in Table 3.
Given the $ 2 nm diameter of the (15,15) CNT and despite the molecular layering that actually occurs within the flow field, as evidenced in Fig. 12b, our enhanced CFD approach can be considered reasonably robust in predicting important averaged fluid properties to the correct order of magnitude. These CFD results are obtained with negligible cost in comparison with full MD simulations.

Discussion and conclusions
A new procedure for solving nanoscale flows using CFD has been presented. The state, constitutive, and boundary condition information for the CFD solver is extracted from MD pre-simulations. We have demonstrated that this enhanced solver can then provide good predictions for a range of nanoscale flow geometries. A number of questions and possibilities now arise. What happens when CFD is applied far beyond the limits of its applicability? For example, how robust is the predictive performance of CFD at the nanoscale? These questions have been addressed, at least to some extent, by our results for water flow along a CNT. Our answer is CFD can be more robust than perhaps is often implied in the literature. A deeper investigation into how CFD and the continuum-fluid model perform at the limits of their applicability is needed; unfortunately, this may be restricted by the need to assess their accuracy by comparison with expensive MD simulations.
Another natural question relates to how these CFD simulations should be used given that, by comparison with MD, they are computationally cheap. It is not the case that CFD can replace MD for nanoscale simulations (in the same way that MD cannot replace experiment). However, there are a number of situations in which CFD enhanced with MD pre-simulations can be an invaluable affordable alternative or addition to MD. For example, in iterative conceptual design, where multiple simulations with slightly modified geometries are required; in initialising full-scale MD simulations, which would otherwise need to be simulated for a much longer time in order for the flow to develop from a stationary to a steady state (Kalweit 2008); and in helping to locate far-field and symmetry boundaries in full-scale MD simulations, such that their influence is not felt in the flow region of interest.
There is also the possibility that enhanced CFD could be used in some cases to produce more realistic predictions than MD. Quite often MD simulations are performed at much higher velocities (orders of magnitude higher) than would be seen in reality, solely for the purpose of increasing the signal to noise ratio (as is the case for the CNT simulations in Sect. 4). These simulations rely on the assumption that the system behaves linearly up to the extreme condition. While a full MD simulation at realistic velocities is currently extremely challenging [although not intractable (Wang et al. 2012)], a single MD pre-simulation for one fluid/solid combination is far less so. This would enable enhanced CFD to be used to investigate whether the linear response assumption is likely to be valid for a particular configuration.
A criticism of this CFD approach is that it requires an assumption beforehand about the flow and fluid behaviour. For example, the viscosity coefficient for a fluid may depend on a multitude of fluid variables with a variety of functional forms; in the examples of this paper, we have assumed these functional forms based on our experience. This is true also of hybrid particle/continuum methods in The percentage difference (error) between the mass flow rates predicted by the CFD models and the full MD results are presented in parentheses general and is not unique to the method we propose. The molecular-based simulations of HMM (Ren and E 2005), for example, have to be 'constrained' by the overall continuum model; the choice of how the constraint is performed (i.e. what variables are to be imposed on the MD subdomain) requires similar suppositions about how the fluid will behave. Also, in the present paper, we have assumed that the channels are homogeneously filled with fluid before the simulation begins. This neglects the multiphase, transient problems that occur in the fill-up process of a CNT or the flow through nano-pores, such as an aquaporin. Currently, these types of problems could not be solved by our enhanced CFD. An additional advantage of the approach in this paper is that it can be deployed recursively (and not necessarily for the same CFD simulation, or by the same researcher/designer). For example, a basic pre-simulation could be used to make a first-estimate CFD prediction; subsequent MD simulations could then be used to refine and finesse the fluid and interface models, thereby producing successively more accurate CFD predictions. Users would need to approach this refinement and finessing of fluid property models for enhanced CFD in the same way as they would a conventional mesh-dependency study.