Continuum simulations of water flow in carbon nanotube membranes

We propose the use of the Navier–Stokes equations subject to partial-slip boundary conditions to simulate water flows in Carbon NanoTube (CNT) membranes. The finite volume discretizations of the Navier–Stokes equations are combined with slip lengths extracted from molecular dynamics (MD) simulations to predict the pressure losses at the CNT entrance as well as the enhancement of the flow rate in the CNT. The flow quantities calculated from the present hybrid approach are in excellent agreement with pure MD results while they are obtained at a fraction of the computational cost. The method enables simulations of system sizes and times well beyond the present capabilities of MD simulations. Our simulations provide an asymptotic flow rate enhancement and indicate that the pressure losses at the CNT ends can be reduced by reducing their curvature. More importantly, our results suggest that flows at nanoscale channels can be described by continuum solvers with proper boundary conditions that reflect the molecular interactions of the liquid with the walls of the nanochannel.


Introduction
In recent years, detailed experiments and simulations have demonstrated that flows in natural and artificial nanochannels exhibit remarkably high flow rates [1][2][3][4][5][6]. These findings have spurred great interest in applications ranging from nanosyringes for drug delivery to new technologies for water desalination [7,8].
The two main lines of investigation for these flows involve experiments and simulations using molecular dynamics (MD). These two approaches are considered complementary, however, there are large discrepancies in the spatiotemporal scales that each can resolve. These discrepancies have fueled arguments with regards to the large differences observed in flow enhancement rates calculated in experiments and MD simulations [4,9,10].
The computational cost of the MD simulations and the relatively short time scales that they resolve has motivated the development of a continuum model [11] with partial-slip boundary conditions to model fluid flows through carbon nanotube (CNT) membranes. Such partial slip boundary conditions are consistent with the observation that at the nanoscale the wall-fluid interaction plays an important role on the overall behavior of the flow and the standard empirical law of the no-slip boundary condition is not applicable [12].

Governing equations
In this work, we present a continuum, numerical study of water flow through CNT membranes using finite volume discretization of the steady-state incompressible Navier-Stokes equations: 2 subject to partial-slip boundary condition at the fluid-wall interface: where v t is the tangential to the wall component of the fluid velocity and ∂ ∂ v n t its gradient in the direction normal to the wall [13], r is the radius of the curvature of the wall in the flow direction with the radius being negative for concave boundaries and l s is the slip length. The slip length is defined as the ratio between the fluid viscosity and the friction coefficient between the wall and the fluid and can be extracted from MD simulations [9,[14][15][16]. The derivation of the partial-slip boundary condition is given in the supplementary material, available from stacks. iop.org/njp/16/082001/mmedia.
The equations are written in non-dimensional form, where * ≡ U v v and μ * = p pR U ( ) are the non-dimensional fluid velocity and pressure, respectively. ρ is the fluid density, and U the characteristic velocity of the flow through the CNT. The Reynolds number is , where μ the fluid viscosity, and R the radius of the CNT.

Methodology
The CNT membrane is modeled as a collection of circular pipes connecting two water reservoirs. The radius = R 1.017 nm is identical to the CNT membrane used in the recent MD study of Walther et al [4]. The length (L) of the CNT varies between 3 and 7000 nm. The CNT membrane is set up so that the CNTs form a rectangular lattice with the displacements between adjacent CNTs of 18.17 nm in one direction and 17.87 nm in the other direction. Using the symmetries of the CNT membrane we reduce the computational domain to a quarter of a single CNT (see figure 1). On either side of the CNT we place a reservoir extending 50 nm away from the CNT membrane. At the inlet side of the reservoir we apply a uniform velocity boundary condition maintaining a constant volumetric flow in the range of μ − 0.8 m s 3 1 μ − 10 m s 3 1 through a single CNT. At the outlet we impose a uniform pressure boundary condition and the CNT entrance and exit have edges with a radius of curvature (r f ) ranging from 0.0 nm to 0.6 nm. Unless otherwise stated we consider sharp edges (the radius of curvature of the fillet is 0 nm). We use μ = × − − 7.2 10 kg ms 4 1 and ρ = − 997 kg m 3 for water viscosity and density, respectively, corresponding to the SPC/E water model at ambient conditions [17,18]. Here the Reynolds number is  = − Re (10 ) 3 . The incompressible Navier-Stokes equation (1) are solved using a second order finite volume method. The flow is assumed steady and we use the SIMPLE algorithm [19] implemented in the STAR-CD simulation package, into which we implemented the partial-slip boundary conditions [20]. The geometry of the computational system. The computational domain is reduced to a quarter due to the symmetry of the problem. The CNT of radius R and length L is connecting two large water reservoirs of length L r . The edges at the CNT ends are rounded with a fillet of constant curvature r f (see insert). We apply the partialslip boundary conditions to the solid walls of the membrane and CNT surfaces. The origin of the coordinate system is placed in the center of the CNT and the x axis points in the direction of the flow.

Results
We verify the accuracy of the partial-slip model by comparing the present continuum results with those obtained by related MD simulations [4]  ) obtained in this recent study. Our results on flow rates are within the error bars of the values obtained with MD simulations. We note that recently it has been reported [21,22] that in highly confined nanoscale channels an extended Navier-Stokes equation coupled to the microscopic molecular spin angular velocity describes the water flow better than the standard Navier-Stokes equation. According to our results the standard Navier-Stokes equation still gives the correct description of flow in highly confined channels when the slip length is large. Large slip lengths, corresponding to hydrophobic walls, decrease the effect of the coupling of the Navier-Stokes equation to the microscopic molecular spin angular velocity. Our slip-length = l 63 nm s is more than an order of magnitude larger than the slip length in Hansen et al [21].
The introduction of the partial-slip boundary condition has a significant effect on the pressure profile along the CNT and the water flow through the CNT. Figure 3 shows the pressure loss between the inlet and the outlet dependence on the slip length. The figure shows the results from our simulations as well as the results predicted by the model proposed by Sisan and Lichter [11] Δ π * = * s where C = 3, Δ * p is the non-dimensional pressure loss between the inlet and the outlet and * = L L R is the non-dimensional CNT length. The above equation attributes the overall pressure loss along the CNT to two contributions: pressure loss due to a fully developed flow inside the CNT given by the slip enhanced Hagen-Poiseuille equation and pressure loss at the CNT ends approximated by the pressure loss through a thin orifice [23,24]. It is interesting to observe from figure 3 and equation (3)   pressure loss does not vanish [11,25]. Figure 4 shows the pressure profile along the axis of the CNT. We observe that in the case of the standard no-slip boundary condition the pressure steadily drops through the CNT. In the case of the full slip boundary condition the pressure loss occurs solely at the CNT ends. For the partial-slip boundary condition pressure loss occurs both at the CNT ends and along the length of the CNT. The pressure loss at the CNT ends is present in the no-slip case as well. It is, however, negligible in comparison to the pressure loss inside the CNT. In the case of a large slip on the other hand the pressure loss at the CNT ends plays a significant role. We note from figure 3 that in the no-slip case the non-dimensional pressure loss The red X symbols represent pressure loss obtained by CFD simulations and the green dashed line represents pressure loss obtained from equation (3). Our error estimate due to discretization is roughly 0.2. is approximately 250 and in the full-slip case about 10. Because of large slip effects in CNT membranes we study in detail the pressure loss at the CNT ends [26,27]. The pressure loss at the CNT ends occurs due to viscosity, where a non-zero viscous part of the stress tensor can be significant at this low Reynolds number flow [28]. In the case without the partial-slip boundary condition the pressure loss at CNT ends could be estimated by dimensional analysis up to a multiplicative factor C [25], This result is in agreement to the one obtained by Weissberg [29], who studied pressure losses at tube ends analytically and obtained an upper bound for the constant, ⩽ C 3.47. A similar expression was obtained by Samson [23] and Roscoe [24] for the pressure loss through a thin orifice, effectively a tube with a vanishing length. In this case the pressure loss is known exactly and the constant C = 3.
The partial-slip boundary condition introduces a new length scale into the problem * l ( ) s that characterizes the pressure loss at CNT ends. To this end, we measure the pressure loss at CNT ends dependence on the slip length. We adopt the same expression for the pressure loss as suggested by Sisan and Lichter [11] (equation (3)) with the difference that we view C as a function of * l s and * r f . The flow upon entering the CNT fully develops into the slip enhanced Hagen-Poiseuille flow within the length of the order of the CNT radius [11]. The pressure loss due to the CNT ends is therefore confined to their vicinity. We evaluate the pressure loss due to the CNT ends by calculating the pressure loss dependence of the CNT length. We estimate C by fitting equation (3) to the obtained results. We do this for varying * l s to obtain = * C C l ( ) s and the results are shown in figure 5. The pressure loss at CNT ends clearly exhibits a dependence on * l s . Furthermore the transition of the pressure loss from the no-slip case to the full-slip case is not monotonous. Around ∼ * l 0.1 s we observe a local minimum for the pressure loss at the CNT ends and the range of C is roughly between 3.0 and 3.6. In the large slip scenario, where the majority of the pressure loss along the CNT stems from the CNT ends, an accurate description of pressure loss at the CNT ends is vital and the choice of C = 3 might not be adequate. For ≫ * l 1 s we observe C to be at its maximum. For = * l 1000 s we measure = ± C 3.583 0.005. Gravelle et al [5] in contrast give C = 3.75 for full-slip ). Gravelle et al also show that as the length of a tube approaches zero, C approaches 3. This is in agreement with analytical derivations for pressure loss through a thin orifice [23,24]. For small slip = * − l 10 s 4 on the other hand we measure = ± C 3.19 0.05, which is in agreement to Weissbergs ⩽ C 3.47. We explain the dependence of pressure loss at the CNT ends on * l s by considering the energy dissipation rate [28] ⎛ We divide the energy dissipation into two contributions:    * =˙+* * 1 2 corresponding to the energy loss outside ( * 1 ) and inside ( * 2 ) of the CNT, respectively. We assume a sharp edge at the CNT ends (i.e. = * r 0 f at the CNT ends). The sharp edges at the CNT ends imposes an effective no-slip boundary condition at the point of the corner. We assume that the velocity profile at the CNT end is independent from * l s and is equal to the velocity profile of the fluid flow through a thin orifice [29] The pressure loss outside the CNT is then equal to the pressure loss for a flow through a thin orifice and is given by equation (4) . Inside the CNT we restrict ourselves only to the axial component of the velocity. The pressure loss due to dependence of the velocity from the radial coordinate is included in the slip enhanced Hagen-Poiseuille equation (first term in equation (3)) although in the vicinity of the CNT ends not completely accurate. We thus restrict ourselves to the velocity dependence from the axial coordinate. We assume a linear transition from the velocity profile in equation (6) s s 2 and we assume that the transition happens in the length proportional to the CNT radius. Energy dissipation is then approximated by . At this * l s the velocity profile at the CNT ends (equation (6)) bears the closest resemblance to the fully developed slip enhanced Hagen-Poiseuille flow (equation (7)). This results in the least dissipated energy during the development of the flow. We also check the effect of slip at the outer walls of the CNT membrane (see supplementary material). We test three cases. One where there is no slippage present at the outer walls, one with full-slip at the outer walls and one where the boundary condition at the outer walls is equal to the boundary condition at the CNT wall. The results show that slip at the outer walls has a negligible effect on the pressure loss. This is in agreement to Gravelle et al [5] who also found the boundary condition at the outer walls to have negligible impact.
The pressure loss at the CNT ends is due to viscous dissipation of energy. In the vicinity of the CNT openings the streamlines sharply curve into and out of the CNT. This results in a nonvanishing viscous part of the stress tensor and consequentially viscous energy dissipation. This kind of viscous energy dissipation can be reduced by softening the streamline curvature. We therefore test the effect of softer corners at the CNT ends. We round the corners by introducing a fillet of constant curvature ≡ This introduces yet another length scale and we now view the constant C as a function of two variables: the slip length * l ( ) s , and the curvature radius of the fillet * r ( ) f . We measure the dependence of the constant C on the slip length as shown in figure 5. We also test our model for pressure loss at the CNT ends. The assumptions we made when developing the model are not valid any more due to the rounded corners. Therefore, we do not fit equation (10) only for C 2 but also for C 1 . We obtain  (6)) upon which we developed the model are no longer valid. We find that the rounded corners at the CNT ends reduce the energy dissipation. Figure 6 shows energy dissipation rates per unit volume for three different boundary conditions: the full-slip boundary condition, the partial-slip boundary condition with = * l 1 6 s and the no-slip boundary condition. For the full-slip case the energy dissipation rate is highest near the corner at the CNT entrance. For the no-slip case on the other hand it is highest near the walls of the CNT. For the partial slip case with = * l 1 6 s it is comparable in both regions. We also observe high pressure gradients at the wall near the transitions from the rounded corner to the outer membrane wall and the inner CNT wall cf figure 6. The high pressure gradients can be understood by considering the partial-slip boundary condition cf equation (2) and recognizing that the shear stress at the wall is linearly proportional to the flow velocity at the wall (see supplementary material). At the transition from the rounded corner to the straight wall the radius of curvature changes discontinuously from * r f to ∞ corresponding to the curvature radius of a straight wall. The discontinuous curvature results in a discontinuous wall shear stress and surface pressure gradient. The flow velocity is higher at the fillet-CNT connection than it is at the fillet-outer membrane wall connection with correspondingly higher pressure gradient at the fillet-CNT transition than at the fillet-outer membrane wall transition. In the limiting case of no-slip ( = * l 0 s ), however, the flow velocity at the wall is zero which results in a vanishing curvature term (equation (2)) and hence a continuous stress at the wall. Indeed, the pressure gradient is significantly smaller for the system with no-slip boundary condition than for the system with full or partial slip.
Finally, we perform simulations of CNT membranes with lengths varying from 3 nm to 7000 nm exceeding the current state-of-the art MD simulations [4] at a fraction of the computational cost. We measure the flow rate enhancement (E) as the ratio of the flow rate ( π * = Q ) predicted using the continuum model to the corresponding no-slip Hagen-Poiseuille value HP Combining equations (3) and (11)   independence from the CNT length [4]. The asymptotic value corresponds to the situation where the pressure loss is dominated by the friction losses inside the CNT.

Conclusion
In summary, the present study demonstrates that continuum flow models with appropriate boundary conditions can be used to accurately predict the slip enhanced flow through CNT membranes. Here we have performed finite volume simulations of the steady-state incompressible Navier-Stokes equations subject to partial-slip boundary conditions. Our results are in agreement with corresponding MD simulations and confirm the importance of pressure loss at the CNT ends. Furthermore the present simulations enable studies at several μm long CNTS and confirm the existence of an asymptotic value for the flow rate enhancement of water inside a CNT as the CNT length increases. A novel model is proposed to describe the influence of the slip length on the pressure loss at the CNT ends, which shows the existence of an optimal slip length for the pressure losses at the CNT ends. We show that the pressure loss at the CNT ends can be significantly reduced by reducing the curvature of the edges of the CNT ends. . The triangles show the flow rate enhancement for a recent MD study of water flow through CNT membranes [4]. The black triangles correspond to MD simulations using FASTTUBE and blue and green triangles to simulations using NAMD with pressure difference 200 and 20 bar respectively [30,31]. The dashed line is flow enhancement rate obtained from equation (12)