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.

Carbon nanotubes (CNTs) possess many excellent characteristics such as high thermal conductance, high strength, and chemical stability. Moreover, they are widely applied to electrical123, structural456, and biometric789 materials. Among them, vertically aligned CNT forests (VACNFs10) (Fig. 1) have attracted great attention because they can be produced on a large scale at low cost by chemical vapor deposition11, and their mesoporous structures can potentially be used in nanofluidic applications such as nanofilters, biosensors, and catalysts121314. In the design of nanofluidic equipment based on the CNTs, it is essential to understand and control the interaction between the CNTs and the fluid12151617. Therefore, many analytical18192021 and experimental2223242526 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.

Figure 1

VACNFs.

(a) General and (b) top views of a VACNF.

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 model28 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 flow29, which has a high velocity in the liquid/solid interface30. 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 = 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 as V = (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 molecules in the cell. Then, we investigated the impregnation of water through capillarity and did not apply any force to the impregnation.

Figure 2

Analytical models and permeating behavior.

(a) Analytical model of impregnation outside a CNT. (b) Impregnation is performed by applying the periodic boundary condition. (c) Behavior of permeating water outside CNTs.

Fluid permeation outside CNTs

In the evaluation formula of fluid permeation of the actual state calculated by the MD simulation, we used the following permeability in the Darcy’s law31 [equation (1)]:

where u is the velocity of the fluid permeation into porous media; K 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 equation28 [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 is the permeability calculated by the Gebart equation, and R is the fiber radius. Then, we considered K to be the actual permeability and K 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.

Figure 3

Pressure gradient.

Axial pressure gradient between the CNTs at intervals of 0.2 nm. dP/dz is calculated by evaluating the pressure within several sub-volumes along the CNT axis and performing a linear regression analysis.

where V is the volume, N is the number of atoms, k is the Boltzmann constant, T is the temperature, r is the distance between atoms i and j, and F is the interaction force between atoms i and j.

In the case of a liquid slip at the solid/liquid boundary, the actual permeability (Kactual, 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 (Kno-slip)19. This increase in K leads to the definition of flow enhancement (ε) given by equation (4)19.

We consider K as K and K as K 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, 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 = 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 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 = 0.106 varies. In Fig. 4(c) Kn is higher than 0.01 for V 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

Results of the MD simulation and calculation.

(a) Velocity (V = 0.106) of the CNT axial direction in position r from the CNT center, and its velocity profile is calculated in the domain from z = 0.5 nm to z = 1.0 nm. The vertical dashed line at y = h represents the CNT surface. (b) V versus velocity of water along the impregnation direction. Its velocity is calculated in the computational domain from z = 0.5 nm to z = 1.0 nm. V versus (c) Kn, (d) permeability, and (e) enhancement.

Figure 4(d) shows the relationship between V 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 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 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 103 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 Rossi33, 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. 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 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 where the permeability is more than 103 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, where λ is the mean free path in a gas and L is the characteristic channel dimension. Kn is used to determine whether a flow field is a continuum35. 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 as37

where is the molar volume and N is Avogadro’s number. For liquid water, because is 18 cm3, the lattice spacing is 0.3 nm. We used equation (6) to determine Kn for the permeation outside the CNT.

L corresponds to the hydraulic diameter28 that considers the control volume [Fig. 5(b)] of the model where the CNTs stand close together and straight on the substrate [Fig. 5(a)]. Here, many cases exist where the flow is approximated by reflecting the slip condition of the Navier–Stokes equations for nanoscale liquid flow30. 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.

Figure 5

Definition of the control volume to calculate the hydraulic diameter.

(a) Top view of the interface in a two-dimensional patterning surface with hydrophobic posts organized on a square lattice. (b) Definition of the control volume.

Molecular simulation

In the investigation, we conducted an MD simulation using the LAMMPS39 and permeated water outside a CNT in the VACNF. In the simulation, we used the TIP3P model404142 for the water and AMBER96 43 for the potential function. The viscosity of water was μ = 0.321 mPa·s at 300 K30. The long-range Coulomb forces were computed using the particle–particle particle-mesh method44, and the mean square error was 10−4. The SHAKE method45 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 motion46 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.

Table 1

Parameters of AMBER96 potential function.

Atoms(ij)	εij [kcal/mol]	σij [Å]
CC	0.086	3.407
OO	0.1521	3.151
HH	0.0	0.0
CO	0.114	3.275
CH	0.0	0.0
Bonds(ij)	KijR [kcal/mol/Å2]	Rijeq [Å2]
OH	1000.0	1.0
Angles(ijk)	Kijkθ[kcal/mol/rad2]	θijkeq [°]
HOH	1,000.0	109.47

i, j, and k represent an atom.

Table 2

Computational condition of simulation.

Relaxation	Time step [fs]	0.01
	Running time [fs]	20
	Damping constant	100.0
Simulation	Time step [fs]	1.0
	Running time [fs]	8.0 × 105
	Temperature [K]	300

Additional Information

How to cite this article: Tajiri, T. et al. Simulation of water impregnation through vertically aligned CNT forests using a molecular dynamics method. Sci. Rep. 6, 32262; doi: 10.1038/srep32262 (2016).