Nanofluidics and NanofluidsView this Special Issue
Molecular Dynamics Simulation of Nanoscale Channel Flows with Rough Wall Using the Virtual-Wall Model
Molecular dynamics simulation is adopted in the present study to investigate the nanoscale gas flow characteristics in rough channels. The virtual-wall model for the rough wall is proposed and validated. The computational efficiency can be improved greatly by using this model, especially for the low-density gas flow in nanoscale channels. The effect of roughness element geometry on flow behaviors is then studied in detail. The fluid velocity decreases with the increase of roughness element height, while it increases with the increases of element width and spacing.
Micro/nano-electromechanical systems (MEMS/NEMS) have received considerable attentions over the past two decades. Fluid flows are usually encountered in these systems [1–3]. Fluid transport and interaction with these systems serve an important function in system operations . Understanding the behaviors and manipulations of fluids within nanoscale confinements is significant for a large number of applications [5–7].
The effect of the wall serves as a distinct feature of fluid flow in micro/nanoscale-confined devices [8–10]. The wall plays an increasing role in fluid flow when decreasing the flow characteristic length scale. Barisik and Beskok found that, in a channel with 5 nm in height, 40% of the channel is immersed in the wall force field . Therefore, the fluid transport characteristics, such as momentum and energy, significantly deviate from predictions of kinetic theory . Therefore, the effect of this near-wall force field on the nanoscale channel flow must be understood and evaluated.
Molecular dynamics simulation (MD) investigates the interactions and movements of atoms and molecules, using N-body simulation . This method has been employed by many researchers in the past to study the liquid flow in nanochannels [13–16]. Recently, the MD simulation is also adopted to investigate the gaseous flow in nanoscale-confined channels [11, 17–19]. Barisik and Beskok [11, 17] investigated shear-driven gas flows in nanoscale channels to reveal the gas-wall interaction effects for flows in the transition and free molecular regimes. Hui and Chao  studied the gas flows in nanochannels with the Janus interface and found that the temperature has a significant influence on the particle number near the hydrophilic surface. Recently, Babac and Reese  investigated classical thermosize effects by applying a temperature gradient within the different-sized domains.
In some MD simulations, idealized-wall models are considered. The interactions of fluid-wall atoms are usually considered as functions, for example, the diffuse and specular reflections, Maxwell's scattering kernel , or Cercignani–Lampis model . These idealized-wall models are feasible in some specific situations. However, when we study the detailed flow behaviors in the rear-wall region, the atomic-wall model must be considered. But the atomic-wall model is expensive both in computational time and memory. In confined channel flows, most atoms are requisite to describe the atomic wall. The number of wall atoms is much larger than that of fluid molecules. This drawback is particularly fatal for the gas flow. For example, Barisik et al.  studied a nanoscale Couette flow at Kn = 10. The simulation box is 162 nm × 3.24 nm × 162 nm. In their study, the number of gas molecule is 4900, while the number of wall atom is 903003. As a result, most of the computational resource is consumed on the computation of wall atoms.
Recently, Qian et al.  proposed a virtual-wall model for the MD simulation to reduce the computing time. The unit cell structures are infinite repetitive in the atomic wall. As a result, the force on a fluid molecule from wall molecules is periodical. This force was first calculated and stored in memory. During the simulation, when a fluid molecule moves into the near-wall region, the force on this fluid molecule from wall molecules can be determined directly, according to the position of the molecule relative to the wall. The near-wall region here refers to the region near the wall with distance smaller than the cutoff radius. Excessive calculations of fluid-wall interactions can be avoided, and the computing time can be reduced drastically. The time reduction is more significant for lower fluid density in nanoscale channels.
In present study, the virtual-wall model is adopted to describe the rough wall. The remainder of this paper is organized as follows. Section 2 introduces the MD simulation and the virtual-wall model. Section 3 describes the application of this model to the rough wall. Finally, Section 4 elaborates the conclusions of the study.
2. MD Simulation and Virtual-Wall Model
In the present MD simulation, interactions between fluid-fluid atoms and fluid-wall atoms are both described using the truncated and shifted Lennard–Jones (LJ) 12-6 potential given as follows:where rij is the intermolecular distance between atoms i and j, ε is the potential well depth, σ is the atomic diameter, and rc is the cutoff radius. Lorentz–Berthelot mixing rule  is employed to calculate the LJ parameters between fluid-wall atoms.
In the virtual-wall model, the force on a fluid atom from wall atoms can be expressed aswhere N is the number of wall atoms which interact with the fluid atom. The atomic wall is composed of FCC lattices and the unit cell structures in repetition. When wall atoms are fixed to their lattice point, the force on the fluid atom is periodic in both x and z directions. For example, the force of a fluid molecule located at x, y, and z is exactly the same as the force of the same molecule located at x + iL, y, and z + kL, where i and k are integers and L is the lattice constant. If the force distributions in the unit cuboid domain (L × rc × L) are known, then the force can be determined anywhere else. This is the core concept of the virtual-wall model.
The virtual-wall model for the smooth wall is first examined. Without losing generality, gas argon flow confined between FCC platinum walls is considered. The walls are along the xz plane and the simulation box are periodic in both x and z directions. For argon-argon interactions, σAr and εAr are 0.3405 nm and 119.8kB, respectively. For argon-platinum interactions, σAr-Pt is 0.3085 nm and εAr-Pt is 64.8kB, according to the Lorentz–Berthelot mixing rule . In this study, rc is set as 0.851 nm, which is approximately equal to 2.5σAr. The masses of argon and platinum atoms are 6.64 × 10−26 kg and 3.24 × 10−25 kg, respectively. These parameters have been validated in previous studies [25, 26].
The simulation box is set to be 40.9 nm × 17.1 nm × 40.9 nm in x, y, and z directions. A force of 0.008εAr/σAr is acted on each gas molecule as an external force  to drive the gas to flow in the nanoscale channel. The atomic-wall model is also carried out here to make a comparison. The thickness of the wall is 1.18 nm, which is larger than the cutoff radius. The lattice constant of the FCC platinum lattice is 0.393 nm.
In the MD simulation, the neighbor-list method is used to calculate the force between atoms while the velocity-Verlet algorithm is adopted to integrate the equations of motion . The timestep in the simulation is set to be 10.8 fs. The first 1 million steps are used to equilibrate the system, and another 5 million steps are used to accumulate properties in the y direction, with the bin size to be 0.0614 nm. The Langevin thermostat method  is employed to control the gas temperature before equilibrium. Only thermal velocities are used to compute the temperature and pressure. The above parameters and techniques are adopted in all simulations.
The open-source MD code called large-scale atomic/molecular massively parallel simulator (LAMMPS) , developed by Sandia National Laboratories, is adopted to carry out the MD simulations.
The density and velocity profiles across the nanoscale channel calculated using the atomic- and virtual-wall models are compared in Figure 1. Perfect agreement between these two models can be found, which indicates that the virtual-wall model works well in the MD simulation.
In order to compare the computational time, these two simulations are performed on a single Inter i7-4790K CPU processor. The computational time for the virtual-wall model is 0.4 h, while for the atomic-wall model, the time is 67.5 h. The virtual-wall model is much more efficient in the present case.
3. Rough Wall Simulations
3.1. Virtual-Wall Model for the Rough Wall
From the micropoint of view, all walls are rough. Surface roughness plays an important role in fluid flow and heat transfer . So, in the present study, the virtual-wall model is adopted to describe the rough wall.
In the present study, platinum atom cuboids on the smooth atomic wall are used to represent the roughness element, as illustrated in Figure 2. The roughness element is periodic in both x and z directions. The geometry of the roughness element is shown in Figure 2(b). The height of the roughness element is h, and the widths in x and z directions are both l. The spaces between two elements in x and z directions are both L.
In order to perform the virtual-wall model, a unit cuboid is first introduced, as shown in Figure 2. The rough wall can be considered as the close-packed array of this unit cuboid. The size of the unit cuboid is L × H × L, where H = h + rc. Fluid molecules interact with wall atoms only when they are located within these cuboids. When fluid molecules are outside these cuboids, the distances are larger than rc and no interactions between fluid and wall atoms are needed.
The cuboid is periodic in both x and z directions. Therefore, the force of a fluid molecule located at x, y, and z is exactly the same as the force of the same molecule located at x + iL, y, and z + kL, where i and k are integers. If the force distribution in the unit cuboid domain (L × H × L) is known, then the force on a molecule anywhere else can be deduced. The unit cuboid domain is then divided into MX × MY × MZ bins, and the forces in each bin are calculated and stored in the memory . During the simulation, the corresponding force of a fluid molecule located in the near-wall region is called directly from memory according to its position.
The virtual-wall model for the rough wall is first validated. Argon molecules are supposed to flow between nanoscale rough platinum walls. The simulation setup is the same as in Section 2. For the roughness element, h = l = 2a and L = 4a, where a is the lattice constant of the FCC platinum lattice, which is 0.393 nm. In the simulation, gas density is set to be 7.17 kg/m3. The Knudsen number, which is defined as the ratio of gas mean free path to the channel height, is 0.95, and the flow is in transition regime. In order to make a comparison, the atomic-wall model is also carried out here. In the simulation, 3087 gas argon atoms and 218406 wall platinum atoms are used.
The density and velocity profiles of the virtual-wall model are shown in Figure 3. These profiles are compared with the corresponding atomic wall simulation. Perfect agreement between these two models can be found, which indicates that the virtual-wall model works well for the gas flows in rough wall channels.
The gaseous flows in nanoscale channels with smooth and rough walls are first compared. The schematic diagram of channel geometry is shown in Figure 4(a). Three channels are investigated. The outer channel and the inner channel are both smooth, with channel heights equal to H′ and H′ − 2h, respectively. Here, h is the height of the roughness element. The third channel is rough, with the channel height equal to H′, and the roughness element height is equal to h. In the simulation, H′ is 15.35 nm and h is 0.786 nm. So, the height of the inner channel is 13.67 nm. The other parameters are kept the same as in Section 2.
The velocity profiles for these three channels are shown in Figure 4(b). It can be found that the velocity of the rough channel is much smaller than those of smooth channels. It is well known that, in nanoscale channel flows, the wall plays an extremely important role in the fluid flow. Here, in the rough channel, the total surface area is much larger than those in smooth channels because of the existence of roughness elements. As a result, the collision probability between fluid-wall atoms is larger and more fluid molecules are affected by the wall in the rough channel. So, the fluid velocity of gas in the rough channel is smaller. The effect of roughness is of great importance to nanoscale channel flows.
3.2. Roughness Element with Different Heights
The influences of roughness element geometry on flow behaviors are then studied. Roughness elements with different heights are first studied. The widths l and the spacing L of the roughness element are kept the same, while the element height h is variable. Three element heights (h = a, 2a, and 3a) are considered.
The velocity profiles of the rough wall with different element heights are shown in Figure 5. It can be found from the figure that the fluid velocity decreases with the increase of element height. This is because the total surface area is larger at higher element height. According to the explanation in Section 3.1, the wall effect is larger at higher element height. So, the fluid velocity is smaller.
Fitting curves are obtained for each velocity profile at different roughness element heights, based on the gas velocity in the central part of the channel. From the fitting curves, we can deduce the slip velocity on the wall conveniently. It can be found from Figure 5 that the slip velocity also decreases with the increase of element height.
3.3. Roughness Element with Different Widths
Roughness elements with different widths are then studied. The height h and spacing L of the roughness element are kept the same, while the width l is variable. Three roughness element widths (l = a, 2a, and 3a) are considered.
The velocity profiles at different roughness element widths are shown in Figure 6. It can be found from the figure that the element width has a great influence on the velocity profile. The fluid velocity increases with the increase of element width. The total surface areas are the same in these three cases, so are the wall effects, according to Section 3.1. However, at large roughness width, for example, l = 3a, the gap between two roughness elements is small. As a result, it is hard for the gas molecules to enter into the gap, because of the repulsive force between fluid-wall atoms, according to (1). That is to say, the effective surface area diminishes. So, the fluid velocity increases in the rough channel with the increase of the element width.
The fitting curves obtained for each velocity profile at different roughness element widths are also shown in Figure 6. It can be found that the slip velocity increases with the increase of the element width.
3.4. Roughness Element with Different Spacings
Roughness elements with different spacings are studied at last. The height h and width l of the roughness element are kept the same, while the spacing L is variable. Three roughness element spacings (L = 4a, 6a, and 8a) are considered. Other parameters are kept the same as introduced above.
The velocity profiles of rough walls with different element spacings are shown in Figure 7. It can be found from the figure that the fluid velocity increases with the increase of element spacing. This is because the total surface area is smaller at larger element spacing. According to the explanation in Section 3.1, the wall effect is smaller at larger element spacing, so the fluid velocity is larger.
The corresponding fitting curves for each velocity profile at different roughness element spacings are also shown in Figure 7. The results show that the greater the spacing, the larger the velocity slip.
The wall plays an extremely important role in the nanoscale channel flows. In the present study, MD simulation is carried out to investigate the nanoscale gas flows in rough channels. The virtual-wall model for the rough wall is proposed, and its validity is confirmed. The computational efficiency can be improved greatly by using this model, especially for the low-density gas flow in nanoscale channels. The effects of roughness element geometry on flow behaviors are then studied in detail.
From the simulations, we found that the total surface area is of great importance in nanoscale channel flows. The fluid velocity is inversely proportional to the total surface area. The fluid velocity and velocity slip decrease with the increase of roughness element height, while they increase with the increase of element width and spacing.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This work was supported by the National Key R&D Program of China (Grant no. 2017YFB0603701) and National Natural Science Foundation of China (Grant nos. 11672284 and 11602266).
M. W. Collins and C. S. König, Micro and Nano Flow Systems for Bioanalysis, Springer, New York, NY, USA, 2012.
M. Gad-el-Hak, MEMS: Introduction and Fundamentals, CRC Press, Boca Raton, FL, USA, 2010.
G. Karniadakis, A. Beskok, and N. R. Aluru, Microflows and Nanoflows: Fundamentals and Simulation, Springer, New York, NY, USA, 2006.
S. Colin, Microfluidics, ISTE, Ltd. and John Wiley & Sons Inc., London, UK, 2013.
D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, London, UK, 2002.
G. Karniadakis, A. Beskok, and N. R. Aluru, Microflows and Nanoflows: Fundamentals and Simulation, Springer, New York, NY, USA, 2006.
M. Barisik, B. Kim, and A. Beskok, “Smart wall model for molecular dynamics simulations of nanoscale gas flows,” Computer Physics Communications, vol. 7, pp. 977–993, 2010.View at: Google Scholar
D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, New York, NY, USA, 2004.