Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 842837, 12 pages

http://dx.doi.org/10.1155/2015/842837

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

^{1}Department of Energy Systems Engineering, Muğla Sıtkı Koçman University, 48100 Muğla, Turkey^{2}Department of Mechanical Engineering, Gebze Technical University, Gebze, 41400 Kocaeli, Turkey

Received 20 July 2015; Accepted 1 December 2015

Academic Editor: Manfred Krafczyk

Copyright © 2015 Volkan Ramazan Akkaya and Ilyas Kandemir. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Classical solution of Navier-Stokes equations with 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 high Knudsen numbers (Kn) seen at dilute gases or narrow channels. In this study, an event-driven molecular dynamics (EDMD) simulation is proposed to estimate properties of hard-sphere gas flows. Considering molecules 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.

#### 1. 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, (the ratio of the mean free path of gas molecules, , to the characteristic length of the channel, ), increases. Magnitude of determines the flow regime: continuum regime for , slip regime for , transitional regime for , and free molecular flow for . In many macroscale applications, is negligible compared to the channel size 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–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 micro- or 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–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 event-driven 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 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.

#### 2. 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 was selected as the monatomic gas of this study. In order to reduce numerical errors, all magnitudes were nondimensionalized by considering unit Boltzmann constant () and scaling with the reference variables temperature , mass ( kg), and diameter ( 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 . 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, molecule-boundary 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 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.

##### 2.2. 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:whereFlux is calculated for each cell face of flow boundary cell . Here, streamwise velocity and local temperature are denoted by and , respectively. Value of is 0 for upstream and for downstream boundary. is the number density of molecules in the cell and is the most probable speed. For a time interval , the number of molecules () entering domain from surface area of cell isIncoming 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 ( and ) of an incoming molecule are independent of streamwise velocity and generated as follows:Here, and are mean tangential velocity components. and are random numbers generated from normal distribution with zero mean and unit variance.

If streamwise inlet velocity is zero normal component is generated as follows:Here, is a random number generated from uniform distribution of interval . For nonzero streamwise inlet velocities , Garcia and Wagner [24] introduced several efficient acceptance-rejection methods with the general form ofHere, is a random number selected through an acceptance-rejection method. In this study, the recommended method for low speed flows [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, , is obtained from ideal gas equation:Transverse velocities, and , are set to zero. Streamwise velocity, , is extrapolated from neighbour boundary cell where and are density and speed of sound, respectively [22]:At downstream boundary only pressure is known. Thus, number density , temperature , and velocity components , , and are extrapolated from neighbour cell:

#### 3. 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.