1 Introduction

The recent decade has witnessed the huge increase in the number of studies for the simulation of molecular flows in medicine, physiology, nano scale devices, etc.^[1-4]. According to the researches of Travis et al.^[5] and Bird^[6], Navier-Stokes equations are not applicable for domains on the molecular scale. Published reviews of nano scale flow problems indicate that a molecular dynamics (MD) approach has been widely used to simulate the nano scale flow problems^[7-8].

The simulation of MD systems was first introduced by Alder and Wainwright^[9]. In their work, atoms were assumed to be hard spheres. Later on, Rahman^[10] used Lennard-Jones potential to simulate soft spheres. Since then, MD simulations using Lennard-Jones potential have become prevalent among other researchers. Progressive usage of MD led to simulations of complex molecular systems such as liquid water^[10-11], liquid n-butane^[12], and proteins^[13]. These are just examples of a few early works of numerous molecular simulations in a wide range of applications to analyze nano scale fluid flow systems. Although the full particle MD simulation would provide a more realistic simulation for the nature of the flow^[14], this approach carries high computational costs with it. The most important affecting parameter on the computational cost is the number of simulated particles. The number of particles depends on different flow parameters such as the flow density and the domain length. Even with the rapid development of computer processors, which has contributed greatly to the improvement of the efficiency of large scale particle system simulations, the time scale of the simulations is limited to the order of nanoseconds^[15]. This improvement is inadequate for many problems of nano particle systems. Therefore, researchers have always been looking for ways to reduce the computational cost.

Since about 90% of the CPU time in the most of MD simulations^[15] is used for the computation of inter-molecular forces, Taylor et al.^[16], Plimpton^[17], and Murty and Okunbor^[18] developed parallel algorithms using the force decomposition technique to speed up the MD simulations. The growing application of biomolecule flow simulations has brought researchers' attention to the limitations of computational cost in MD simulations. To solve this problem, they have developed the coarse graining technique to provide more time and length scale information. In this approach, a limited number of atoms will be considered as one particle. Therefore, the number of particles in the simulation will be decreased significantly^[19].

Authors of this article believe that one of the effective techniques to reduce the computational order of MD simulations is the symmetric solution. In the numerical solution of continuum flows such as computational fluid dynamics (CFD), for domains with the symmetric nature, a symmetric solution is applicable by implementing the appropriate boundary condition on the symmetry plane, which decreases the simulation run time by up to 50% with regard to the run time of the whole domain solution^[20]. Review of the literature shows that, other than the previous works^[21], the symmetric solution has not been studied in nano scale fluid flows. To investigate possibility of solving the half symmetry of the whole domain, Karimian and Amani^[21] proposed a new boundary condition for the symmetry plane. This boundary condition was first used to solve the Poiseuille flow. Very good results have been obtained with the symmetric solution in comparison with the whole domain solution for the Poiseuille flow, which encourages us to apply this new symmetry boundary condition to other molecular channel flows.

For this purpose, two general geometries containing converging-diverging and diverging-converging nano scale ducts are selected. The periodic boundary condition is applied in the flow direction and in the direction normal to flow, and the wall boundary condition is applied in the other direction. To the best of our knowledge, no experimental data exist in the literature studying the velocity profile either in a converging-diverging or diverging-converging nano scale duct. Therefore, for the validation purpose, velocity profiles of the symmetric solution with the proposed boundary condition on the symmetry plane and those of the whole domain solution are compared with each other. Properties of a molecular flow can be extracted by sampling within a control region in the domain, called the bin^[2]. As a matter of fact, extracted velocity profiles are highly influenced by the bin size and the employed averaging method^[22]. The appropriate bin size is determined via simulations which will be presented in Subsection 2.3. The averaging method based on implementations of sample-averaged measurement (SAM) and cumulative average measurement (CAM) called SAM-Modified-CAM (SMC)^[23] is used in this paper. In Section 3, principals of the proposed symmetric solution will be explained. In Section 4, symmetric solutions of molecular flows in two geometries are validated with their whole domain solutions. Then, to demonstrate that the validation of the proposed symmetry boundary condition is held for a wide range of cases in these two types of geometries, simulations with different values of effective parameters including the driving external force of the flow, the density of the flow, and the domain length are conducted.

2 Simulation 2.1 Solution domain and initial setup

Due to the successful implementation of the proposed symmetry boundary condition to the constant area channel, two other channels with the varying area are selected to examine the applicability of the proposed boundary condition for more general cases. Two selected domains including converging-diverging and diverging-converging nano scale ducts are shown in Fig. 1. Boundary conditions applied in x-and y-directions are periodic. However, the wall atoms are located in the z-direction. The height and the width of the both geometries are L_z=L_y=3.615 nm. The length of the domain is not fixed since applicability of the proposed symmetry boundary condition is studied for different channel lengths. Atoms with argon's properties are used to generate the fluid within the channel. In addition to the channel length, the effect of density has also been studied in this paper. The flow density is defined as ρ =n/V^*, where n is the number of atoms, and V^* is the non-dimensional volume, defined as V^*=V/σ³ (σ is the characteristic length of argon atoms). Argon fluid atoms are initially setup with a face-centered cubic (FCC) lattice form^[23]. Initial thermal velocities of atoms are calculated according to the designated initial temperature of T_ini=300 K with directions determined by the Maxwell-Boltzmann distribution function^[24-25]. Two rows of atom layers shape the walls in both geometries. Physical properties of wall atoms including the atomic mass and the atomic diameter are the same as those of argon fluid atoms. For the first row of wall atoms, atoms are centered every 0.4016 nm in the x-direction and 0.45187 nm in the y-direction. The second row of wall atoms is a copy of the first row which is shifted 0.2008 nm in the x-direction, 0.22595 nm in the y-direction, and 0.1795 nm in the z-direction. According to this arrangement, the oblique distance between two atoms in the z-direction would be 0.35158 nm. Although this value is slightly larger than the argon atomic characteristic length σ (0.34 nm), our experience suggests that fluid atoms cannot escape through the wall.

To constrain the wall atoms in their positions and let them only oscillate around their initial positions, a spring force with the stiffness of 500ε/σ² (ε is the energy scale of argon atoms equal to 1.67× 10^-21 J) is applied to each of the wall atoms^[21]. The wall temperature is kept constant at 300K by rescaling the velocity of the wall molecules at each time step^[26]. To ensure continuous flowing of fluid, an external force is applied to all of the fluid atoms^[27]. Note that all the units in this paper are in the international system of units.

2.2 MD simulation

The motion of atoms in the time is simulated using Newton's second law. The acting inter-molecular force is determined by , in which U_ij is the Lennard-Jones potential function^[27], which is described as

(1)

where r_ij is the distance between two interacting atoms. Based on the position of atoms in the solution domain, each iteration potential function between every two atoms is determined. For each atom, all of the inter-molecular forces applied from the neighboring atoms are calculated. Summation of these forces is used with Newton's second law to determine the new position of the atom. Since the neighboring atoms far from the atom under consideration in the solution domain have very few effects on the atom motion, only forces exerted by molecules within a cut-off distance of 2.5σ to the atom will be considered^[23]. The Verlet algorithm is used to integrate the equation of motion in the time. In this paper, a time step of dt_MD=0.0023 τ is used, where τ is the characteristic time of Lennard-Jones potential defined as τ≡(mσ²/ε)^0.5[28]. m and σ are equal to 39.948 g/mol and 0.34 nm, respectively. Based on these values, the time step in the present study is set to 0.01 ps. In all of the simulations presented here, 2000000 time steps will be taken, while only the last 1500000 time steps are used for the sampling, as the system has well reached the equilibrium by then. It is worth mentioning that MD simulations in the present study are carried out using an in-house-developed MD code parallelized by MPI library.

2.3 Domain binning

Flow properties such as the mean velocity should be obtained from the solution with the averaging method. The averaging method used in this paper is the SMC, which is thoroughly described in Ref.[23]. The bin size is also an important parameter in the averaging process. Large bins would miss the local accuracy, and small bins would not include enough sampling data which cause undesirable fluctuations in the extracted profiles^[22]. Therefore, the first step is to decide on the appropriate bin size in which reliable profiles with the local accuracy can be obtained. Since in a constant area duct, the velocity profile changes only in the z-direction (between two walls), the solution domain is merely divided into 10 cubic bins in the z-direction. In other words, the bin height (h_bin) would be set to the 10% of domain height, and the bin width and length would be equal to the width and length of the domain. For a varying-area duct, the velocity profile will change along the duct, and the bin length should not be equal to the length of domain. In this case, the velocity profile in the middle of the duct will be used for the comparison purpose. Therefore, the length of bin is gradually decreased from its initial length that is equal to the length of duct, in order to reach an appropriate size. For a length of one σ, discrepancies have been observed in the velocity profile due to lack of sufficient data. After some tests, a length of 1.5σ or 0.51 nm is appropriate for extracting the velocity profile in the middle of the solution domain, as shown in Fig. 2. Note that 10 bins used in the z-direction for the whole domain simulation mean that there are 5 bins in the z-direction for symmetric simulations. The center of the fifth bin is half of a bin (h_bin/2) lower than the symmetric plane. Thus in Figs. 5-10, there is a distance gap between the maximum bin position of symmetric simulation profiles and the symmetric plane which is equal to (h_bin/2).

3 Proposed symmetry boundary condition

Domains with the symmetric nature could be divided by a plane into two equal parts. In the case that just one part is solved, the computational cost would be halved. The most important point, however, is to define an appropriate boundary condition along the symmetry plane. Our experience has shown that the conventional symmetry boundary condition applied in the numerical solution of continuum flows such as CFD is not applicable in MD simulations, at least for the Poiseuille flow case^[21]. Therefore, it is necessary to search for an appropriate symmetry boundary condition for micro scale flows.

As stated before, the solution domain is divided into two vertical parts, conveniently named as the lower and upper parts. Assume that the lower part is going to be solved. Therefore, we have to replace the effect of atoms of the upper part on the lower part with the symmetry boundary condition. Obviously, atoms of upper part may approach the symmetry plane. These atoms can affect atoms of the lower part which are within their cut-off distance. As a result, atoms of the lower part that appear within a distance equal to or less than the cut-off distance (2.5σ) from the symmetry plane will be affected by the atoms of the upper part. This region of lower part is called the "boundary zone". This boundary zone is shown in Fig. 3. The symmetry boundary condition should properly represent the interaction between atoms of the upper part and atoms of the boundary zone. The effect of this interaction will eventually spread out in the lower part through the inter-molecular actions between atoms of the lower part. In Ref.[21], a number of assumptions for the symmetry boundary condition were numerically tested for Poiseuille flows by the authors of the present paper. Among them, one boundary condition which worked very well was proposed for the symmetry boundary condition.

If the proposed symmetry boundary condition is applied, then the only part that will be solved is the lower part and the atoms within it. This means that we will not have any information from the atoms of the upper part. However, the effect of the upper atoms on the boundary zone in the lower part should be considered with the applied boundary condition. To follow the general concept of symmetry, in each iteration, a "reflected zone" which is a symmetric copy of the described boundary zone of the lower part is generated above the symmetry plane, as illustrated in Fig. 3. Note that all the atoms in the reflected zone are just a symmetric copy of atoms in the boundary zone. The position of these virtual atoms, in each iteration, is updated based on the position of their relative boundary zone atoms in that iteration. Therefore, we do not need to calculate the acting inter-molecular forces on them or include them in the time integration process. We just take into account the effect of these atoms in the virtual reflected zone on the atoms in the boundary zone. Now, besides the effect of atoms in the reflected zone on atoms in the boundary zone, simulations of MD in the lower part are carried out.

One may assume that the approaching atoms reaching the symmetry plane will meet with their symmetric copy atoms of the reflected zone. Therefore, they will bounce back and will not pass the symmetry plane. It seems that a symmetry boundary condition based on this approach results in a very good prediction of flow variables including the velocity profile in comparison with those of the whole domain solution. Quite interestingly, our experience reported in Ref.[21] did not prove that. Based on the extensive numerical testing of the authors, it is found that very good results will be achieved if the atoms in the reflected zone are shifted along the flow direction by half of the domain length, as shown in Fig. 4. It should be noticed that the periodic boundary condition is applied in the flow direction. Therefore, by shifting the reflected zone, its atoms which pass the periodic boundary in the flow direction, would be set in the first half. Therefore, in the proposed boundary condition, in each iteration, a symmetrical copy of the boundary zone, termed as the reflected zone, is shifted into the flow direction by a length which is equal to half of the domain length. The positions of these atoms are not independent as they are just a copy of atoms in the boundary zone. These atoms will simulate the influences of the atoms in the upper part on atoms which are in the boundary zone and in the lower part, respectively.

Before the implementation of the proposed symmetry boundary condition in the micro nozzle flow, this symmetry boundary condition is applied to a simple channel flow (Poiseuille flow) of argon with the density of 0.641 in a channel with L_x=14.46 nm and L_y=L_z=3.615 nm. The external driven force is equal to 1.8 pN. The results could be seen in Fig. 5 that, other than a small difference close to the symmetry plane, very good agreement is reached between the velocity profile of the symmetric solution and that of the whole domain solution. A previous comprehensive study by the authors on this case presented in Ref.[21] confirmed that a symmetric solution with the proposed boundary condition was valid in the Poiseuille flow.

Having tested successfully the proposed boundary condition in a Poiseuille flow, now we can implement it to solve flows in variable-area ducts.

4 Results and discussion

In this section, results obtained with the MD simulation of flows in two types of ducts are presented to investigate the applicability of the symmetry boundary condition proposed in Ref.[21] to molecular flows in channels with varying areas. In each case, the results obtained with the symmetric solution are compared with the results obtained with the whole domain solution in the same problem. Therefore, in each case, both the symmetric solution and the whole domain solution are calculated. From the results, the mean velocity profile in the channel is selected for the purpose of comparison. Molecular flows of argon in geometries shown in Fig. 1 are solved for the density of 0.55 and domain lengths of 9.81 nm for the first case and 6.4 nm for the second case. The driving external force is equal to 0.72 pN for both cases. Velocity profiles obtained with the whole domain solution and the symmetric solution are compared with each other in Fig. 6. In both cases, especially for the second case, very good agreement has been observed. These results encourage us to investigate the generality of the proposed idea of symmetric solution. For this purpose, a series of test cases are defined to examine if the proposed method is valid under variation of different parameters including the flow driving external force, the fluid density, and the domain length.

4.1 Effect of flow driving external forces

For a fixed fluid density of 0.48, external forces of 0.09 pN, 0.18 pN, 0.36 pN, and 0.72 pN are applied to the fluid atoms to drive the flow. In all of the simulations of this section, the domain lengths are equal to 9.81 nm for the first case and 6.4 nm for the second case. Mean flow velocity profiles obtained in the middle of channel with both the symmetric solution and the whole domain solution are compared with each other in Figs. 7 and 8 for the first and second cases shown in Fig. 1, respectively.

All of the results shown in Fig. 7 show very good agreement between velocity profiles of the symmetric solution and the whole domain solution. Although it seems that under lower external forces, the symmetric solution slightly over-predicts the velocity profile, while under higher external forces, it slightly under-predicts the velocity profile.

Results of Fig. 8 which are for the second case of Fig. 1 however show better agreement in comparison with those shown in Fig. 7. Excellent agreement is seen for results with higher external forces, although small discrepancies are close to the symmetry plane.

Although a wider range of external forces can be examined, for the present range of external forces which change the maximum velocity in the channel by a factor of 10, the proposed symmetry boundary condition has worked very well.

4.2 Effect of fluid density

The fluid density is another important parameter in MD simulation whose effect on the proposed symmetric boundary condition is considered in this paper. For a fixed external force of 0.36 pN, molecular flows in both geometries are solved for fluid densities of 0.48, 0.55, and 0.65. In all of the simulations, the domain length is equal to 9.81 nm for the first case and 6.4 nm for the second case. Mean flow velocity profiles of both the symmetric solution and the whole domain solution are compared with each other in Fig. 9.

As shown in Fig. 9, all of the results show excellent agreement between the velocity profiles obtained with the symmetric solution and the whole domain solution for a variety of fluid densities. It seems that the density variation does not have a significant effect on the results obtained from the proposed method.

4.3 Effect of domain length

Since in the proposed method, the reflected zone is shifted by a length equal to the half length of domain, one question that arises is how significant the length of domain would be in this method. All of the above results obtained for the first geometry (converging-diverging nozzle) are with the length of 9.81nm. For the fixed external force of 0.36 pN and the density of 0.55, molecular flows are solved in this case for three channel lengths of 4.81 nm, 9.81 nm, and 19.86 nm. For the second geometry, the same flow defined above is solved within domains with lengths of 3.39 nm, 6.4 nm, and 13.28 nm. Mean flow velocity profiles obtained in these cases with both the symmetric solution and the whole domain solution are compared with each other in Fig. 10.

According to Fig. 10, for both the first and second cases, symmetric solutions for relatively short domains of 4.81 nm and 3.39 nm both under-predict the velocity profile relative to the velocity profile obtained with the whole domain solutions. By increasing the domain length to 9.81 nm in the first case and 6.4 nm in the second case, very good results of velocity profile of the symmetric solution are obtained. However, both relatively longer domains of 19.86 nm in the first case and 13.28 nm in the second case slightly over-predict the velocity profile with respect to that of the whole domain solution. However, the overall agreement among all of the results is very good. Note that an exact agreement in MD simulations where chaotic nature of collision between particles is dominant is not generally expected.

From the discussion above, one can conclude that the proposed boundary condition for the symmetric solution in the constant and variable-area channels such as those presented in Fig. 1 is an appropriate and fairly accurate approach. In all of the simulations, the CPU time of the symmetric solution is about 50%~55% of the related whole domain solution.

5 Conclusions

A new boundary condition is proposed to be applied for the solution of symmetric molecular flows in channels. This boundary condition is implemented to solve molecular flows inside a channel (the Poiseuille flow) and two different variable-area channels. The results obtained with the symmetric solution are in very good agreement with the results obtained with the whole domain solution. Generality of the proposed boundary condition is examined under various combinations of important flow parameters which govern the flow, namely, the external force driving the flow, the density of the flow, and the channel length. In each case, the results are compared with those of the whole domain solution. In all of the symmetric solutions, the CPU time is reduced by about 50%~55% of the CPU time of the corresponding whole domain solution. From the results of these cases, it is concluded that the proposed boundary condition and the predicted flow field are valid for micro scale fluid flows in domains with symmetric nature, at least for domains described in this paper. The symmetric solution saves the time and computational cost substantially.