On the validity of the Navier-Stokes equations for nanoscale liquid flows: The role of channel size

In this work, we investigate the validity of the Navier-Stokes (NS) equations for nanoscale liquid ﬂows through molecular dynamics simulations. We focus on the role of channel size by considering the ﬂuid-wall interaction. Liquid ﬂows between two planar parallel walls driven by an external force with channel size ranging from 2 to 80 nm are studied. The volumetric ﬂux is computed and the dependence of the volumetric ﬂux on the channel size is explained both qualitatively and quantitatively. It is found that the ﬂow is sensitive to the ﬂuid-wall binding energy and the classical ﬂuid mechanics falls apart in small nanochannels. However, the wall effects become insigniﬁcant and the NS equations are valid when the channel size is larger than about 150 molecular diameters ( ∼ 50 nm). Copyright 2011 Author(s). This article is distributed under a Creative Commons Attribution 3.0 Unported License .


I. INTRODUCTION
Nanoscale liquid flows are of great importance in both science and engineering. [1][2][3] They have attracted much attention in the past decade due to their potential applications in a variety of areas, including molecular separation and chip-level cooling. 1,4 A deep understanding of nanochannel flows is also critical for the further miniaturization of microfluidic devices. In small nanochannels, due to the strong confinement, the wall may affect the flow greatly through some microscopic parameters, such as the fluid-wall molecular binding energy, which are not considered in the classical fluid mechanics. 4,5 These microscopic parameters usually couple with macroscopic parameters and lead to new flow phenomena. For instance, in small nanochannels, the fluid density and viscosity can be inhomogeneous [6][7][8] and the boundary condition at the fluid-wall interface can vary from slip to multilayer sticking with the change of fluid-wall binding energy. [9][10][11][12][13][14][15][16][17][18] As a consequence, the classical Navier-Stokes (NS) equations may not be valid under strong confinements. 19,20 On the other hand, when the channel size is sufficiently large, the wall effects are negligible and the flow behavior can be well described by the NS equations. Therefore, as the channel size is increased, there exists a transition from the regime, where the wall effects are important, to the continuum regime, where the classical fluid mechanics works. Experimentally, it has been observed that liquids show noncontinuum behaviors in 5 nm-diameter carbon nanotubs (CNTs), while the flow is sort of continuum in CNTs of 50 nm in diameter. 21,22 However, at what channel size the microscopic parameters are essential and the classical NS equations become invalid remains unanswered.
In recent years, the influence of channel size on nanoscale liquid flows has been investigated. [23][24][25][26] Although many interesting flow fashions have been observed, the findings about the role of channel size and the validity of the NS equations for nanoscale flows seem to be controversial. Some work concludes that the flow motion can be predicted by the NS equations even when the channel size is of 10 molecular diameters (∼ 3 nm), 19 while some other work finds that the flow shows strong dependence on the fluid-wall binding energy in channels of width about 20 molecular diameters. 23 The confusion about the validity of the classical fluid mechanics is a consequence of the unavailability 2158-3226/2011/1(3)/032108/10 C Author(s) 2011 1, 032108-1 of the channel size, where the NS equations fall apart. In most of the previous studies, the channel size is smaller than 20 molecular diameters (∼ 6 nm), where the confinement is still strong and the flow strongly depends on the fluid-wall interactions. In small nanochannels, molecular dynamics (MD) simulations have shown that the NS equations appear to be valid in estimating the flow rate if the fluid-wall binding energy ε fw is well controlled. However, if ε fw is varied, the flow rate can dramatically deviate from the prediction of the NS equations. 27 The results of our recent studies also indicate that nanochannel flows are affected by ε fw and the influence of ε fw weakens with increasing channel size. 4,5 Therefore, the controversial results about the validity of the NS equations, to some extent, confirm that there is a transition in nanochannel flows and there exists a regime-division channel size, above which the NS equations apply no matter how the fluid-wall interaction is changed. Unfortunately, this channel size is unknown.
In this work, through MD simulations, we investigate the regime-division channel size by studying liquid flows in nanochannels with a wide range of channel sizes under different fluid-wall binding energies ε fw . We focus on Poiseuille flows, where the fluid is driven by an external force. The transition of the flow from ε fw dependent to the flow regime where the classical fluid mechanics works is observed and it is found that the regime-division channel size, above which the NS equations are valid regardless of fluid-wall binding, is about 150 molecular diameters (∼ 50 nm).

II. SIMULATION METHOD
The Poiseuille flow systems consist of two planar parallel solid walls with a liquid confined in between. A typical system is shown in the inset of Fig. 1. The initial positions of the wall atoms are obtained by truncating a rectangular portion from a face-centered cubic (fcc) structure with lattice constant equal to 4.086 Å. The tight-binding potential function is used to calculate the interactions among the wall atoms such that the thermal motion of the wall is considered, which affects the energy accommodation of the wall and the dynamics of the fluid molecules at the fluidwall interface. 28,29 The parameters for Ag are used to simulate the wall. 30 The fluid is modeled by the Lennard-Jones (LJ) potential, U (r ) = 4ε ff [(σ ff /r ) 12 − (σ ff /r ) 6 ], where r is the separation between a pair of interacting fluid molecules, ε ff is the fluid-fluid self-binding energy, and σ ff is the collision diameter of the fluid molecules. The values of ε ff and σ ff are set to be 114 K and 3.47 Å, which can be used to describe liquid Ar. 31 The reason why liquid Ar is used is because the binding energy of Ar is ∼ 100 K, which can represent many popular liquids, such as water. 31,32 The fluid-wall interaction is also calculated by the LJ potential and the Lorentz-Berthelot mixing rule is used to calculate the interaction parameters. In order to implement the mixing rule for the interaction, the binding energy ε ww = 3995.4 K and collision diameter σ ww = 2.54 Å are used for the wall (Ag) atoms. 33 In the simulations, to explore the wall effects on the flow, the fluid-wall binding energy ε fw is arbitrarily varied in a wide range, which can cover most of the liquid-solid interactions. 31,34 The planar walls are perpendicular to the y axis. Four layers of atoms are used to model the walls (Fig. 1). The atoms in the outmost layer are fixed to maintain a stable system whereas the atoms in the other three layers are free to vibrate to consider the flexibility of the wall. In order to consider the channel size effect, the distance between the two innermost layers is varied from 2 to 80 nm, and the lengths of the system are 19.5 and 4.9 nm in the x and z directions respectively. Liquid molecules are also initially on fcc lattice sites and uniformly distributed in the channel. The initial velocities of the wall atoms and fluid molecules follow the Maxwell-Boltzmann distribution and correspond to the desired temperature. Periodic boundary conditions are employed in the x and z directions only. An external force is applied to each fluid molecule in the positive x direction. The Berendsen thermostat is used to maintain the temperature of the wall at the desired value. 35 The potential is truncated at 10.21 Å and Newton's equations are integrated with Beeman's leapfrog algorithm. 36 The time step is equal to 1 fs. At the very beginning, the system is relaxed for 100 ps to reach equilibrium. After the equilibration, an external force is applied to each fluid molecule. This external force is linearly increased with time until it reaches the desired value in 50 ps. Afterward, the system is allowed to equilibrate for another 3 ns before the flow properties are calculated.

III. RESULTS AND DISCUSSIONS
In the simulations, the channel size is varied in a wide range. For large channels, the temperature of the fluid can be much higher than the thermostat applied to the wall due to viscous frictions. 37,38 To maintain a uniform temperature in the whole system, the external force applied on the fluid molecules are carefully controlled. In nanoscale flows, fluid adsorption on solid surfaces is a popular phenomenon, which leads to new flow behaviors. 5 The fluid adsorption is favored by strong fluidwall binding energy ε fw and low temperature T . 5, 28 Therefore the dimensionless number ε fw /kT is usually used to characterize the wall effect, where k is the Boltzmann constant. In this work, the temperature T is maintained at 300 K. To understand the wall effects, ε fw /kT is varied by changing the fluid-wall binding energy ε fw from 30 to 3000 K, which can cover the interactions between most fluids and solids.
The volumetric flux Q of the liquid is chosen as the major flow property for comparing MD simulations and the NS equations. In nanochannels, the fluid properties, such as density and viscosity, may fluctuate across the channel and their effects on the flow motion are sort of considered if the flux Q is used. The velocity and density profiles of the liquid are also used for the comparison. Figure 1 shows the volumetric flux ratio Q/Q NS of the liquid as a function of the scaled channel size D/σ ff under different ε fw /kT , where Q NS is the volumetric flux predicted by the NS equations and D is the channel size. It is seen that the flux strongly depends on the fluid-wall binding energy ε fw /kT and deviates greatly from unity in small channels. As the channel size is increased, it converges to 1 no matter how ε fw /kT is changed. It is also seen that the flux is always larger/smaller than unity for weak/strong fluid-wall interactions. Figure 1 demonstrates that the flow is sensitive to the wall effects when the channel is small, where whether the NS equations apply or not depends on the fluid-wall interaction. This explains the controversial results in the literature. Figure 1 also indicates that the NS equations can be viewed as definitely valid when the channel size is larger than 150 molecular diameters (∼ 50 nm) if the relative error of 5% is acceptable for practical applications.
To qualitatively understand the fashion of the volumetric flux in Fig. 1, it is necessary to obtain the velocity V and density ρ profiles of the fluid in the channels, as shown in Figs. 2 and 3. In Fig. 2, it is seen that the velocity profiles are generally parabolic in relatively large channels. However, their dependence on the fluid-wall binding energy is different in different channels. In small nanochannels (e.g. 4 nm channel in Fig. 2(a)), they are sensitive to the fluid-wall binding energy. For weak fluidwall interactions (small ε fw /kT ), due to small interfacial friction, the viscous resistance becomes unimportant and there is a clear velocity slip at the interface, which enhances the flow flux. The velocity slip is reduced as the fluid-wall interaction grows strong and the flow resistance increases. For sufficiently large ε fw /kT , a layer or more liquid molecules are adsorbed by the wall, as indicated by the velocity distribution in Fig. 2(a) for ε fw /kT = 10 and the density profiles in Fig. 3. With the adsorption, the effective channel size available for the flow decreases and the flow flux drops, as will be discussed later.
As the channel size is increased, the velocity distribution becomes less sensitive to the fluid-wall interaction and converges to the prediction of the NS equations, as confirmed by the velocity profiles in Fig. 2(f). This is because the fluid-wall interaction is short-range and the wall can only affect the liquid molecules very close to the wall. Figure 4 depicts the potential distribution in the vicinity of the wall for different values of ε fw /kT . It is found that the potential changes greatly near the wall as ε fw /kT is increased by a factor of 100, from 0.1 to 10. However, the potential distribution does not change too much a couple of molecular layers away from the interface. Therefore, the flow field in the channel can be divided into two regions, i.e., the parameter sensitive region near the wall and the main stream region, where the wall effect is insignificant. In the parameter sensitive region, the potential wells deepen with increasing fluid-wall binding energy and fluid molecules are apt to be trapped by the wall, which will increase the flow resistance and reduce the flow flux. In the main stream region, however, the wall barely affects the fluid and the flow motion is dominated by the fluid-fluid interactions (viscosity) and can be well described by the NS equations on the basis of the fact that the fluid density in this region is uniform (Fig. 3). The volumetric ratio of the two regions can be used to qualitatively explain the dependence of the flux on the channel size in Fig. 1. In small channels, the volume of the parameter sensitive region is significant compared with that of the main stream region and therefore the flow depends on the fluid-wall interaction. As the channel size increases, the percentage of the main stream region increases while that of the parameter sensitive region decreases. Consequently, the fluid-fluid viscous interaction tends to be dominant and the NS equations become valid.
Since the velocity profiles are generally parabolic, a quantitative analysis can be performed based on the flow velocity for small ε fw /kT and effective channel size for large ε fw /kT by considering the velocity slip and the size of parameter sensitive region, respectively. For weak fluid-wall interactions, as shown in the inset of Fig. 5(a), the volumetric flux can be estimated by the parabolic velocity and the velocity slip as Q = V slip +V NS , where V slip is the slip velocity, andV NS is the mean velocity (volumetric flux) predicted by the NS equations. Thus, the volumetric flux ratio reads, Therefore, the discrepancy between the real volumetric flux and the prediction of the NS equations is determined by the velocity ratio V slip /V NS . SinceV NS scales with D 2 ,V NS ∝ D 2 , as the channel size increases, V slip /V NS approaches zero given that V slip is weakly dependent on D compared with V NS , 39 as shown in Fig. 5(a). It is noted that V slip depends on the competition between the fluid-wall binding force and the external driving force. Even if these conditions are in favor of a large velocity slip, the viscous friction at the interface can heat up the fluid and stop V slip from increasing. Hence, there is a self-control mechanism in nanoscale flows to maintain V slip as finite. 37,39 For strong fluid-wall interactions, the velocity slip becomes negative due to fluid adsorption on the wall. In this case, the effective channel size D e offers a better way to evaluate Q/Q NS by considering the area of the parameter sensitive region near the wall. The effective channel size can be defined as D e = D − 2λ, where λ is the thickness of the parameter sensitive region (adsorbed liquid layer), as depicted in the inset of Fig. 5(a). Within the effective channel, it is seen that the velocity distribution can be well fitted by a parabolic profile such that Q/Q NS can be calculated as where A is the cross-sectional area of the channel and the subscript "e" denotes the effective channel of channel size D e . Since the thickness of the parameter sensitive region λ is finite and independent of channel size D, Q/Q NS (D e /D) changes from a value smaller than 1 and converges to unity as the channel size increases, as shown in Fig. 5. It is worth mentioning that λ depends on ε fw /kT . However, λ varies only in a short range from zero to a few angstroms even if ε fw /kT is changed by a factor of 100. This information can help offer a conservative estimation of the channel size where the classical fluid mechanics is valid based on Eq. (2). Equations (1) and (2) also explain why Q/Q NS is always positive/negative for weak/strong fluid-wall interaction as the channel size varies. In Fig. 5(b), the predictions of Eqs. (1) and (2) are compared with MD simulations. It is seen that they are in good agreement when the channel size is larger than 12 molecular diameters (∼ 4 nm). For very small channels, the discrepancy is caused by the large numerical error due to the fluctuations in fluid density and non-parabolic velocity distributions. 40,41 Alternatively, the fluid density profiles can also be used to understand the validity of NS equations. As seen in Fig. 3, the density fluctuates greatly in the area near the wall, while it is quite uniform in the center area of the channels. As the channel size increases, the density in the center area ρ c approaches the average density ρ o and the NS equations tend to work. In Fig. 6(a), the density ratio ρ c /ρ o is depicted as a function of D/σ ff for different values of ε fw /kT . It is seen that ρ c /ρ o → 1 as D increases and the deviation of ρ c from ρ o is below 2% when the channel size is larger than 50 nm.
It is noted that the average fluid density was maintained constant (ρ 0 = 21.14/nm 3 ) in the above studies. To understand the effect of fluid density on the channel size (50 nm), above which the NS equations apply, we examine the density distribution in a 60 nm channel, as shown in Fig. 6(b). It is found that the fashion of density fluctuation is similar to that in Fig. 3. The density varies near the wall and becomes quite uniform several molecular layers away from the wall. Furthermore, ρ c is not significantly changed even if ρ o is reduced by about 50%. Therefore, it is expected, based on the density ratio shown in Fig. 6(a), that the fluid density ρ o does not greatly affect the general rule of 50 nm obtained from Fig. 1.

IV. SUMMARY
In summary, the role of channel size in the validity of the NS equations has been investigated through MD simulations by considering the wall effects. It has been shown that the fluid-wall interaction plays an important role in nanoscale flows when the channel size is small. However, the