Introduction

With the fast development of micro/nanoscale technologies, an increasing number of two-dimensional (2D) nanomaterials have been developed in laboratories and used in practical applications. [1][2][3][4] As a typical 2D nanomaterial, black phosphorus (BP) has attracted enormous attention because of its direct bandgap, 5 high carrier mobility, 6 and anisotropic properties. 7 The outstanding physical, 8 electronic, 9 and optical properties 10 of BP have driven the development of micro/ nanofluidic devices 11 based on the material. However, current research on BP-based micro/nanofluidic devices is still at the exploratory stage, focusing on fluid flow characteristics in nanochannels 12,13 and solid-liquid interface mechanisms. 14,15 The interaction mechanisms at the solid-liquid interface in BP nanochannels remain elusive. Therefore, unveiling the interplay between the BP-based nanochannel and fluids is crucial to understanding the transport properties of fluids in such nanochannels for practical applications.

In recent years, molecular dynamics (MD) simulations have emerged as a widely-used tool to investigate fluid transport at the molecular scale and have been extensively applied to study the pressure-driven flow of liquid hydrocarbons in nanochannels. 7,[16][17][18][19][20][21][22][23][24][25][26][27][28] The liquid molecules confined in a nanochannel behave differently compared to the bulk state because of the more dominant role of interaction between the molecules and confined nanochannel walls. [29][30][31] The flow behavior of molecules in BP nanochannels becomes more complicated because of the solid-liquid interaction with the anisotropic lattice structure of BP. 32 In a recent study, water flowing in BP nanochannels exhibited a discrepancy (maximum 2 times) of flow velocity along the zigzag/armchair directions. 7 However, whether the simulation results in water as liquid phase can be generalized to other liquids (e.g., hydrocarbons) is not clear. Particularly, the flow behavior of long-chain hydrocarbons in nanochannels is highly related to their molecular orientation. 28 Therefore, investigating the flow behavior of hydrocarbon molecules confined in BP nanochannels with anisotropy is warranted.

In this paper, MD simulations are conducted to investigate the anisotropic flow behavior of liquid oil confined in BP nanochannels. As a representative oil model n-dodecane (C 12 ) is employed in the simulations. The velocity profiles of oil molecules are analyzed by changing the external driving force and nanochannel width. The average velocities of C 12 molecules flowing along the armchair and zigzag directions differ by a factor of B9 under the same conditions due to the anisotropic lattice nature of BP, which is significantly higher compared to the reported water flow in BP channels with a velocity difference of 2 times along the armchair and zigzag directions. Furthermore, we investigate the friction factor at the interface between the nanochannel and C 12 molecules and the orientation of the C 12 molecules near the BP walls. Results show that oil flow along the armchair and zigzag exhibit a large velocity discrepancy, which can be influenced by the external driving force and width of the nanochannels. These findings provide a better understanding of the transport properties of the C 12 molecules in BP nanochannels under different conditions.

Methods and simulation model

In this work, all simulations are performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 33 and visualized by Open Visualization Tool (OVITO). 34 Fig. 1a shows a schematic of the simulation setup. Two BP layers with a dimension of B5 nm (x) Â 5.5 nm (z) are used as the nanochannel walls with 180 C 12 molecules sandwiched between the two parallel BP walls while the atoms number of different conditions is listed in Table S1 (ESI †). Periodic boundary conditions are applied along the x and z directions. The time step size is set to 1 fs. The cutoff distance is set to 12 Å. The C 12 molecules are modeled using the polymer consistent force field (PCFF). 35 The BP molecules are modeled using the Stillinger-Weber (SW) potential. 7 The non-bond interactions are described by the Lennard-Jones (L-J) potential, 36,37 which has included the impact of van der Waals interactions: 36

where e and s are the energy and length constants, respectively, and r ij is the distance between two atoms i and j. For cross-species pairwise L-J interactions, the Lorentz-Berthelot rule is used: 38,39 e ij ¼ ffiffiffiffiffiffiffiffi e ii e jj p

;

where e ij and s ij are the energy constants and distance parameters of the L-J potential between type i and type j atoms, respectively. All adopted L-J potential parameters 7 are listed in Table 1. The longrange electrostatic interactions in the entire system are computed by the particle-particle particle-mesh (PPPM) approach 40 with an accuracy of 1 Â 10 À5 .

The system is relaxed by the following processes: first, the system is run in the microcanonical ensemble (i.e., NVE) for 10 ps. Then, a canonical ensemble (i.e., NVT) is applied at a temperature of 300 K for 2 ns. Subsequently, the system is equilibrated by applying an isobaric-isothermal ensemble (i.e., NPT) for 10 ns at 1 atm with one of the BP nanochannel walls fixed by a spring force of 2 Â 10 À2 kcal mol À1 Å À1 while the other is free to move. Finally, a following NVT for 5 ns is applied. The density profiles of C 12 molecules along the y direction are collected and averaged during the last 5 ns of the NVT process. As shown in Fig. 1b, the average density of C 12 in the middle region of the BP nanochannels is calculated to be 0.708 g cm À3 , which is within 94.53% to the experimental density data of its bulk phase (0.749 g cm À3 ) 41 at the same pressure and temperature. The density of C 12 is calculated with the same relaxed process without BP wall as shown in Fig. S3 (ESI †).

The flow behavior of C 12 molecules confined in the BP nanochannel is studied by nonequilibrium molecular dynamics (NEMD). During the simulations, both BP nanochannel walls are fixed, and constant external driving forces along the z direction (i.e., armchair direction) or x direction ´-Hoover thermostat applied along the directions perpendicular to the C 12 molecules flow 28,[42][43][44] since employing the velocity perpendicular to the flow direction and subtracting the center of mass velocity help accurately control the temperature. The interactions between BP nanochannel and C12 molecules under different conditions are quantitatively analyzed using the group force computation method (Table S1, ESI †).

Results and discussion

For the simulated model (Fig. 1a), the zigzag direction of the BP is in the x direction, and the armchair direction is along the z direction. External driving forces ranging from 5 Â 10 À4 kcal mol À1 Å À1 to 2 Â 10 À3 kcal mol À1 Å À1 are applied to each atom in the C 12 molecules along either armchair or zigzag direction of the BP nanochannel. The group force computation method applied to calculate the interaction between BP nanochannel and C 12 molecules. This interaction influences the pressure gradience (Table S1, ESI †) and velocity profile of the C 12 molecules under various conditions. The velocity profiles of C 12 molecules in the relaxed system along the y direction are averaged every 20,000 time steps (i.e., 2 Â 10 À2 ns) based on the data from the last 10 ns of NVT. Fig. 2 shows the velocity profiles of C 12 molecules confined in the BP nanochannels when the C 12 molecules are flowing along the armchair (Fig. 2a) and zigzag (Fig. 2b) directions under different driving forces. Under the same external driving force, the average velocity of C 12 molecules flowing along the armchair direction is much lower than that flowing along the zigzag direction. When C 12 molecules are flowing along the zigzag direction of a 2-nm nanochannel under an external driving force of 5 Â 10 À4 kcal mol À1 Å À1 , the average velocity of C 12 molecules is 9 times higher than that flowing along the armchair direction. Moreover, the velocity profiles of C 12 molecules flowing along the armchair direction exhibit larger curvatures than those along the zigzag direction, which can be attributed to the larger friction 32 between BP walls and C 12 molecules along the armchair direction. This can also explain why the average velocity along the armchair direction is smaller than that along the zigzag direction. In addition, with the increase of the external driving force, the average velocity change along the armchair direction is higher than that along the zigzag direction (Fig. S1, ESI †). It can be inferred that the friction between the C 12 molecules and the BP walls changes as the external driving force increases.

To further explain the different flow behaviors of the C 12 molecules along the armchair and zigzag directions, the Darcy-Weisbach friction factor is calculated based on the data collected from the MD simulations by the following formula: [45][46][47]

where r is the mass density, f ext is the external driving force added on every fluid atom while N is the number of the fluid atoms, A is the sectional area of the nanochannel and L is the length of nanochannel segment, v is the average fluid velocity at a nanochannel cross-section, and D h is the hydraulic diameter which is twice the width of the nanochannel. Fig. 3 depicts the calculated friction factors at the interfaces between C 12 molecules and BP nanochannels with a width of 2 nm under different external driving forces. The magnitude friction factor varies from 1 to 10 4 , which agrees well with the simulation results of other nanochannels in prior work. 46,47 Results show that the friction factor along the zigzag direction is consistently smaller than that along the armchair direction. On the other hand, the friction factor decreases as the external driving force increases along both directions. Consequently, the disparity between friction factor along the armchair and zigzag directions reduces as the external driving force increases (Fig. S2, ESI †). The friction factor results are consistent with the velocity profile results shown in Fig. 2.

To further elucidate the mechanisms of the drastically different flow behaviors along the armchair and zigzag directions, the orientation of the C 12 molecules near the BP wall is investigated because it can affect the molecules flowing at the interface. 28 The orientation parameter (S) of each C 12 molecule is calculated by: 48,49

where y(x) is the angle between the flowing direction and the headto-tail vector of the C 12 molecule when flowing along the zigzag direction, and y(z) is that of the C 12 molecules flowing along the armchair direction. The value of S ranges from À0.5 to 1.

Particularly, S = 1 indicates that the C 12 molecules are aligned parallel to the flow direction, while S = À0.5 indicates that the C 12 molecules are aligned perpendicular to the flow direction.

Fig. 4a exhibits the x-z cross-sectional view of the onelayered oil molecular structure close to the BP wall. The chain-like C 12 molecules are parallel to the BP wall and aligned in different directions. As Fig. 4b shows, under an external driving force of 1 Â 10 À3 kcal mol À1 Å À1 , the orientation parameter S of C 12 molecules flowing along the zigzag direction is higher than that along the armchair direction under the same external force. The result implies that along the zigzag direction, most C 12 molecules are oriented parallel to the flow direction. On the contrary, most C 12 molecules are oriented perpendicular to the flow direction along the armchair direction.

As such, when the C 12 molecules flow along the zigzag direction, the friction between the C 12 molecules near the BP walls is low, leading to a high average velocity. In contrast, when the C 12 molecules flow along the armchair direction, the molecules in the neighboring layers need to overcome the corrugated landscape, resulting in a lower average velocity. Therefore, we can conclude that for C 12 molecules flowing along armchair and zigzag directions, the alignment of oil molecules near the wall is different, leading to substantial divergences in the friction factor and average velocities.

Additionally, Fig. 5 and Fig. S4 (ESI †) illustrate that as the width of the nanochannel increases, the velocity profiles along both the armchair and zigzag directions retain their highcurvature and low-curvature configurations, respectively, 32 indicating that the friction between the BP walls and C 12 molecules flowing along the armchair direction is consistently higher than that along the zigzag direction. 32 The mechanism of the different velocity profiles will be analyzed by following the friction factor and orientation data with different width of nanochannels.

The comparison between Fig. 5a and b with Fig. 5c and d and further results (Fig. S4, ESI †) demonstrate that the curvature of the velocity profiles declines as the width of the nanochannel increases. This suggests that the velocity difference between C 12 molecules in the middle region and the nearwall region of the nanochannel is gradually decreasing because the relative importance of the friction from the BP walls is Fig. 4 (a) x-z cross-sectional view of the one-layered structure, in which the molecular chains of C 12 are aligned parallel to the BP wall. (b) The orientation parameter S value versus time along armchair and zigzag directions in 2-nm BP nanochannels.

progressively diminishing. As the width of the nanochannel increases, the interaction between the BP walls and the C 12 molecules in the middle region of nanochannels gradually weakens. On the other hand, when the nanochannel widens, the velocity gradient difference along the y direction decreases due to the increase in the total number of atoms, making the velocity profiles appear smoother.

In Fig. 6a and b, we present the logarithm of friction factor at the interface between BP nanochannel surface and C 12 molecules versus external forces for BP nanochannels with widths of B4 nm and B6 nm. Similar to the results for the 2-nm-wide nanochannel, the friction factors along the armchair direction are significantly higher than those along the zigzag direction. As the external driving force increases, the friction factors along both directions decrease. Similarly, the friction factors along the zigzag direction experience a smaller reduction than along the armchair direction with the increment of the external driving force. This reduces the disparity between the friction factors along the two directions in each nanochannel, and these findings are consistent with the velocity distribution results. Moreover, the friction factors along different directions decrease as the width of the nanochannel increases under identical driving forces. For instance, under the same external driving force of 5 Â 10 À4 kcal mol À1 Å À1 , the friction factor of the 4-nm-wide nanochannel (Fig. 6a) is 2 times higher than that of the 6-nm-wide nanochannel (Fig. 6b). The discrepancy between the friction factors along different directions diminishes when the width of the nanochannel increase.

To further reveal the influence of nanochannel width on the flow behavior, we calculate the orientation parameters of the C 12 molecules close to the BP walls in 4-nm and 6-nm nanochannels under identical external driving forces. Fig. 7a and b shows that the orientation parameters S calculated from eqn (5) along the zigzag direction are larger than those along the armchair direction in the 4-nm and 6-nm nanochannels. This implies that along the zigzag direction, most C 12 molecules are oriented parallel to the flow direction, whereas along the armchair direction, most C 12 molecules are oriented perpendicular to the flow direction. This can also explain why the frictional resistance along the zigzag direction is lower than that along the armchair direction under the same conditions. By comparing the orientation parameters of C 12 molecular in 4-nm (Fig. 7a) to 6-nm (Fig. 7b) nanochannels, the difference gradually diminishes between the orientation parameters S along the armchair and zigzag directions. This indicates that the disparity in molecular orientations of the C 12 one-layered structure near the BP wall along the direction of motion is decreasing along different directions, which is consistent with the results of the friction factor and velocity distribution differences.

Conclusion

systematic molecular dynamics simulations have been performed to study the flow behavior of C 12 molecules in BP nanochannels. Average velocities of C 12 molecules along armchair and zigzag directions of BP walls under different conditions were calculated and compared. Results show that the C 12 molecules exhibit prominent anisotropic flow behavior in the BP nanochannels. When C 12 molecules are flowing along the zigzag direction of a 2-nm nanochannel under an external driving force of 5 Â 10 À4 kcal mol À1 Å À1 , the average velocity of C 12 molecules is 9 times higher than that flowing along the armchair direction. The difference in average velocities of oil flowing along armchair and zigzag directions depends on the external driving force applied to the C 12 molecules and the flow resistance, while the difference in friction is related to the orientation of C 12 molecules in the near-wall region. Along the zigzag flow direction, most C 12 molecules in the layered structure are aligned parallel to the flow direction, resulting in smaller friction and higher velocities. In contrast, along the armchair flow direction, most C 12 molecules in the one-layered structure are aligned perpendicularly to the flow direction, leading to larger friction and therefore lower velocities. Moreover, the molecular orientation also explains the velocity profiles of the C 12 molecules in the BP nanochannels under varying width of nanochannels. These findings facilitate a better understanding of the flow behavior of fluid molecules under nanoconfinement in BP nanochannels and may shed light on the solid-liquid interactions in nanochannels.