Simulation of water impregnation through vertically aligned CNT forests using a molecular dynamics method

The flow rate of water through carbon nanotube (CNT) membranes is considerably large. Hence, CNT membranes can be used in nanofluidic applications. In this work, we performed a molecular dynamics (MD) simulation of the introduction of water into CNTs in the CNT membranes, especially in vertically aligned CNT forests. The results showed that the Knudsen number (Kn) increased with an increasing volume fraction of CNT (VC) and was greater than 10−3 for each VC. Beyond this value, the flow became a slip flow. Further, the permeability increased as VC increased in the actual state calculated by the MD simulation, whereas the permeability in the no-slip state predicted by the Hagen–Poiseuille relationship decreased. Thus, a clear divergence in the permeability trend existed between the states. Finally, the flow enhancement ranged from 0.1 to 23,800, and the results show that water easily permeates as VC increases.


Simulation of water impregnation through vertically aligned CNT forests using a molecular dynamics method
Tomohiro Tajiri 1 , Ryosuke Matsuzaki 1 & Yoshinobu Shimamura 2 The flow rate of water through carbon nanotube (CNT) membranes is considerably large. Hence, CNT membranes can be used in nanofluidic applications. In this work, we performed a molecular dynamics (MD) simulation of the introduction of water into CNTs in the CNT membranes, especially in vertically aligned CNT forests. The results showed that the Knudsen number (Kn) increased with an increasing volume fraction of CNT (V C ) and was greater than 10 −3 for each V C . Beyond this value, the flow became a slip flow. Further, the permeability increased as V C increased in the actual state calculated by the MD simulation, whereas the permeability in the no-slip state predicted by the Hagen-Poiseuille relationship decreased. Thus, a clear divergence in the permeability trend existed between the states. Finally, the flow enhancement ranged from 0.1 to 23,800, and the results show that water easily permeates as V C increases.
Carbon nanotubes (CNTs) possess many excellent characteristics such as high thermal conductance, high strength, and chemical stability. Moreover, they are widely applied to electrical 1-3 , structural [4][5][6] , and biometric [7][8][9] materials. Among them, vertically aligned CNT forests (VACNFs 10 ) ( Fig. 1) have attracted great attention because they can be produced on a large scale at low cost by chemical vapor deposition 11 , and their mesoporous structures can potentially be used in nanofluidic applications such as nanofilters, biosensors, and catalysts [12][13][14] . In the design of nanofluidic equipment based on the CNTs, it is essential to understand and control the interaction between the CNTs and the fluid 12,[15][16][17] . Therefore, many analytical [18][19][20][21] and experimental [22][23][24][25][26] investigations on water permeation inside a CNT have been carried out. The flow rates of water through a CNT have been reported to be from one to five orders of magnitude greater than those predicted by the continuum-based no-slip Hagen-Poiseuille relationship. Furthermore, these flow rates increase as the area of the permeation region decreases. A similar trend was observed in the experimental results of Byeongho et al. 27 , who investigated water permeability outside a CNT in a VACNF. The similarity of the water flow inside and outside CNTs is that the flow rate and penetration coefficient increase, and the fluid flow becomes easier as the area where the fluid is impregnated becomes narrower. On the other hand, the difference is that the water flow inside a CNT is more restricted compared with that outside a CNT because the flow passes through a circular pipe. In addition, we succeed in defining the slip distance for water flow inside a CNT as well as derive the penetration coefficient formula considering the slip distance. For the water flow outside a CNT, the complex shape of the impregnation area prevents the definition of slip distance. Determining the penetration coefficient formula that considers slip is also difficult.
Although analytical and experimental flow investigations inside a CNT have been performed, no analytical investigations have been conducted outside a CNT in a VACNF. Thus, in the present study, we simulated the permeation of water outside a CNT in a VACNF and investigated the water flow. Here, we used the Gebart model 28 as the permeability model of a porous medium, which is known for its consistency with experiments using macro porous medium and continuum numerical analyses. However, the fluid inside a minute porous medium at a nanolevel scale involves a high Knudsen number (Kn), and the flow turns into a slip flow 29 , which has a high velocity in the liquid/solid interface 30 . Thus, the Gebart equation cannot be applied to a fluid in a nanopore; however, this inapplicability has not been confirmed yet. Our research goal is to verify the application of the Gebart equation by comparing the permeabilities of the Gebart equation (no-slip state) (derived using a hypothesizing Hagen-Poiseuille flow) with those in an actual state that uses molecular dynamics (MD) results and to verify the flow tendency in a nanopore.
In this study, we consider a fluid water flow with a flow front; the flow is induced by a capillary force outside a CNT. Moreover, we present our results in the following aspects: (1) investigation of the Knudsen number (Kn), (2) verification of the application of permeability predicted by a no-slip state, and (3) investigation of the divergence between the permeability in the slip and actual states by calculating the flow enhancement in a nanopore.

Results
Analytical model. The analytical model of the permeation outside the CNT is shown in Fig. 2(a). We set the CNT with a cap at the center of a graphene and arranged water molecules above the CNT. The CNT was fixed because the CNT and graphene were bonded between carbon atoms. The CNT diameter was 2.16 nm, and its length was 3.19 nm. We varied the size of the graphene to obtain V C = 0.077, 0.090, 0.106, 0.188, 0.311, 0.424, 0.505, 0.611, 0.706, and 0.780. Here, we defined V C as V C = (area of the CNT)/(area of graphene) in the model shown in Fig. 2(a). By applying periodic boundary conditions to the analytical model shown in Fig. 2(a), we created the VACNF shown in Fig. 2(b) and performed permeation outside the CNT, as shown in Fig. 2(c). Periodic boundary conditions were applied along every boundary of the unit cell. We established a sufficiently long distance with respect to the CNT axis so that the model duplicated in each direction does not affect the water where u MD is the velocity of the fluid permeation into porous media; K MD is the permeability, which is calculated by the MD simulation; μ is the viscosity of the fluid, and dP/dz is the pressure gradient in the flow direction. The Gebart equation 28 [equation (2)] is widely used to verify the macroscopic-pore permeation behavior in case the permeation direction is parallel to the fibrous direction under the assumption of a no-slip flow.
K G is the permeability calculated by the Gebart equation, and R is the fiber radius. Then, we considered K MD to be the actual permeability and K G to be the permeability of the no-slip state. We compared the permeability obtained from equation (1) using the MD simulation with that from equation (2). Furthermore, as shown in Fig. 3, we calculated dP/dz by evaluating the pressure within several sub-volumes along the CNT axis and performing a linear regression analysis. Here, we calculated the pressure in the MD using the virial theorem [equation (3)] 32 .
where V is the volume, N is the number of atoms, k B is the Boltzmann constant, T is the temperature, r ij is the distance between atoms i and j, and F ij is the interaction force between atoms i and j.
In the case of a liquid slip at the solid/liquid boundary, the actual permeability (K actual , measured from the experiment or predicted from the MD simulation) becomes larger than the calculated value from the Hagen-Poiseuille relationship under the assumption of no-slip flow (K no-slip ) 19 . This increase in K leads to the definition of flow enhancement (ε) given by equation (4) 19 . We consider K actual as K MD and K no-slip as K G in this paper.

Results of the MD simulation and calculation.
We show the results of the MD simulation and the calculation in Fig. 4. We make three analysis attempts for each V c , and the data shown in each of the figures are its average values. First, Fig. 4(a) shows the result of the mean velocity of water in the CNT axial direction at position r from the CNT center. The velocity profile is calculated in the domain from z = 0.5 nm to z = 1.0 nm at V C = 0.106. In this figure, the velocity becomes maximum near the CNT surface and minimum near the center between the CNTs. We find that a recess-shaped meniscus is formed. The same tendency is observed for other V C values. Figure 4(b) shows the result of the average permeation velocity for all the water molecules in the CNT axial direction. The figure shows that the permeation velocity does not monotonically changed as the permeation area varies, and the flow tendency across V C = 0.106 varies. In Fig. 4(c) Kn is higher than 0.01 for V C ≥ 0.188, and in Fig. 4(a), velocities exist at the CNT surface. Thus, we find that the water flow becomes a slip flow. Furthermore, Kn increases as the permeation area becomes narrower. Figure 4(d) shows the relationship between V C and permeability (log K). In this graph, although the actual permeability using the MD results increases, that of the no-slip state decreases as V C increases or the permeation area becomes narrower. On the other hand, at Vc < 0.188, the permeability of the actual state using the MD results corresponds with that of the no-slip state. Finally, Fig. 4(e) shows the relationship between V C and the flow enhancement. The range of the enhancement obtained in this research is 0.1-23800. In some areas, the permeability is more than 10 3 compared with that in the no-slip state. Furthermore, we find that the enhancement increases, and water more easily permeates.

Discussion
First, Fig. 4(a) shows that the velocity profile formed a recess-shaped meniscus. Similarly, the experimental results of Rossi 33 , who investigated the flow front of water inside a CNT using environmental scanning electron microscopy, also showed that water formed a meniscus. In the present research, we considered the water flow with a flow front, and in this flow, the surface tension exerted a strong influence. The velocity distribution was meniscus, although this does not necessarily mean that the shape is meniscus. However, because the MD model in this study assumes ordinary temperatures, the model assumes an atmosphere rather than a vacuum in areas containing no molecules. In other words, because the liquid and vapor layers co-exist, we believe that meniscus formation can occur, and the result of the mean velocity can accurately represent the permeation behavior. Moreover, Kenneth et al. 34 experimentally showed that the surface of the VACNF, which has an array of CNTs of not more than tens of micrometers, is not sufficiently hydrophobic by itself, and water droplets deposited on the surface permeate into the forest. On the other hand, in the CNT array of more than hundreds of micrometers, the top surface of the CNTs is superhydrophobic, and water does not permeate into the forest. Because we used a 3.19-nm-long CNT in the VACNF, we deduce that no contradiction occurs in the result of water permeation into the VACNF in the MD simulation.
Second, from the trend shown in Fig. 4(d), we can state that a clear divergence in the permeability tendency occurs between the slip and no-slip states, and the permeabilities show an opposite trend with V C . The trend in which the permeability increases and the fluid flows more easily as the permeation area becomes narrower corresponds to the experimental results of Byeongho et al. 27 who investigated water permeability in a VACNF. On the other hand, the permeability of the slip state corresponds to that of the no-slip state as the permeation area becomes wider. According to the above results, in the permeation in nanopore-based CNTs, a wide region exists in V C where permeability in the no-slip state cannot be applied; therefore, considering a new permeability equation with slip is necessary. Then, a wide area is necessary to prove the pressure statistics, but in the present study, we averaged over a very narrow calculation area. To clarify the effect of the variation due to the limited calculation area, we conducted an MD analysis three times under each condition. As a result, we showed that the penetration coefficient was not affected, which is the topic of this paper.
Finally, from the result of the enhancement shown in Fig. 4(e), a wide region exists in V C where the permeability is more than 10 3 compared with that in the no-slip state because the superhydrophobicity of the CNT generated almost zero friction. Moreover, weak interfacial forces exist between the CNT and water molecules, resulting in a very large enhancement. This phenomenon was believed to occur due to the superhydrophobic  surface property of CNT. However, to clarify this phenomenon even further, we believe an indicator is necessary to quantitatively evaluate hydrophobism and hydrophilism. However, this index is still currently under consideration and can be studied in the future. On the other hand, we believe that extremely valuable elements are present in the current discussion that show the applicability of the Gebart formula assuming a no-slip state derived by the Navier-Stokes equations as well as the trend in liquidity outside the CNT in the model that simulates the VACNF.

Methods
Knudsen number. For gases, Kn is defined as Kn = λ/L s , where λ is the mean free path in a gas and L s is the characteristic channel dimension. Kn is used to determine whether a flow field is a continuum 35 . In this study, the equivalent Knudsen number is calculated using lattice spacing δ according to ref. 36 instead of using the mean free path. This method of determining the Knudsen number is used because liquid molecules do not exhibit a mean free path. Here, the lattice spacing δ is defined as 37 where V 1 is the molar volume and N A is Avogadro's number. For liquid water, because V 1 is 18 cm 3 , the lattice spacing is 0.3 nm. We used equation (6) to determine Kn for the permeation outside the CNT. Here, many cases exist where the flow is approximated by reflecting the slip condition of the Navier-Stokes equations for nanoscale liquid flow 30 . The Gebart model used here utilizes the Navier-Stokes equation assuming no-slip conditions. We expect that the Gebart model will no longer be applicable once Kn > 0.01 38 because differences will appear under these conditions. Considering these points, we utilize Kn as a condition for determination.

Molecular simulation.
In the investigation, we conducted an MD simulation using the LAMMPS 39 and permeated water outside a CNT in the VACNF. In the simulation, we used the TIP3P model [40][41][42] for the water and AMBER96 43 for the potential function. The viscosity of water was μ = 0.321 mPa·s at 300 K 30 . The long-range Coulomb forces were computed using the particle-particle particle-mesh method 44 , and the mean square error was 10 −4 . The SHAKE method 45 was used to solve the equation of motion under a constrained condition. We used the NVT ensemble as the statistics ensemble and controlled the system temperature by solving the Langevin equation of motion 46 in the relaxation calculation and permeation simulation. We used the values listed in Table 1 for each parameter in the potential function and analysis conditions in the relaxation calculation. The permeation simulation conditions are listed in Table 2.