Molecular Dynamics Simulations of Nanochannel Flows at Low Reynolds Numbers

In this paper we use molecular dynamics (MD) simulations to study nano- channel flows at low Reynolds numbers and present some new interesting results. We investigated a simple fluid flowing through channels of different shapes at the nano level. The Weeks-Chandler-Anderson potentials with different interaction strength factors are adopted for the interaction forces between fluid-fluid and fluid-wall molecules. In order to keep the temperature at the required level, a Gaussian thermostat is employed in our MD simulations. Comparing velocities and other flow parameters obtained from the MD simulations with those predicted by the classical Navier-Stokes equations at same Reynolds numbers, we find that both results agree with each other qualitatively in the central area of a nanochannel. However, large deviation usually exists in areas far from the core. For certain complex nanochannel flow geometry, the MD simulations reveal the generation and development of nano-size vortices due to the large momenta of molecules in the near-wall region while the traditional Navier-Stokes equations with the non-slip boundary condition at low Reynolds numbers cannot predict the similar phenomena. It is shown that although the Navier-Stokes equations are still partially valid, they fail to give whole details for nanochannel flows.


I. INTRODUCTION
The behavior of a flow at the nanometer scale has been a subject of interest in recent years. As a typical flow in the reduced-size fluid mechanics system, the nanochannel flow embodies a series of special properties and has attracted much attention. The understanding of the physical properties and dynamical behavior of nanochannel flows has great importance on the theoretical study of fluid dynamics and many engineering applications in physics, chemistry, medicine and electronics.
It is obvious that when the system length reduces to the nano scale, the behavior of the flow is mainly affected by the movements and structure of many discrete particles. The molecular dynamics (MD) simulation, based on the statistical mechanics of nonequilibrium liquids [1], is an effective way to describe the details of a flow at the nano scale. At the same time, MD simulations also calculate physical properties of nanofluids by solving the equations of molecular motion. The role of replacing the continuum description for the nano level flows makes MD a powerful tool to study many fundamental nanofluid problems which may be extremely difficult to be implemented in the laboratory at the present time. MD simulation has been used successfully in studying aspects of molecular hydrodynamics, such as the lubrication, wetting and coating problems. Simple or complex channel flows have been investigated by several researchers, e.g. Todd et al. [2] and Jabbarzadeh et al. [3].
Couette and Poiseuille flows of simple Lennard-Jones liquids or polymers in a nanochannel of several molecular sizes demonstrate some new features and have also attracted much attention at the present. The salient features of such simple flows are the departure from the Navier-Stokes (NS) theory [4] and some different dynamic properties compared with the macroscopic systems [5]. It is well known that the quadratic velocity profile in a Poiseuille flow may be successfully predicted in the Newtonian fluid mechanics framework when the fluid density is assumed to be a temperature dependent function and does not vary appreciably over the length and time scale comparable to the molecular size and molecular relaxation time. The transport coefficients are also considered to be space and time independent in the classical continuum theory. However, as we focus attention on flows at nano level, the density becomes a function related to the position and time of fluid particles, that means, the density cannot be uniform in space and time. The irregular density distribution also affects other physical properties and hence changes the dynamical behavior of a nanofluid system. Travis and Gubbins [6] reported that the density oscillates along the flow direction with a wavelength of the order of a molecular diameter in a channel of about 4 molecular-diameter width. They also verified that the quadratic velocity profile is approximately obtained for a channel of 10 molecular-diameter width in a planar Poiseuille flow, but this phenomenon disappears when the channel width is less than 10 molecular diameters. They even indicated that the disagreement with the NS prediction occurs when the channel width is of the order of 5 molecular diameters. Koplik and Banavar [7] further pointed out that the classical fluid mechanics theory can be applied safely for Couette and Poiseuille flows in channels of width greater than 10 molecular diameters with no-slip boundary conditions. Based on such assumptions, Fan et al. [5] carried out simulations of a periodic nozzle flow of a Weeks-Chandler-Anderson (WCA) liquid.
The flow that occurs in a practical nanochannel may be more complicated than a simple Poiseuille flow in that many factors such as the channel geometry, flow type and boundary conditions may be involved altogether. The faced challenge in this area is t o find the deviation between the MD simulation and the classical NS prediction, and to obtain the qualitative difference between them with physical explanation. In this paper, we will focus mainly on flows of a simple liquid in various channels. We emplo y the MD simulations to explore the dynamic details of nanochannel flows and try to examine the limitations of the NS solutions for a liquid moving along nanochannels.

II. PHYSICAL MODELS
In this paper, we assume the liquids are constituted of quantities of spherical molecules. As a rigorous quantum-mechanical approach is not feasible at the present for a system of more than a few molecules, Newton's second law is regarded as valid to describe the molecular motion in the system. For a fluid, the momentum p i of molecule i should satisfy Newton's second law, where F i is the intermolecular force on molecule i by other molecules, and F e is the external force. In simulation, the external force may be the body force, e.g. the gravity force. The velocity v i of molecule i is related to the fluid molecule mass m By The molecular interaction forces usually depend on the physical properties and space phase structures of the fluid and wall molecules. A complex model taking into account both the fluid/fluid and fluid/wall intermolecular interactions is presented in this paper to obtain the velocity and stress distributions in a nanochannel flow. For the fluid/fluid interaction, both repulsion and dispersion effects of molecules are considered. The WCA potential, which is a modification of the Lennard-Jones potential, is adopted in the MD simulations in this paper where s is the diameter of a molecule, e is the energy parameter characterizing the molecular interaction strength and r c stands for the cutoff distance. For the fluid/wall interaction, noting the existence of the fluid -solid melting line and probable substance interchange on the boundary, a realistic sixth-power soft sphere potential is used to simulate the nanofluid system.
where ε fw is called the fluid -wall energy parameter. For the weak fluid -wall interaction, we take ε fw = ε, and for the strong fluid -wall interaction, ε fw = 3.5e. The factor 3.5 is obtained by matching the argon potential well on a smooth wall (Heinbuch & Fischer [8]). In MD simulations, the solid walls may be represented by several layers of wall molecules which are located at the sites of a planar face-centered lattice or a body-centred cubic lattice. Each wall molecule is assumed to be anchored at its lattice site by a stiff spring. A Hookean spring force is introduced here as the external force for wall molecules. For simplicity, we adopt the two-layer linear spring model. The potential may be written as where δ the displacement of the wall molecule from its lattice site, and C is the stiffness of the spring.
To keep the molecular oscillation small from its site and to overcome the slip problem, a relatively soft wall spring with a stiffness of C = 75e/ σ 2 is adopted. From the viewpoint of a microscopic system of particles, the stress tensor is just a result of the space phase distribution and the dynamic properties generated by the movements of molecules in the system. According to the Irving-Kirkwood method [9], the contribution of each particle to the stress tensor may be divided into two parts, a kinetic component related to the velocity distribution and a configuration component related to the position distribution of particles. The stress tensor is written as The angular brackets denote the ensemble average, n is the density. The first sum on the right-hand side indicates the contribution from the momentum transfer where m is the molecular mass, u a and u ß are the peculiar velocity components of molecule i in the a and ß directions. The second sum represents the potential contribution where r a and F ß are the a and ß components of the distance vector and the potential force vector between two molecules. For a simple channel flow, the viscosity of the fluid is calculated by the formula µ = τ/γ, the shear rate γ is usually acquired by the finite difference method.

III. SIMULATION METHOD
We have outlined in the previous section the molecular models to be used in our MD simulations. The simulation details of a liquid flowing through various nanochannels will be given here. For brevity, we use a reduced system to perform the simulations. T he units of the characteristic quantities are listed in Table 1 for constructing the dimensionless parameters used in the simulations.
As the work done by the external force may partially be converted into heat in the system, the wall/fluid interaction also contributes to the increase of temperatures of the fluid and the wall. It is important to keep the nano-scale fluid system at a fixed temperature during the simulation. There are usually two ways to keep the temperature of the system constant, using a thermostat or adjusting the temperature by re-scaling the molecular peculiar velocity during the simulation. The latter method requires more calculation work, as the approach is to adjust the wall temperature to a given constant at each time step [8]. It produces satisfactory results in many cases as long as the wall and fluid temperatures link appropriately. As regards to the use of thermostats, some researchers have argued that although it imposes additional constrains on the equations of motion, it improves computation efficiency when applied to a more complex system and sometimes even enhances the accuracy of the simulation [3]. In order to avoid some unphysical effects in the simulation, a Gaussian thermostat method [1], which has been approved to be a good thermostat method, is adopted here to obtain reasonable simulation results. Actually, the Gaussian Method is to add an external force with the following term (7) where Pi is the peculiar momentum of molecule i and coefficient ξ represents the thermostatting multiplier for keeping the kinetic temperature fixed in the system, where γ i is the local shear rate of flow field at r i and P ix and P iy are the x and y components of the peculiar momentum of molecule i.
In this simulation, we have developed a FORTRAN code program called MDFluid. MDFluid consists of three modules: data preprocessor, MD simulation and data output. The first part is used to generate the initial configuration of the computation domain. The second part is the core module and its main function is to realize the MD simulation for a nanofluid and to acquire the flow field via some initial and boundary conditions. It is open to adopt different algorithms for the computation of velocity and temperature. The last part is the data output interface. Currently MDFluid is only suitable for several simple nano flows, it is under further development by us. For simplicity, the mass of wall molecules is usually assumed to be identical to that of fluid molecules. The initial configuration of wall molecules are generated separately by the preprocessing module and read in as input data. The total number of molecules depends on the size and geometry of the computational domain and the densities of the fluid and wall materials. All fluid molecules are initially located at the sites of a face-centered lattice. The initial velocities of inner fluid and wall molecules are set randomly according to the given temperature and the inlet and outlet velocities are assigned the same velocity profile due to the periodic boundary condition. At the beginning of simulation, fluid molecules are allowed to move without applying the external force until a thermodynamic equilibrium state is reached. Then the external force field is switched on and the non-equilibrium simulation starts. During the course of simulation, the equations of motion are solved by the leap-frog method. Periodic boundary conditions are applied on the fluid boundary of the computational domain.
To verify and validate the NS equations in a nano scale system is obviously an important and serious task at the present. In this paper, we try to check the validity of the NS equations in a nanochannel flow. Considering a body force e xerting on the system, the governing equations in dimensionless form may be written as where the Reynolds number is defined as Re = ρUL/µ, the fluid density ρ and viscosity µ are assumed to be constants and U and L are t he characteristic velocity and length of the domain field. The dimensionless body force is F body = ReLg/U 2 e x , where g is the gravity constant and e x is the unit vector in the x direction. We use the finite element method (FEM) to obtain the flow field of a nanochannel flow. The same initial and periodic boundary conditions are applied in the FEM simulations. To compare the MD and FEM results, we use the Reynolds number as an adjustable parameter in the FEM to match the flux of the MD simulation. In the f ramework of continuum fluid mechanics, the Reynolds number plays an important role in analyzing the flow field. Any two solutions with the same Reynolds number and initial/boundary conditions must be similar. In nano flows, the characteristic length reduces to the nano size and the Reynolds number becomes extremely small. In this sense, all nano flows may be classified as low Reynolds number flows and the knowledge on low Reynolds number flows can be inherited in some ways. We wish to find similarities and differences between the NS equations and the MD simulations when they are applied to the same nano fluid system in this paper.

IV. RESULTS AND DISCUSSION
Let us examine five complex channels with different geometry. The first one is a simple plane channel whereas the other four have concave or convex surfaces. The computational domain for the planar flow of a WCA liquid is shown in Figure 1. The system is surrounded by periodic images of itself in x dimension. The computational domain is 0 < x < 200 and -10 < y < 10, the fluid density and viscosity may be determined during simulation.
The periodic dynamic conditions are also used in the course of simulation. The spring constant C of wall molecules is C = 75 and the gravity constant g = 0.5. We assume the fluid/fluid interaction energy parameter e = 1 and the wall/fluid interaction energy parameter is characterized by ε wf = ε or ε fw = 3.5ε according to the interaction intensity. For convenience, we compute the five nanochannel flows with the Reynolds number based on the width of the channel being 3 (first three cases) and 5 (last two cases).  The x -direction velocity profile at position x = 100 is plotted in Figure 3. We find from Figure 3 that the velocities obtained by the MD simulation agree very well with the FEM results. This verification in a simple Poiseuille flow encourages us to apply MDFluid to more complex cases.  This observation can further be confirmed by the x -direction velocity profile in Figure 5, from which we notice that there are obvious negative velocity molecules near the walls (y = 7.5 and y = -7.5). The reverse momenta induce the generation and development of nano-size vortices under the non-slip boundary conditions. Case 3 is a similar channel with convex surfaces. We also set 1600 fluid molecules including 250 wall molecules and strong fluid/wall interaction in the MD simulation. In this case, both predictions on the flow field are analogous (see Figure 6). We cannot see any vortices in both results. One possible reason is that the stream velocities maybe somewhat small compared with those in Case 2 (see Figures  5 and 7). In case 2, the concave surfaces make the middle channel like a contraction nozzle and increase the stream velocity. In this case, on the other hand, the corresponding momentum cannot form vortices in the cavity. We next consider the channel with semi-circular surfaces. In order to watch the flow more clearly, we increase the system to 2100 fluid molecules and 360 wall molecules in the course of MD simulation. The Reynolds number is adjusted to 5. The nano vortices reappear in the MD simulation, while the FEM still predicts no any abrupt velocity changes after the middle nozzle area. In Figure 9, we also find the negative velocities near the walls. Case 5 is a simple alteration of Case 4, the concave surfaces are replaced by the semi-circular convex surfaces, and all other conditions are kept the same. The vortices occur in the dents due to the accurate simulation by large number of molecules in the MD simulation, while the FEM presents no striking flow attribution (see Figure 10). The velocity profile shows negative velocities near the wall region. From the above analysis, we find that: (1) The NS equations and the MD simulations are both applicable for nanochannel flows. The velocity profiles still maintain quadratic for nano Poiseuille flows in the MD simulations (see Figure 12). (2) The MD simulation can determine the nano vortices for flows at low Reynolds numbers. The nano vortices may be generated near the walls in complex flows due to the large molecular momenta. In the meantime, the continuum NS equations cannot predict nano flows well near the walls, as they fail to give flow details in the whole field. To this point, we can conclude that large deviations between the NS equations and the MD simulations exist in areas far from the core of the flow.
(3) A guarantee for a good MD simulation is to include as many molecules as possible into the system and to deal appropriately with wall molecules and the boundary conditions.

V. SUMMARY
In this paper we adopt MD simulations to investigate various nanochannel flows at low Reynolds numbers and to present some new interesting results. The WCA potentials with different interaction strength factors are employed for the interaction forces between fluid -fluid and fluid -wall molecules. In order to keep the temperature at the required level, a Gaussian thermostat is introduced in our MD simulations. Comparing the flow fields obtained from the MD simulations with the predictions by the classical NS equations at the same Reynolds number, we find that both results agree with each other qualitatively in the central area of a nanochannel. However, large deviations usually occur in areas far from the core. The nano-size vortices due to the large momenta of molecules in the near-wall region may be found in some complex nanochannel flows. But the traditional NS equations with the no-slip boundary condition at low Reynolds numbers cannot predict the similar phenomenon. It is shown that although the NS equations are still partially valid for nano flows, they fail to give whole details of the flow fields.