Event-Driven Molecular Dynamics Simulation of Hard-Sphere Gas Flows in Microchannels Volkan

Classical solution ofNavier-Stokes equationswith nonslip boundary condition leads to inaccurate predictions of flow characteristics of rarefied gases confined in micro/nanochannels. Therefore, molecular interaction based simulations are often used to properly express velocity and temperature slips at highKnudsen numbers (Kn) seen at dilute gases or narrow channels. In this study, an eventdrivenmolecular dynamics (EDMD) simulation is proposed to estimate properties of hard-sphere gas flows. Consideringmolecules as hard-spheres, trajectories of the molecules, collision partners, corresponding interaction times, and postcollision velocities are computed deterministically using discrete interaction potentials. On the other hand, boundary interactions are handled stochastically. Added to that, in order to create a pressure gradient along the channel, an implicit treatment for flow boundaries is adapted for EDMD simulations. Shear-Driven (Couette) and Pressure-Driven flows for various channel configurations are simulated to demonstrate the validity of suggested treatment. Results agree well with DSMC method and solution of linearized Boltzmann equation. At low Kn, EDMD produces similar velocity profiles with Navier-Stokes (N-S) equations and slip boundary conditions, but as Kn increases, N-S slip models overestimate slip velocities.


Introduction
Over the last decades, micro-and nanoelectromechanical systems (MEMS, NEMS) have been rapidly developed.Today, many applications like flow sensors or miniaturized jet engines are taking place in industry.Therefore, investigation of gas flows in micro-and nanochannels is essential for design and operation of microdevices such as micropumps, microvalves, and microturbines [1].As devices get smaller, Knudsen number, Kn (the ratio of the mean free path of gas molecules, , to the characteristic length of the channel,   ), increases.Magnitude of Kn determines the flow regime: continuum regime for Kn ≤ 0.01, slip regime for 0.01 < Kn ≤ 0.1, transitional regime for 0.1 < Kn ≤ 10, and free molecular flow for Kn ≥ 10.In many macroscale applications,  is negligible compared to the channel size (Kn ≪ 0.01) and the continuum assumption in Navier-Stokes equations is valid.For smaller scales, channel size is in the order of ; flow is either in slip or in transitional regime, and compressibility and rarefaction effects are present.Since Navier-Stokes equations are hardly valid in these regimes, fluid can be treated as an ensemble of particles interacting with each other and the boundaries.Molecular Dynamics (MD) and Direct Simulation Monte Carlo (DSMC) methods are two simulation methods based on the computation of motions and interactions of particles involved in such ensembles.In molecular simulations, computational difficulties arise as the number of molecules gets larger simply because interacting particles and their positions need to be computed for each interaction.
DSMC method [2] uses a number of representative molecules to simulate larger number of real molecules.Motions of molecules are exact but collisions are generated in a probabilistic manner.On the other hand, MD simulations are much more realistic and accurate since each particle represents a real molecule and its position and velocity are exactly known.Classically, trajectories of molecules are calculated from integration of Newton's equations of motion over a time step regarding interaction potential.Hence, this process is time driven.Although the potential is smooth and continuous for real interactions, approximate approaches are used especially in dilute gas simulations for the sake of computational efficiency.Lennard-Jones potential is the most common simplified interaction potential model used in MD simulations.Standard MD simulation based on Lennard-Jones interaction potential was first introduced by Rahman [3] and Verlet [4].This model has been used in numerous studies up to present time.The major drawbacks of this approach are limitations of simulation time and size [5].Integration time step is so small that, during a simulation of ten microseconds, the enormous number of time steps is required even for ten thousand molecules representing a very small volume.
Modelling the system as an ensemble of hard particles is an alternative to continuous interaction potential approach.Interaction potential between hard particles is discrete; that is, no force is exerted on the molecules except impacts.Assuming no external force field, resulting trajectories are linear.Intermolecular collisions are instantaneous and binary.When a collision occurs, postcollisional velocities are analytically determined via conservation of energy and momentum.Analytically predicted collision times yield the decomposition of simulation into a series of asynchronous events.This process is event driven since simulation time advances directly from one event to the next.Due to its nature, event-driven approach eases simulation of larger systems for longer periods compared to time-driven one.Since the first introduction of event-driven molecular dynamics (EDMD) simulations [6] developments of more effective algorithms have further enhanced performance of EDMD [7][8][9][10][11][12].Simulation of millions of particles for longer periods is possible with current computational power of a desktop computer [13].
In order to simulate a confined gas flow using molecular simulations, system boundaries and interactions with the gas molecules must be modelled appropriately.Three types of boundary are usually sufficient to represent a flow in a microor nanochannel.These well-known boundaries are periodic, wall, and flow boundaries.
Since the performance of MD simulations is highly affected by the number of simulated molecules, it is desired to represent the computational domain with the smallest possible sample.Periodicity is the most common boundary model for such representations.It also removes surface effects in the simulations of infinite systems.In the presence of periodic boundaries, computational domain is repeated in the direction of periodicity to form an infinite lattice.When a molecule reaches a periodic boundary it continues its course on the opposite face with exactly the same velocity vector.
On the wall boundary, a molecule hits the wall if the distance between the wall and molecular centre is a half molecular diameter.The molecule is reflected back with a new velocity determined according to the wall model.Specular reflection is the most basic wall model in which surface is perfectly smooth at molecular level.When a molecule undergoes specular reflection, normal velocity component is inverted while tangential components remain unchanged, conserving tangential momentum.However, in reality, wall surface is rough and the molecules are reflected at random angles from the wall.Diffuse reflection is the most common model to represent such surfaces.In diffuse model, postreflection velocity of the molecule is substantially independent of the incoming velocity and determined stochastically from a distribution based on wall temperature and velocity.
Several treatments for flow boundaries have been developed to induce a flow inside the channel.These are selective reflection treatment [14][15][16][17][18], acceleration body force treatment [19,20], and the implicit treatment for flow boundaries [21,22].In selective reflection and acceleration body force treatments, total energy and the number of molecules in the system are conserved since no new molecule or energy is introduced.In the implicit treatment for flow boundaries, molecules reaching either end (upstream or downstream) of the system leave the computational domain permanently.Moreover, new molecules are introduced to the computational domain according to local domain properties.Hence, unlike aforementioned methods, a number of molecules in the computational domain and total energy are not preserved.Although these treatments are convenient to use both in time and event-driven simulations, the implicit treatment in DSMC method has not been applied to MD simulations until this study.
Investigation of Shear-and Pressure-Driven flows of a hard-sphere gas in microchannels is beyond the capabilities of MD simulation methods using continuous interaction potentials.Additionally, implicit treatment of flow boundaries is not easily applicable since the method is time driven.Therefore, simulations were carried out with the use of eventdriven molecular dynamics simulations in this study.EDMD has similarities in implementation with both classical MD and DSMC methods.In contrast to DSMC method, collisional dynamics such as movements of molecules, determination of collision partners, and calculation of postcollisional velocities are deterministic in the hard-sphere EDMD.All collisions are real and computed in MD simulations.In this study, simulations were conducted with real number of molecules taking advantage of an event scheduling algorithm with (1) computational complexity.Added to this, a cell division method increased computational speed greatly without missing any event.On the other hand, as a major contribution of this study, wall and flow boundaries were handled stochastically as in DSMC method which created an opportunity to adapt implicit treatment for flow boundaries to EDMD simulations for the first time.
This study is organized as follows: EDMD algorithm and application of boundary treatments are described in Section 2. In Sections 3 and 4, theory and simulation results of Shear-(Couette) and Pressure-Driven flows are presented.Concluding remarks are given in Section 5.

Methodology and Model
2.1.Computational Model.In this study, EDMD simulation of monatomic molecules was conducted.Computational domain was a rectangular parallelepiped region.Molecules were modelled as hard-spheres; thus only translational energies were considered.Only binary collisions, as a valid assumption for dilute gases, were employed.Postcollisional velocities were calculated according to mass, diameter, and precollisional velocity information of colliding pair by satisfying conservation equations of energy and momentum.
Argon (Ar) was selected as the monatomic gas of this study.In order to reduce numerical errors, all magnitudes were nondimensionalized by considering unit Boltzmann constant (  = 1.0) and scaling with the reference variables temperature ( ref = 300 K), mass ( = 6.62 × 10 −26 kg), and diameter ( = 3.41×10 −10 m).It is worth mentioning that the results of the study can be extended to any monatomic gas by selecting proper scaling variables.
Due to the nature of MD simulations, immense computational resources are required.Simulation speed decreases almost linearly as the number of molecules in the system increases.In order to speed up the simulations, computational domain was partitioned into cubical cells using the method described in detail by Kandemir [9].A cell has 26 neighbours at maximum, fewer if it is at an edge of computational domain.In this method, molecules are assigned to the cells according to their current positions.An intermolecular collision is probable only if the distance between two molecules is less than cell width.Otherwise molecules would encounter cell crossing events beforehand.Thus, only the molecules of same or neighbour cells are considered as candidate collisional pairs at a given time.Implementation of cell partitioning method reduces the computational effort to (  ) for one molecule where   is the average number of molecules in the neighbourhood.Since there are  molecules in the system, total computational effort is about () assuming that (  ) ≈  (1).Simulations of the same system on the basis of single and multicell approaches yield exactly identical numerical results using the same seed and random number generator.Hence, one can conclude that the multicell method can guarantee that no collision is missed.
Three types of events are possible in event-driven simulations.These are intermolecular collisions, moleculeboundary interactions, and cell crossings.A molecule can be involved in many future events but only the earliest one is registered as the candidate event for the molecule.Thus, the number of candidate events is always equal to the number of simulated molecules.The goals are determination and execution of the earliest event.Following the execution, candidate events of related molecules are invalidated and possible new candidate events are calculated.This process repeats throughout the simulation.
Determination of earliest event by linear search has a complexity of () and makes EDMD simulations practically impossible for large number of molecules.In this study, a priority queue method suggested by Paul [11] was implemented.In his method, all events are distributed among time interval bins and events belonging to earliest bin are sorted using a complete binary tree (CBT) structure.This approach reduces computational effort approximately to (1) on average.
Macroscopic thermodynamic properties of a simulated system can be estimated by ensemble averaging over particles' behaviour.However, this would be expensive in terms of computational effort because of the requirement of large number of molecules to obtain smooth profiles.The alternative approach is the time averaging for a smaller number of particles since ensemble average and time average of a phase variable yield equal values in an ergodic system [23].Since this study deals with large number of molecules, an ensemble averaging scheme was performed by averaging a series of snapshots of computational domain taken in equal intervals.For the sake of profiling, computational domain was partitioned into small subdomains (bins).Each molecule is associated with a certain bin based on its molecular centre.Although the selection of bin size is arbitrary, the average number of molecules inside a bin is an important parameter for an accurate profiling of system properties.Small number of molecules may demonstrate large fluctuations due to statistical nature of EDMD simulations.On the other hand, resolution of distribution is reduced for larger bins.In order to satisfy both criteria, computational domain was divided into 51 spatial bins along a direction (51 × 51 for 2D distributions) in this study.Snapshotting was started after at least 100 collisions per particle (CPP) with a frequency of at least 1 CPP.In order to obtain smoother profiles, longer CPP values were used for averaging.

Implicit Treatment for Flow Boundaries.
From microscopic viewpoint, in addition to bulk velocity, gas molecules are translated by thermal velocity, also.In high speed (hypersonic) flows, the magnitude of the bulk velocity is higher than that of the thermal velocity.For such flows, conventional approach is the use of vacuum boundary at the exit of the channel (when a molecule reaches the vacuum boundary it leaves the computational domain permanently and there is no inward flux at that boundary).
In contrast to hypersonic flows, the magnitude of thermal motion is comparable with bulk velocity in low speed flow and inward flux because thermal motion of molecules should be taken into consideration.In order to cover a wide range of flow regimes, Liou and Fang [21] introduced an implicit treatment for flow boundaries in DSMC simulations.This treatment was adapted for EDMD simulations in this study.
In this approach, the number of incoming molecules and their corresponding velocities are estimated implicitly by the local flow properties.For either upstream or downstream boundary, molecular flux () into the computational domain can be determined using the Maxwellian distribution function: where Flux is calculated for each cell face of flow boundary cell .
Here, streamwise velocity and local temperature are denoted Mathematical Problems in Engineering by  and , respectively.Value of  is 0 for upstream and  for downstream boundary.  is the number density of molecules in the cell and  mp is the most probable speed.For a time interval , the number of molecules (  ) entering domain from surface area   of cell  is Incoming molecules are introduced at random positions within the corresponding cells, touching the outside of the cell surface.In case of any overlap (i.e., centre-to-centre distance of the molecules being less than sum of the radii) molecule is repositioned.Tangential velocity components (V and ) of an incoming molecule are independent of streamwise velocity and generated as follows: Here,  and  are mean tangential velocity components.R  and R   are random numbers generated from normal distribution with zero mean and unit variance.
If streamwise inlet velocity is zero ( = 0) normal component is generated as follows: Here, R  is a random number generated from uniform distribution of interval [0, 1).For nonzero streamwise inlet velocities ( ̸ = 0), Garcia and Wagner [24] introduced several efficient acceptance-rejection methods with the general form of Here, (/ mp ) is a random number selected through an acceptance-rejection method.In this study, the recommended method for low speed flows (−0.4 mp <  < 1.3 mp ) [24] was used.Flow properties of both upstream and downstream boundaries must be determined in order to be used in ( 1)- (6).For upstream boundary, temperature and pressure are known.The number density,  in , is obtained from ideal gas equation: Transverse velocities,  in and  in , are set to zero.Streamwise velocity, ( in )  , is extrapolated from neighbour boundary cell  where  and  are density and speed of sound, respectively [22]:

Verifications
Shear-Driven and Pressure-Driven flows in confined microchannels are well-studied problems.They are frequently used as test cases for the robustness and validations of flow simulations.In this section, the theoretical background and the simulation setup for such flows are presented.A new implementation of EDMD simulation including cell partitioning, event scheduling, and implicit boundary treatment was developed in C# with serial processing.Benefiting from object oriented programming (OOP), code can be further extended to cover different molecular types and boundary treatments.On the other hand, considering the computational performance the parallelization of the code is an essential issue in molecular simulations.However, this deserves a detailed study of its own.All simulations were carried out on a fourth-generation Intel i7 4700 CPU computer with 16 GB of RAM.
3.1.Shear-Driven (Couette) Flow.Shear-Driven flow in a confined channel is obtained by moving the top plate at a constant velocity (Figure 1(a)).Pressure gradient is zero between the upstream and downstream boundaries.Initial temperature of the argon inside the channel and at the upstream boundary is 300 K. Upper wall moves with a velocity of 100 m/s while lower wall is stationary.Periodicity is applied along axis.Spacing ratios (the ratio of average distance between molecules, , to average diameter of molecules, ) for this study were selected as 4.0 and 7.0.Flow configurations were chosen according to Knudsen numbers of 0.1, 0.5, and 1.0 in three cases (Table 1).The number of simulated particles was around 200,000 for each case.Simulation results were compared with the continuumbased model in which velocity gradient between two walls of the channel is affected by the Knudsen number and has a linear nature in Shear-Driven flows.For a fully diffusive wall, normalized velocity profile is shown as follows, where   is wall velocity and  is channel height [25]:  Boundary conditions Fully diffusive walls @300 K along -direction Periodicity along -direction 3.2.Pressure-Driven Flows.Pressure-Driven flow is generated by setting different pressure levels on upstream and downstream boundaries (Figure 1(b)).Two Pressure-Driven flows of different pressure ratios (PR =   /  ) in a microchannel were investigated.Flow properties are tabulated in Table 2.The number of simulated molecules was about 300,000 for both cases.
For Pressure-Driven flows, Arkilic et al. defined the pressure distribution along the fully diffusive channel as a function of Knudsen number at the exit of the channel (Kn e ) of length  as [26]  () For fully diffusive walls, analytical solution of Navier-Stokes equations with second-order slip boundary condition yields a normalized velocity profile [27]: where   is centreline velocity.Slip-to-centreline velocity ratio along the channel then becomes For second-order slip boundary condition, velocity profile and corresponding slip-to-centreline ratio are If slip velocity is normalized with the local average velocity, ( 12) and ( 15) become [28]

Results and Discussions
In this section, we present simulation results of Shear-Driven and Pressure-Driven flows for aforementioned cases.
Results were compared with analytical solutions derived from Navier-Stokes equations with slip flow boundary conditions and DSMC results.Simulation configuration (number of simulated molecules, spacing ratio, boundary types, etc.) greatly affects computational speed.As an example, simulation of one million argon molecules in a cubical region when the spacing ratio is 5.0 takes 450 seconds to reach one collision per particle on the aforementioned PC configuration.

Shear-Driven Flow.
In Cases A-C, six Couette flow configurations at different Knudsen numbers were simulated.
Comparison of the predicted slip velocities between simulation results and the analytical solution is given in Table 3.On the upper and lower walls, slip velocities were estimated by linearly extrapolating average velocities of neighbour bins.Knudsen number was 0.1 in Case A and the flow was in the slip-to-transition regime.Consequently, deviation from the analytical solution is reasonably small, less than 10%.In Cases B and C, flows were in the transitional regime (0.1 < Kn ≤ 10).Therefore, simulated slip velocities were predicted noticeably lower than the analytical solution.Deviation is between 21.3% and 25.6% in these cases.
Predicted velocity profiles and deviations from the analytical solution are presented in Figures 2-4.For low Kn flow (Case A), velocity profile is almost linear and corresponding deviation is insignificant.As rarefaction intensifies (Cases B and C), velocity profile deviates particularly at the near wall regions.S shaped velocity profile is more apparent in Case C.
Another finding is that nearly identical velocity profiles were obtained for different spacing ratios at the same Knudsen number.As the spacing ratio reduces from 7.0 to 4.0, computational speed quadruples.Utilization of this piece of information can ease the computational cost of larger systems; that is, EDMD performs faster in smaller spacing ratios [9].respectively.The results are compared with linear pressure distribution and slip formulation given by Arkilic et al. [26].Figures 5(b) and 6(b) show deviation between MD simulation results and analytical solutions.Relative pressure difference is less than 5% for both cases.This indicates that adapted boundary treatment is appropriate for EDMD simulations.
Deviation from the linear pressure distribution is given in Figures 5(c) and 6(c).In both cases, nonlinearity is slightly lower in EDMD predictions.
Local mean free path increases along the channel due to pressure gradient.This yields an increase in local Knudsen number (Figure 7).In Case D, Kn was between 0.080 and  0.211 (i.e., the flow was in slip-to-transition regime).In Case E, flow was in the transitional regime with corresponding Kn up to 1.2.Ratio of the slip velocity to the local average velocity is given in Figure 8. Knudsen number values in the abscissa correspond to local values formed in Figure 7 along the channel.Slip velocities were linearly extrapolated from outermost bins similar to Couette flow calculations.Simulation results were compared with analytical solutions with first-and secondorder slip boundary conditions (16), linearized Boltzmann solution, and DSMC results [28].It is seen that EDMD results are in good agreement with DSMC and linearized Boltzmann  solution.On the other hand, breakdown of slip boundary models starts around Kn = 0.1.Both first-and secondorder models overestimate the velocity slip.Note that, for both Cases D and E, upstream and downstream slip velocity ratios are smaller than in other regions of the computational domain.Use of Maxwell distribution to sample number density and velocities in pressure boundary treatment is the main reason of this behaviour.Such implementation of boundary treatment works better when the rarefaction effects are small.For highly rarefied (high Kn) flows, Ahmadzadegan et al. [29] advise using Chapman-Enskog velocity distribution to improve pressure boundary treatment.Normalized velocity profiles (/  where   is local centreline velocity) in cross-stream direction at three locations along the channel are shown in Figures 9 and 10 for Case D and Case E, respectively.Note that flow velocity along the channel increases as pressure decreases.Thus, for a given position along the channel, normalization was carried out by corresponding centreline velocity.Here, results of EDMD simulations were compared with second-order analytical solution.In Case D, where Kn is small, curves agree well whereas in Case E estimation of higher slip velocity in analytical solution becomes more apparent as Kn increases.

Conclusion
This paper presents an event-driven molecular dynamics (EDMD) simulation of monatomic gas flows in microchannels.In this study, molecules are considered as hard-spheres with discrete interaction potentials and interaction times are analytically predictable.Therefore, each collision is real and takes place in its calculated position.Postcollisional velocities are also analytically calculated according to the conservation of momentum and energy.By its deterministic nature, EDMD differs from DSMC method in which simulated system is represented by a smaller group of molecules and collision pairs are selected stochastically.Additionally, by using cell partitioning technique and use of priority queues for event sorting, EDMD is now a strong alternative to DSMC method in terms of not only accuracy but also computational speed.
In this study, an implicit treatment for flow boundaries is adapted from DSMC method to EDMD for the first time in order to induce Shear-and Pressure-Driven flows in confined channels.The number of introduced molecules through upstream and downstream boundaries and their corresponding velocities and locations are determined implicitly from local number density, mean flow velocity, and temperature.
Velocity slip in Pressure-Driven flows is compared with DSMC results and solution of linearized Boltzmann equation.The good agreement confirms the robustness of the treatment and its applicability to EDMD simulations.For both Shearand Pressure-Driven flows, velocity slip is also compared with the one derived from Navier-Stokes equations of slip boundary conditions.For low Knudsen number (Kn) flows, velocity profiles and calculated slip velocities agree well.But as the flow rarefies, slip boundary models estimate higher slip velocities.
Additionally, simulations of two flows having the same Kn but different spacing ratios produce almost identical velocity profiles in Shear-Driven cases for the investigated range of Kn.Since simulations with lower spacing ratios run faster, this finding can be further extended to scale rarefied flows and to reduce the computational time.

Figure 5 :
Figure 5: Pressure distribution in Pressure-Driven flow, Case D: (a) comparison of pressure distribution along -axis; (b) deviation between simulation and N-S results; (c) deviation from linear pressure distribution.

Figure 6 :
Figure 6: Pressure distribution in Pressure-Driven flow, Case E: (a) comparison of pressure distribution along -axis; (b) deviation between simulation and N-S results; (c) deviation from linear pressure distribution.

Figure 7 :
Figure 7: Local Kn distribution in Pressure-Driven flow.

Figure 8 :
Figure 8: Normalized slip velocity for slip-to-transition region.

Figure 9 :Figure 10 :
Figure 9: Velocity profiles in -direction of Pressure-Driven flow, Case D.

Table 2 :
Flow configurations in Pressure-Driven flow.

Table 3 :
Comparison of slip velocities.