Challenges in Nanofluidics — Beyond Navier – Stokes at the Molecular Scale

The fluid dynamics of macroscopic and microscopic systems is well developed and has been extensively validated. Its extraordinary success makes it tempting to apply Navier–Stokes fluid dynamics without modification to systems of ever decreasing dimensions as studies of nanofluidics become more prevalent. However, this can result in serious error. In this paper, we discuss several ways in which nanoconfined fluid flow differs from macroscopic flow. We give particular attention to several topics that have recently received attention in the literature: slip, spin angular momentum coupling, nonlocal stress response and density inhomogeneity. In principle, all of these effects can now be accurately modelled using validated theories. Although the basic principles are now fairly well understood, much work remains to be done in their application.


Introduction
The classical Navier-Stokes theory describing flow of Newtonian fluids has been remarkably successful, but it is inadequate under certain conditions.If the rate of deformation is high, nonlinear effects such as a shear rate dependent viscosity and normal stress differences may become apparent.If the frequency of oscillatory deformation is high, we may observe viscoelastic effects related to the elastic storage and release of energy.These effects are now quite well understood and can be described using standard treatments of non-Newtonian fluid mechanics [1].Generally speaking, the theory of fluid flow at macroscopic and microscopic scales is successful and well developed.However, when Navier-Stokes fluid dynamics and its extensions to shear rate dependent and frequency dependent constitutive relations are applied to nano-confined flows, the theory can fail.New physical effects become important and serious errors and inconsistencies can arise if the methods of macroscopic fluid mechanics are used without modification.
The reason for this is that several effects that are negligible or absent in macroscopic flows may become significant or even dominant in nanoflows.Wall slip, spin angular momentum coupling, spatially nonlocal response, and nonlinear, nonlocal coupling of the shear pressure and velocity gradient to a strongly inhomogeneous density field can all modify the flow of a fluid under strong confinement.Recent studies have shown that all of these effects can be successfully modelled.Experimental studies of fluid flow are steadily being extended to smaller and smaller system sizes, but computer simulations have a unique advantage in studies of nanoconfined fluid flow.They are ideally suited to studies of nanoflows because the system sizes and timescales that are accessible in computer simulations match the relevant size and timescales.In addition, computer simulations allow us to study velocity, density and shear pressure profiles at resolutions that would be impossible to achieve experimentally.
In this paper, we review some recent advances in the fluid dynamics of nanoconfined liquids, placing them in context and discussing the conditions under which they become important.In this work, we focus on single component fluids.To study multicomponent systems, we would need to include the composition as an additional variable and consider the many ways that the composition couples to the effects that are already present for pure fluids.Likewise, we have not considered electrical effects, which are so important in microfluidic and lab on a chip systems.These topics have been discussed in detail by other authors [2][3][4][5], and we urge interested readers to consult their work.We assume that all flows discussed are time independent, so we restrict our attention to steady states, but it is worth mentioning that nanoscale viscoelasticity remains largely unexplored.

Slip
In most macroscopic flow situations, a stick boundary condition, where the fluid velocity at the wall is taken to be equal to the wall velocity, is assumed.At the macroscopic scale, it is only in extreme cases that we must allow for slip, for example when we have plug flow of a paste or polymer melt, or when trapped gas bubbles in superhydrophobic surfaces lead to extreme slip [6].However, even simple liquids can experience some slip.Strong slip is expected for water flowing near hydrophobic surfaces such as graphene or carbon nanotubes.Given a suitable surface structure, slip can also be observed for relatively hydrophilic surfaces [7].For Poiseuille flow through a wide channel, the slip velocity, defined as the difference between the wall velocity and the fluid velocity at the wall, is only a very small fraction of the maximum velocity of the fluid in the channel, and so it is safely approximated as zero.On the other hand, for very narrow channels of nanometre scale, the slip velocity can be a significant fraction of the maximum velocity in the channel.Predictions of flow rates based on the assumption of the no-slip boundary condition can therefore be significantly in error, underestimating the true flow rate.
Here, we focus on the slip that occurs at the atomic scale on atomically smooth surfaces, sometimes known as intrinsic slip.Slip that occurs on a larger length scale, involving structured or patterned surfaces, roughness and chemical heterogeneity has been discussed by other authors [8].
The slip velocity depends on the strain rate at the wall and is not a material property.The relevant material property describing slip is the slip friction coefficient, defined below.The ratio of the fluid viscosity to the slip friction coefficient gives us another material property (really a property of the interface between the two materials), called the slip length.Very high spatial resolution measurements of the flow velocity of aqueous solutions near a hydrophobic wall, for example, have shown slip lengths of 80-100 nm [9].The slip length of water confined by highly hydrophobic graphene and carbon nanotube surfaces has been computed to be around 60 nm, but there is enormous variation in both experimental and simulation results [10], due partly to subtle differences in simulation technique [11] and molecular models, but also error prone data analysis [12] and experimental difficulties.The consensus of careful simulation and experimental studies of flow of water on molecularly smooth hydrophobic surfaces is that slip lengths on these surfaces typically vary from nanometres up to tens of nanometres [10].Because it has been difficult to reproduce the much larger experimental values sometimes found, it is possible that measurements of some of the larger values could suffer from uncontrolled experimental factors, such as dissolved or trapped gases [13], surface roughness and imperfections [14] or other yet unidentified factors.
One of the difficulties with computation of the slip length and slip friction coefficient by direct evaluation in non-equilibrium molecular dynamics that simulate flow through a channel is thatk for high slip systems, the velocity profile is almost flat.Extrapolation of such a velocity profile to the position where the velocity is zero is extremely error prone [12].Methods based on equilibrium correlation functions do not suffer from these problems.Several theoretical discussions of correlation function methods for computation of the slip friction coefficient have been published [15][16][17][18][19][20] but there is still some doubt concerning the agreement between the different forms.Here, we describe a simple one [18] that has been extensively validated by comparison with both nonequilibrium simulations and experimental results [10,12,[21][22][23].
Navier's slip friction coefficient ξ is defined as the proportionality coefficient relating the shear pressure P yx at the wall to the velocity difference between the wall and the fluid ∆v x , which we call the slip velocity where we assume that the flow is in the x direction and the velocity gradient is in the y direction.
The flow geometry for flow in the x direction adjacent to a flat solid wall is shown in Figure 1.The viscous pressure in the fluid at the wall is given by Newton's law of viscosity where η is the shear viscosity at zero shear rate.Since the shear pressure must be continuous, both values of the shear pressure must be equal.Eliminating the shear pressure, we find which defines the slip length as To derive a correlation function expression for the slip friction coefficient, we must introduce a generalised constitutive equation that describes the fluctuations [18].It is convenient to formulate this in terms of the wall-fluid shear force, given by where A is the area over which the slip frictional force acts, ζ (t − τ) is the friction kernel and F R is the random component of the shear pressure.Now, we multiply both sides of the constitutive equation by the slip velocity and ensemble average to form the correlation functions where When this equation is Laplace transformed, we find The friction kernel is well approximated by a sum of exponentials [18], which has the Laplace transform In practice, a single exponential is often sufficient [18,21].The amplitudes and decay rates of the exponentials can then be obtained from fits to the Laplace transformed correlation functions.
From the computational point of view, the wall-fluid shear force is easily evaluated in computer simulations.The slip velocity is more problematic.One way to evaluate it would be to fit the instantaneous velocity profile each time the correlation functions are evaluated.However, this requires the assumption of a fitting function, which could be biassed.Another way is to evaluate the instantaneous velocity of the fluid averaged over a region within a distance ∆ of the wall.This has the disadvantage that averaging over a finite region could also result in error, but the calculation of the slip friction coefficient for a given wall-fluid combination can be repeated for different values of ∆, and the most physically meaningful value of the slip friction coefficient chosen.In practice, it is straightforward to choose the most physically meaningful value of the slip friction coefficient because it quickly increases to a broad maximum (plateau) value as ∆ is increased, before steadily decreasing thereafter.This maximum value usually occurs when ∆ is approximately equal to one molecular diameter.Choosing this value of the slip friction coefficient gives excellent agreement with the results of nonequilibrium molecular dynamics simulations [18,21].
This method has been used to evaluate the slip friction coefficient and the slip length of a simple Lennard-Jones type atomic fluid near a solid planar LJ wall [18] and a graphene wall [21] and more complex molecular fluids such as water against both Lennard-Jones atomic walls and graphene [12].It has also been adapted to a cylindrical geometry [22] for studies of the slip friction coefficient of water in carbon nanotubes [10,23].

Spin Angular Momentum Coupling
It is well known that extended molecules spin in a shear field [24].What is less well known is that the spin angular velocity is coupled to the translational velocity.This coupling is usually negligible in macroscopic flows, but it can become significant at the nanoscale.
To describe this effect, we begin with the extended Navier-Stokes equations for stable flow in a planar flow geometry (where the convective terms are zero), which are obtained by inserting the relevant linear transport equations into the balance equations for the translational and angular momentum [25,26], The first of these describes the evolution of the translational fluid streaming velocity in the flow direction v x .A new transport coefficient, the rotational viscosity η r , which governs the relaxation of spin angular momentum appears in addition to the usual shear viscosity, η.The rotational (or vortex) viscosity η r describes the rate of conversion between fluid vorticity ∇ × v and molecular spin angular momentum ω.If η r = 0, then we just have the usual Navier-Stokes equation.In the absence of a pressure gradient (omitted from this equation since it is zero for field driven flows in our flow geometry), flow can be generated by the external body force density ρF e , but if this is also zero, we see that it is also possible to have a flow driven by the gradient of the angular velocity.We will return to this point later.The second equation, which describes the evolution of the molecular spin angular velocity field includes a diffusive term involving the sum of the spin viscosities ζ + ζ rr as well as the rotational viscosity η r .The spin viscosity ζ describes the diffusive flux of spin angular momentum due to the traceless symmetric part of the spin angular velocity gradient, while ζ rr describes the diffusive flux of spin angular momentum due to the antisymmetric part of the spin angular velocity gradient [25].These transport coefficients have been evaluated for some molecular fluids, including water [27].Again, there is a term that accounts for external fields, this time in the form of an external body torque density ρΓ e z .When applied to planar Poiseuille flow through a narrow channel driven by an external body force density ρF e , these equations predict a difference between the velocity field calculated from the Navier-Stokes equation alone compared to the results of the extended Navier-Stokes equations including the spin coupling.Simulations of a molecular fluid consisting of extended linear molecules (buta-triene) show that the flow rate difference is small for channels of width greater than 7 nm, but it grows to around 10% at a channel width of 1 nm [28].For very narrow channels, accurate prediction of flow rates and velocity profiles requires consideration of the extended Navier-Stokes equations [26,28,29].Figure 2 shows the flow rate reduction predicted by taking spin angular momentum coupling into account for nanoconfined flows of a dumbbell fluid, liquid butane and liquid water [26].Copyright 2015 American Chemical Society [26].
As mentioned above, a translational flow can be generated even in the absence of a translational body force if a body torque is applied instead.With symmetric boundary conditions (for example, with a stick boundary condition on both sides of the channel), equal flow is generated in both directions and no net flow results, but if the boundary conditions are asymmetric, with slip on one wall and stick on the other, a net flow results.This means that an external torque that spins the molecules can be used to pump a fluid.It has been demonstrated that a rotating electric field applied to polar molecules (such as water) under these conditions can generate a net flow in a nanochannel or nanotube, without the need for electrolyte or a pressure gradient [30,31].

Nonlocal Response
In macroscopic flows, the most common causes of non-Newtonian behaviour are nonlinear and viscoelastic deviations from Navier-Stokes behaviour.The Pipkin diagram [32] shown in Figure 3 schematically illustrates the different regions of fluid behaviour for typical macroscopic fluids.Here, we are interested in steady (zero frequency) flows so the Deborah number ωτ v is zero.In most macroscopic situations, spatial nonlocality is unimportant, except possibly when the system is near a glass transition [33] or the viscous correlation length is extraordinarily large for some other reason as it is for example in suspensions of long fibres that exhibit shear banding [34,35].By contrast, in nanofluidic flows, spatial nonlocality can be a dominant effect that must be taken into account, even for ordinary molecular or atomic fluids.Deformations at sufficiently high wavenumber may differ strongly from homogeneous deformation, particularly in glassy liquids where dynamic heterogeneity with nanoscale dimensions is observed [33].This means we must consider the concept of a spatially nonlocal response, where the stress at a point becomes a linear functional of the local deformation rate [36].The stress at a point then depends not only on the strain rate at that point but also on the strain rate at nearby points.An alternative interpretation is that the stress depends not only on the velocity gradient but also on its derivatives, and so there is a contribution to the stress resulting from the spatial derivatives of the strain rate.
To describe spatially inhomogeneous steady state flows where the strain rate is independent of time but it varies strongly in space, we can introduce an analogue of the Pipkin diagram to characterise nonlocality, as shown in Figure 4.The viscous correlation length of the fluid, ξ v is measured by the width of the nonlocal viscosity kernel, η (y − y ) defined in Equation (12).For situations where strong density inhomogeneity is also present (discussed below in this section), third and fourth axes could be added, representing the amplitude and spatial frequency of the density inhomogeneity.This diagram does not include these dimensions, so it represents nonlocal response for a uniform density fluid.To account for the nonlocal response of the shear pressure to a steady velocity gradient of high spatial wavenumber, we apply a spatial convolution equation analogous to the linear viscoelastic constitutive equation This is the most general linear relationship that we can postulate between the shear pressure and the velocity gradient or shear rate γ for a fluid with spatially homogeneous density where the strain rate varies only in the y-direction.The viscosity kernel η (y − y ) weights the contribution of the strain rate at different distances from the point at which the shear pressure is evaluated.When y − y = 0, we expect the contribution of the shear rate to be greatest, and at larger values of y − y , we expect the effect of the shear rate to be least, so η should be a decreasing function of y − y .At macroscopic length scales, it is reasonable to approximate η (y − y ) as a constant multiplied by a Dirac delta function, and then the convolution integral simply reduces to Newton's law of viscosity with a purely local response to the strain rate field.
A stringent test of this relationship can be made by applying a spatially sinusoidal transverse force to generate a sinusoidal velocity field.Since the velocity is sinusoidal, if the amplitude of the field is sufficiently small, the response is linear and the shear pressure response will also consist of a single sinusoidal component.When the wavelength of the sinusoidal driving force is sufficiently long (greater than a few molecular diameters) the shear pressure is just given by the Newtonian constitutive equation applied locally.In other words, the shear pressure at a given point depends only on the value of the strain rate at that point.However, when the wavelength of the strain rate oscillations is reduced to the order of a few molecular diameters, this procedure fails and the shear pressure is poorly reproduced.Using the nonlocal integral constitutive equation, it is possible to correctly predict the shear pressure, even when the velocity profile varies rapidly in space as displayed in Figure 5 [36].
In nanofluidic systems, it is possible for the strain rate to vary rapidly, and for a nonlocal response of the shear pressure to the rapidly varying strain rate to be observed.However, the rapid variation of the strain rate is also usually associated with rapid variation of the density of the fluid, due to strong molecular packing effects at the wall-fluid interface that occur even at equilibrium [37].Therefore, we must consider the effect of spatial density variations on the velocity profile and the shear pressure profile. .Shear stress yx ) predictions using the simple Newtonian constitutive equation (dashed line), and the nonlocal constitutive equation (solid line) compared with the exact stress (filled circles) for a simple liquid [36] where the y-position is expressed in units of the Lennard-Jones potential distance parameter σ.

Density Inhomogeneity
When a macroscopic fluid system has long wavelength density inhomogeneities, the assumption of constant transport coefficients in the Navier-Stokes equation is clearly inadequate.We can allow for this by making the linear transport coefficients position dependent through the position dependence of the density as well as, if necessary, the temperature and composition.In the simple case of density variation for an isothermal single component fluid, we would write Newton's law of viscosity as In a nanofluidic system, the fluid is almost always in contact with a solid wall.This induces strong oscillatory density variations near the wall, with local maxima that may exceed typical solid densities.Under these circumstances, it is obviously not viable to use the viscosity at the local density values.Bitsanis and coworkers [38,39] made a significant improvement on the local density model by proposing that the viscosity at a locally averaged value of the density could be used in the Newtonian viscosity equation.Despite the strong density variations near the wall, averaging over a sphere or planar layer of one to two molecular diameters gives an average density that is reasonably close to the liquid value, but still accounts for slow variations in the density.The constitutive equation for the local average density model in a planar geometry is then This introduces an element of nonlocality into the constitutive equation for the shear pressure by making it dependent on the density averaged over a region surrounding the point of interest.The local averaged density model was successfully used by Bitsanis [38] to describe the velocity profile of a highly confined simple fluid.It has been extended by Hoang and Galliero [40][41][42], who investigated the effectiveness of different averaging kernels in the evaluation of the averaged density, and combined with a model for slip to provide a tractable and efficient hydrodynamic model for nanoflows by Bhadauria et al. [43,44].The local average density model is, however, only a partial solution to the problem of describing the velocity profile for flow of a nanoconfined fluid.A major deficiency of this model is that it cannot account for the zeroes and velocity gradient reversals that can be seen in the velocity profiles of strongly confined fluids at high densities [45].Strong molecular packing near the solid-fluid interface results in oscillations in the velocity profile that cannot possibly be described by the local average density model.Any constitutive equation that follows the same functional form as the simple Newtonian one would need to have infinite values of the viscosity at the points where the strain rate goes through zero, which is clearly unphysical.Therefore, we are again led to consider more general, nonlocal constitutive equations.
The fluid density inhomogeneity produced by wall-fluid interactions is uncontrolled.It is not possible to easily and independently control the amplitude and spatial frequency of the density oscillations due to the presence of a confining wall.To study nonlocal constitutive equations for the shear pressure, we need a more flexible framework.By applying a small amplitude sinusoidal longitudinal force (SLF) in a periodic simulation cell, it is possible to induce sinusoidal density variations [46] in a system with periodic boundaries and no explicit walls.When the amplitude of the applied force is increased, a nonlinear density response is generated.By adding together several sinusoidal forces of different frequency, it is possible to Fourier synthesise a spatially periodic density profile consistent with the periodic boundary conditions of the simulation cell that closely resembles the density profile seen in a nanoconfined fluid [47].A sinusoidal transverse force (STF) can also be applied that generates flow.When both the longitudinal and transverse fields are applied together, we can generate flow in the presence of density inhomogeneity that can be controlled and manipulated at will [48,49].
Considering a system where the external fields vary in the y direction, the combined longitudinal and transverse force can be represented as Figure 6 schematically shows the two applied body forces.The density response consists of contributions due to both the longitudinal and transverse forces.This is because the transverse force results in flow with a strain rate that varies in the y-direction.Heat is produced by the viscous dissipation associated with the velocity gradient, resulting in thermal expansion, which changes the density, in addition to the direct effect of the longitudinal force on the density.In general, the density response is given by a functional expansion that depends on both the longitudinal and transverse forces, where α 1 and α 2 can be either x (transverse force) or y (longitudinal force) and the response functions are the functional derivatives The response functions are evaluated at equilibrium.Due to the spatial symmetry of the equilibrium system, the response functions must depend only on even powers of the transverse force, because changing its direction cannot change the sign of its contribution to the density.We assume that truncation of the density response at second order in the forces gives a reasonable approximation.It accounts for the density response due to the longitudinal force with terms that are first and second order in the longitudinal field and also the heating and normal pressure effects of the shearing force, which occur to lowest order as a quadratic function of the transverse field, ρ(y) =ρ 0 + χ By expressing both sides of this equation in terms of their Fourier series representations, it is possible to determine the Fourier coefficients of the response functions.
Applying similar arguments to the truncated functional expansion of the strain rate, we find γ(y) = ξ x (y − y ) F x (y ) dy + ξ This expression has been limited to terms that are linear in the transverse force and at most quadratic in the longitudinal force since we are mainly interested in cases where the shear is weak but the density inhomogeneity is strong.The shear pressure profile can be written in a similar form, as x (y − y ) F x (y ) dy + π (2) xy (y − y , y − y ) F x (y )F y (y ) dy dy + π (3) xyy (y − y , y − y , y − y ) F x (y )F y (y )F y (y ) dy dy dy , (21) With the addition of the transverse force, we can generate flow in a system with strong, and realistically complicated density inhomogeneity.Figure 8 shows the resulting velocity profiles with the predictions obtained from the truncated functional expansion for the velocity gradient.This method is clearly capable of accurately representing the response of a nanoscopic liquid system to a combination of confinement and flow-inducing forces.Copyright (2015) by the American Physical Society [49].

Conclusions
In this paper, we have provided a brief review of three fundamental phenomena that should be included in an accurate continuum treatment of nanofluidics: slip, spin angular momentum coupling, and non-local response.At the outset, we pointed out that additional complexities would need to be included to account for compositional variation (in binary and multicomponent fluids) and electrostatics, such that a full treatment of ionic solutions and liquids or even polar fluids can be treated.This is clearly some years away.The main complexity is in dealing with the full range of couplings and non-local effects that occur at the nanoscale.A major challenge remains in how to simplify the treatments given in References [46][47][48][49][50] and potentially others not discussed in this brief review such that they can be efficiently applied in predictive nanofluidics.Further advances in the modelling of slip will undoubtedly see the application of both NEMD and EMD methods to complex fluids (e.g., multicomponent mixtures and ionic liquids) and more targeted engineering applications, such as lubrication (see, for example, [51]).The coupling of linear and angular momentum to generate flow poses an intriguing potential application in fluid actuation at the nanoscale.While theoretical and simulation studies exist, there has been no experimental verification to our knowledge.

Figure 1 .
Figure 1.Schematic diagram showing the flow geometry for the definition of the slip length.The magnitude of the velocity gradient in the fluid at the wall is equal to v s /L s where L s is the magnitude of the slip length.Here, we have ∆v x = v x − 0 because we assume flow between stationary walls.

Figure 3 .
Figure 3. Diagram showing rheological flow regimes for typical macroscopic fluids.The horizontal axis represents the degree of elasticity exhibited by the flow, which is controlled by the Deborah number, the product of the characteristic frequency of the flow ω and the viscoelastic relaxation time of the fluid τ v .The vertical axis represents the degree of nonlinearity of the flow, controlled by the Weissenberg number, the product of the characteristic strain rate γ and the viscoelastic relaxation time τ v .

Figure 4 .
Figure 4. Diagram showing the rheological flow regimes for fluids exhibiting nonlocal shear pressure response under steady flow conditions.The horizontal axis represents the degree of nonlocality exhibited by the flow, which is controlled by the product of the characteristic spatial frequency (or wavenumber) of the flow k and the viscous correlation length of the fluid ξ v .The vertical axis represents the degree of nonlinearity of the flow, controlled by the product of the characteristic strain rate γ and the viscoelastic relaxation time τ v .

Figure 5
Figure5.Shear stress yx ) predictions using the simple Newtonian constitutive equation (dashed line), and the nonlocal constitutive equation (solid line) compared with the exact stress (filled circles) for a simple liquid[36] where the y-position is expressed in units of the Lennard-Jones potential distance parameter σ.

Figure 6 .
Figure 6.Schematic of: the sinusoidal transverse force (STF) (a); and the sinusoidal longitudinal force (SLF) (b) fields.The arrows show the direction of the forces.The length of the arrows, as well as the sinusoidal line, indicate the strength of the force.The STF is shown for F x (y) = F x 1 sin(k 1 y) and the SLF is shown for F y (y) = F y 2 sin(k 2 y).Reprinted figure with permission from Glavatskiy, K.S.; Dalton, B.A.; Daivis, P.J.; Todd, B.D., Phys.Rev. E 91, 062132 (2015).Copyright (2015) by the American Physical Society [48].

Figure 8 .
Figure 8.Comparison between velocity profiles obtained directly from MD simulations and predictions of the truncated functional expansion using independently calculated response functions for a single sinusoidal component STF and a three component SLF.Bold lines represent MD simulations, while thin dashed lines with symbols represent predictions.The system parameters labels are the same as those used in Figure 7. Velocity profiles are only shown for half of a wave cycle.Reprinted figure with permission from Dalton, B.A.; Glavatskiy, K.S.; Daivis, P.J.; Todd, B.D. Phys.Rev. E 92, 012108 (2015).Copyright (2015) by the American Physical Society[49].