- About this Journal ·
- Abstracting and Indexing ·
- Advance Access ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Table of Contents

Advances in Mechanical Engineering

Volume 2012 (2012), Article ID 580763, 10 pages

http://dx.doi.org/10.1155/2012/580763

## Molecular Dynamics Simulations of Pressure-Driven Flows and Comparison with Acceleration-Driven Flows

Laboratoire Modélisation et Simulation Multi Echelle, UMR-CNRS 8208, Université Paris-Est, 5 boulevard Descartes, 77454 Marne-la-Vallée Cedex 2, France

Received 1 July 2011; Revised 6 December 2011; Accepted 15 December 2011

Academic Editor: Yogesh Jaluria

Copyright © 2012 Quy Dong To et al. 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

We use molecular dynamics to simulate fluid flows between two parallel plates with constant wall temperature. Unlike the usual approach in molecular dynamics, instead of applying an external force on the molecules, the periodic boundary conditions are modified to create a pressure difference between the inlet and the outlet sections of the computational domain. The simulation results include velocity, pressure, density, and temperature profiles obtained by the new method. These results are compared with approximate solutions for nonisothermal Poiseuille flows. The method is also applied to simulate a flow in a rib-roughened channel.

#### 1. Introduction

Over the past decades, nano/microfluidic systems have seen a rapid development and contribute to a revolution in many areas such as biotechnology and medicines. At the same time, calculation methods are also needed to analyze the reliability of the systems. The main challenge of the nano/microfluidics computations is that the studied system has a considerable surface/volume ratio so that the surface phenomena become highly important. These size effects cannot be solved by the traditional computational method for fluid dynamics.

Depending on the nature of the problem, there are different suitable computational methods, such as those based on the Navier-Stokes and energy equations, the Burnett or super-Burnett models [1–3], and the direct simulation Monte Carlo methods [4]. For nanoflows, the molecular dynamics (MD) method is one of the most accurate methods since realistic interactions between particles or between particles and boundaries may be accounted for.

One of the fundamental and practical problems in fluid mechanics is the steady-state flow between two parallel plates, the plane Poiseuille flow. Under the incompressibility assumption, the problem has a straightforward analytical solution. It is also studied extensively by other methods (Karniadakis et al. [5]). Most of the molecular dynamic simulations reported in the literature ([6–9], e.g.) concern flows driven by acceleration body force, rather than pressure driven flow. These limitations are due to the usage of the periodic boundary conditions which cannot generate pressure gradient inside the simulation domain. Several authors (Lupkowski and Van Swol [10], Li et al. [11], and Sun and Ebner [12]) have developed new techniques to mimic flows induced by pressure difference. Lupkowski and Van Swol [10] placed two rigid walls at the inlet and outlet and applied external forces on them. Li et al. [11] used a fictitious membrane which allows atoms to pass from one direction and forces atoms from the other direction to be elastically reflected with a given probability. Sun and Ebner [12] simulated 2D fluid flows in a long box while controlling the temperature at the two ends.

The present paper aims at presenting a new molecular dynamics approach based on pressure differences. Section 2 gives the problem statement and presents some approximate analytical solutions reported in the archival literature for Poiseuille flows. The pressure expression and the pressure difference from the atomistic point of view are also discussed. Section 3 is focused on the implementation of periodic boundary conditions linked to the pressure difference concept and its implementation in an MD code. The results of simulations are then compared with approximate analytical solutions of the Navier-Stokes and energy equations in Section 4. The application of the present method to fluids other than ideal gases and to rib-roughened channels is also discussed. Finally, some conclusions and perspectives are given in Section 5.

#### 2. Problem Description and Theoretical Background

##### 2.1. Problem Description

The objective of the present work is to study the fluid flows through a nanochannel. In the cartesian coordinate system , we assume that the flow direction is along the -axis. The length, width, and height of the channel are denoted by , , and in directions , , and , respectively. The channel width being very large with respect to the other dimensions, it will be treated as infinite. To produce a fluid flow from the inlet to the outlet, we keep the inlet pressure at higher level than the outlet . The fluid can exchange energy with the channel walls which are maintained at constant temperature, .

The Knudsen number, Kn, is a quantity to identify the fluid flow regime: from continuum (small Kn) to molecular flow regime (high Kn). The range of chosen in the present study, , belongs to the slip flow regime according to [5], that is, the fluid continuum mechanics is still useful provided that slip boundary conditions are applied at the walls. The validity of the Navier-Stokes equations, as shown by many authors (see [7], e.g.), can therefore be extended to Poiseuille and Couette flows for the present range of Kn. The Knudsen number is defined as where is the mean free path, and is a characteristic scale of the channel. In this work, is estimated by The quantities and in (2) are the density number and the diameter of the fluid molecules, respectively. For Lennard Jones (LJ) fluids, the molecular diameter value is approximately taken equal to the parameter appearing in the LJ potential formula introduced in Section 4.1.

In what follows, the terms “pressure-driven flow” and “acceleration-driven flow” are employed to define two particular cases: (i) the body force applied on a unit fluid volume is zero, and the pressure decreases along the flow direction, , ; (ii) the inverse situation where , . We shall also consider some available analytical solutions of the Navier-Stokes equations reported in the literature. For both flows, the fluid velocity in a cross-section which accounts for slip effects (Navier boundary conditions) at the wall admits the parabolic form where is a positive dimensionless constant related to pressure gradient (or body force), viscosity, slip length, and so forth. It is noted that (3) is generally obtained for incompressible fluids, but it is also valid for compressible fluid flows in a long channel [13]. In the latter case, and are functions of the streamwise coordinate .

An approximation solution of the energy equation can be obtained by neglecting the convective term [14]. The temperature becomes then a quartic function of the coordinate , Equation (4) is based on the incompressibility assumption. For compressible fluid flows in a long channel, Cai et al. [15] used a perturbation technique to derive the following temperature distribution: The dimensionless coefficients in (4) and (5) are related to pressure gradient (or body force), viscosity, conductivity, specific heat, and so forth. In (5), the coefficients depend on the position . We also note that for acceleration-driven flow, Todd and Evans [6] suggested a correction to the temperature profile given by (4) as a sextic function of , In Section 4, MD simulations are used to reexamine the validity of the velocity and temperature profiles given by (3)–(6) and to determine the coefficients , , , and by curve fitting.

##### 2.2. Pressure Tensor and Pressure Difference

Before proceeding the molecular dynamics simulation, let us look at the definition of pressure tensor from the atomistic point of view. Using statistical mechanic theory, Irving and Kirkwood [16] derived the following pressure tensor decomposition (IK): The kinetic term depends on the molecular square velocity, while the potential term depends on pairwise interactions between molecules. In a system of molecules, and read The term inside the angular bracket denotes the ensemble average, and denotes the dyadic product. The terms , , and are respectively the position vector, the velocity of the particle , and the mean velocity. The distance vector and interaction force vector between two molecules and are denoted, respectively, as and . For LJ fluids, the two vectors and are colinear, and the force magnitude is derived from the interaction potential in function of the distance between the two molecules In (8), the notation is used for the delta Dirac function at , and the expression of the operator reads Equations (10) and (8) show that the derivation of IK pressure tensor involves an infinite sum of high-order derivatives of the delta function and ensemble average, not suitable for MD computations. A more convenient form of the pressure tensor and an associated calculation method, the method of plane (MOP), was proposed by Todd et al. [17] and by Evans and Morriss [18]. When the fluid density is uniform, the operator is reduced to unity, and (8) becomes (see [17, 19–22]) where is the volume of the fluid element located at . In the literature, (11) is sometimes referred to as IK1 pressure [17, 18]. In (8) and (11), the pressure can be decomposed into two scalar quantities and as follows: It is clear that for low-density fluids, high temperature, and small interaction energy , the pressure is reduced to the kinetic part , as for ideal gases. For the other fluids, the potential part cannot be neglected. The pressure difference between two points at distance reads The pressure components , , , number density , and stream velocity are invariant in direction . From a microscopic point of view, taking , we must have As it will be discussed in the following, all these conditions are satisfied in our MD simulation algorithm if we apply periodicity for the boundary conditions in the direction. However, , , and can all vary in the flow direction . We note that in the steady-state regime, all quantities like , , , and are stable with time, and thus their differences , , , and are also stable. In particular, (13) shows that the finite variation of temperature along is related to those of other quantities via the expression If the stream velocity is much smaller than the thermal velocity, we can obtain a simplified expression for the temperature difference as The objective of the MD method discussed in what follows is to maintain the differences in squared velocity and, thus in temperature, in direction . By that way, we can indirectly generate the pressure difference . In the case where density is invariant in -direction (incompressible flow assumption), the temperature difference is proportional, for an ideal gas, to the pressure difference since

#### 3. Numerical Procedure

In molecular dynamics, periodic boundary conditions (PBC) applied to velocities are powerful techniques to reduce the study of a large system to a smaller one far from the edge. Considering a simulation domain as a cube, PBC requires that if a molecule passes through one face, it reappears on the opposite face with the same velocity. Obviously, there is no difference in pressure, density, or temperature between any two opposite faces of the simulation domain. Consequently , as it will be shown in the next section.

##### 3.1. Modified Periodical Boundary Conditions

In our problem, all the virial pressure components , , and do not vary along , which can be satisfied by the traditional PBC applied on the faces . However, to create a pressure difference along the -axis, we must develop a strategy that produces a constant difference between the squared velocity of the molecules crossing the face and those crossing the face .

Usually, to keep unchanged the total number of molecules in the domain whenever a molecule goes out of the domain, we must insert another molecule inside. Thus, if we relax the PBC constraints, there are many ways to assign velocities to ingoing molecules, either independent of or dependent on the outgoing ones. In this paper, we generalize the PBC to account for the pressure difference by maintaining the difference in squared velocity at the two opposite faces and .

From (15) and (16), we know that a constant pressure difference is related to a constant difference in squared velocity. A constant parameter is thus used to create a constant difference in squared velocity between the inlet and outlet faces (see Figure 1). In order to apply a velocity square difference equal to a constant between the faces and , we modify the periodic boundary conditions so that if a molecule goes through the face with a velocity , we insert a molecule at with a velocity satisfying the following conditions: Analogously, if a molecule crosses the face with a velocity , it reappears at with a velocity so that Since is positive, it may happen that . It is then impossible to find satisfying (19). The boundary conditions applied to the simulation domain are summarized in Table 1. However, this effect can be assumed negligible for flows with relatively small speed (of the order of the thermal speed), as seen in most MEMS/NEMS devices. In the present simulations, the number of molecules that do not satisfy (19) were found less than of the outgoing molecules. In order to simulate low-speed flows, small values of were used.

#### 4. Results

##### 4.1. Molecular Dynamic Simulations of Ideal Gas Flows

To model the interaction between the molecules, the 6–12 Lennard Jones potential is used where is the depth of the potential well, is the finite distance at which the interparticle potential is zero, and is the distance between the particles.

The potential well depth parameter is with being the reference potential well depth. Two cutoff distances were used in this study: [17] corresponds to repulsive interactions between molecules, while [21, 23–25] corresponds to attractive-repulsive interactions. In general, the present method is valid for any pair potential. However, the simulation conditions of Table 2, that is, , lead to a very small potential part of the pressure (less than 1% of the kinetic part ) making the behavior of the fluid close to an ideal gas.

The two walls are modeled as thermally diffusive walls at the same and constant temperature. Whenever a molecule arrives at a wall, it is reflected with the velocity corresponding to the wall temperature and with a random direction. In this work, we do not use thermostats, and the fluid can exchange energy with the wall due to the wall model used here.

In our simulations, the global number density is kept fixed, in every configurations. Other geometric parameters like , , , and the number of the molecules are changed as shown in Table 2. The values of are chosen such that the difference of the square of the molecular velocity between the inlet and outlet is small enough so that the stream velocity is much smaller than the thermal speed. From the microscopic viewpoint, the local stream velocity is the local average velocity of the molecules, and from the macroscopic viewpoint, it is the flow speed in the Navier-Stokes equations and thus connected to the pressure gradient (or volume force ). On the other hand, the thermal speed is the root mean square of the molecular velocity and is directly connected to the local temperature . The slope of the simulated total pressure distribution along the channel (see Figures 2 and 3) is used to determine the pressure gradient and the equivalent external acceleration used in the acceleration-driven flow simulations. The quantity is the average mass density of the system, . The computation with Leap-Frog Verlet integration scheme is carried out for time steps from equilibrium, each of which is equal to 0.005 unit time . The height and length of the channel are divided into 50 layers to determine accurately the distribution of local velocity , temperature , and pressure tensor . The local velocity measured in one bin located at is given by being the particle number in the bin. The local pressure (IK1-model) is The pressure components were also computed by using the method of plane [17, 26]. With , being , , or , the pressure component along direction and acting on the area element normal to the axis is defined by In (23), the first sum is for all molecules crossing the area element over the period , and the second sum is for all pairs whose distance vectors cut the area element . The temperature of one bin is then calculated by In what follows, the -distributions of pressure, temperature, and density correspond to an average on and of the local quantities and . They are denoted as and , respectively. Analogously, the profile along of velocity, temperature, and density corresponds to an average on and of local quantities with the notations: , , and , respectively.

From the simulations, the axial pressure distribution is plotted in Figure 2 for the pressure-driven and the acceleration-driven flows with , . As expected, for the acceleration-driven flow case, the pressure is constant. In contrast, the pressure distribution for pressure-driven flow decreases linearly. We note here that the slope of the pressure curve, equal to , is used to compute the value in the acceleration flow simulation according to the formula . Under the simulation conditions , , and , , the kinetic part of the pressure dominates.

In order to know how the pressure difference is distributed in the -, -, -directions, the pressure components , , and along are plotted in Figure 3. All the three pressure components decrease with , and the slopes are almost the same. The values of components , are very close, whereas the value of component is smaller because the fluid is confined in the -direction. Although the difference between the pressure components is small (less than 4%), it implies that the velocity distribution deviates from the equilibrium distribution. Because the pressure is the average value of , , , the -curve lies between the others and has the same slope. It should be noted that (22) and (23) lead to relatively close results when considering the average axial pressure component (see Figure 4). Figure 5 shows that the fluid density distribution across the channel width varies slightly around the average density , except near the walls where it can be as high as . However, the phase diagram for LJ fluid [25] shows that the fluid near the wall is still in gas state. In our simulations, and are chosen to be very small (, ), and the fluid temperature near the wall is relatively large (the temperature of the wall is ). The fluctuation of the density profile shows that the incompressibility assumption is no longer valid, especially near the channel walls. However, the fluctuations being localized, it is interesting to check if the analytical solutions presented in Section 2 still agree with the MD solutions.

Figure 6 shows that acceleration and pressure-driven flows exhibit parabola-like velocity profiles in the center region of the channel, in agreement with (3). Near the wall where the Knudsen layer dominates, the velocity distribution tends to deviate from the solution given by (3), and a velocity slip at the wall is predicted. Based on the kinetic theories, different models have been derived to explain the slippage [5], and generally, the slip effect becomes important when the Knudsen number increases. The little difference in velocity profile between the two types of flow shows that, despite their different microscopic natures, volume forces can be seen as equivalent to pressure gradients at the macroscopic scale.

About the temperature profile across the channel width, there is no visible difference between acceleration and pressure-driven formulations whatever Kn (Figure 7). The temperature profile cannot be fitted by any of the three approximate analytical solutions ((4)–(6)). The temperature is minimal at the center of the channel and increases rapidly near the wall, which does not agree with the approximate solutions obtained with the incompressible assumption ((4)-(5)). With the simulation conditions reported in Table 2, this reverse trend cannot thus be explained by using the incompressible flow equations. This anomaly was also observed by using DSMC method [27] and super-Burnett (SB) method [28]. Note that the flow speeds in our simulations are very small, as in [27, 28].

Next, we study the influence of and on the temperature profile. By increasing the acceleration parameter or the -parameter (squared velocity), the stream velocity is increased. A change in the form of the temperature profile is observed in Figures 8 and 9 as (or ) increases. At high values of (or ), the temperature profiles seem to be closer to the approximate solutions which predict a maximum at the center of the channel and a minimum at the walls. At small values of (or ), the temperature profiles for the two flow types differ just a little. However, at high values of (or ), the discrepancy becomes more important: the curvature at the center of channel for the pressure-driven flow case is higher than that for the acceleration-driven flow case. Using a curve-fitting procedure, we find that the temperature profile of the pressure-driven flow case agrees quite well with the quartic expression (5). The temperature for the acceleration-driven flow case does not fit well with (4). To explain this discrepancy, Todd and Evans [6] argued that the transport coefficients are not constant and that an additional cross-coupling between strain rate and heat flux exists. They proposed a correction to (5) with a -extraterm for the temperature profile (see (6)). The dashed line drawn in Figure 9 shows that the sextic polynomial fits very well the present numerical results.

##### 4.2. Applications to More General Cases and Rib-Roughened Channels

The application of the method to cases for which the fluid particles interact strongly, that is, the potential pressure is of the order of , is considered in that section. The case of rib-roughened channels is also briefly discussed.

With the algorithm developed in the previous section, two cases were considered with the same parameters, except for the cutoff distance and density number . In the first case, , and the cutoff distance is set to (attractive-repulsive interaction). In the second case, the fluid density is increased up to while setting as in the previous section. The well depths are the same in both cases and equal to the reference value, that is, . All the parameters for these two cases are summarized in Table 3.

In the first case (), Figure 10 displays the kinetic part and the virial pressure according to (11) along the channel. Except a slight curvature in the curve, both pressure curves decrease continuously and linearly. The relative difference between the slopes of both curves is less than 1%, which means that the gradient of the potential part is also small in comparison with the kinetic part even if is not negligible (). In the second case (), Figure 11 shows a totally reverse trend: both the virial pressure profile and its potential part decrease with . The variation of the potential part contributes to of the virial pressure variation with and . We conclude that the method is convenient for generating virial pressure profiles which decrease in the axial direction.

Next, we consider the velocity profiles across the channel section for these two cases. The velocity profiles (Figure 12) agree then with (3). The slip length for the case is smaller than for the case because the mean free path is decreased when the density number is increased. When computing the average pressure component along the channel, Figure 13 shows that there is no visible difference between (22) and (23).

It should thus be emphasized that the present method is relevant to generate various fluid flows even if we control only the difference between the inlet- and outlet-squared molecular velocities. The potential part , which depends on the interatomic interaction, can vary along the flow direction. As shown numerically, the variation of contributes considerably to the pressure gradient.

The third case aimed at modeling the flow within a rib-roughened channel. The parameters are the same as in the second case (, ), but the channel height is suddenly reduced in its middle part by inserting on both walls ribs of height and length . The rib-to-channel width ratio is equal to and (Figure 14). Figure 14 exhibits the onset of two vortices close to the rib corners at the upstream section. It is also shown that the fluid flow is highly nonuniform and characterized by wavy streamline patterns within the downstream region. From Figure 15, the axial pressure variations are rather different from those predicted for smooth channels. Although the variations of the axial pressure predicted by (22) and (23) are quite similar, there are considerable differences between them. The MOP pressure is smaller than the IK1 pressure. The changes in pressure at the downstream and upstream sections display also smoother pressure variations. This considerable difference is due to the strong inhomogeneity of the fluid within the channel. Generally, the MOP, based on the hydrodynamics equation, must be used in such situation.

#### 5. Conclusions

In this work, we have used a molecular dynamics method for simulating pressure-driven flows in channels. The main novelty lies in the modification of the periodic inlet and outlet velocity conditions: the difference in squared velocity between the ingoing and outgoing molecules of the simulation domain is maintained constant.

Plane Poiseuille flows were simulated, and the results were compared with approximate analytical solutions of the Navier-Stokes and energy equations reported in the literature. This new method has been proven to be realistic and could be applied to many problems with a nonconstant axial pressure gradient: flows around obstacles or in rough channels, for example. The method appears thus relevant for modeling compressibility effects as well as temperature variations along the flow direction, a domain unfulfilled by the usual molecular dynamics methods.

#### References

- H. Struchtrup,
*Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory*, Springer, New York, NY, USA, 2005. - W. W. Liou and Y. Fang,
*Microfluid Mechanics*, McGraw-Hill, New York, NY, USA, 2003. - S. Chapman, T. G. Cowling, and D. Burnett,
*The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases*, Cambridge University Press, Cambridge, UK, 1990. - G. A. Bird,
*Molecular Gas Dynamics and the Direct Simulation of Gas Flows*, Clarendon Press, 1994. - G. Karniadakis, A. Beskok, and N. Aluru,
*Microflows and Nanoflows: Fundamentals and Simulation*, Springer, New York, NY, USA, 2005. - B. D. Todd and D. J. Evans, “Temperature profile for Poiseuille flow,”
*Physical Review E*, vol. 55, no. 3, pp. 2800–2807, 1997. View at Scopus - D. K. Bhattacharya and G. C. Lie, “Molecular-dynamics simulations of nonequilibrium heat and momentum transport in very dilute gases,”
*Physical Review Letters*, vol. 62, no. 8, pp. 897–900, 1989. View at Publisher · View at Google Scholar · View at Scopus - J. L. Barrat and L. Bocquet, “Large slip effect at a nonwetting fluid-solid interface,”
*Physical Review Letters*, vol. 82, no. 23, pp. 4671–4674, 1999. View at Scopus - Q. D. To, C. Bercegeay, G. Lauriat, C. Léonard, and G. Bonnet, “A slip model for micro/nano gas flows induced by body forces,”
*Microfluidics and Nanofluidics*, vol. 8, no. 3, pp. 417–422, 2010. View at Publisher · View at Google Scholar · View at Scopus - M. Lupkowski and F. Van Swol, “Computer simulation of fluids interacting with fluctuating walls,”
*The Journal of Chemical Physics*, vol. 93, no. 1, pp. 737–745, 1990. View at Scopus - J. Li, D. Liao, and S. Yip, “Coupling continuum to molecular-dynamics simulation: reflecting particle method and the field estimator,”
*Physical Review E*, vol. 57, no. 6, pp. 7259–7267, 1998. View at Scopus - M. Sun and C. Ebner, “Molecular-dynamics simulation of compressible fluid flow in two-dimensional channels,”
*Physical Review A*, vol. 46, no. 8, pp. 4813–4818, 1992. View at Publisher · View at Google Scholar · View at Scopus - E. B. Arkilic, M. A. Schmidt, and K. S. Breuer, “Gaseous slip flow in long microchannels,”
*Journal of Microelectromechanical Systems*, vol. 6, no. 2, pp. 167–178, 1997. View at Scopus - L. D. Landau and E. M. Lifshitz,
*Fluid Mechanics*, Pergamon Press, 1987. - C. Cai, Q. Sun, and I. D. Boyd, “Gas flows in microchannels and microtubes,”
*Journal of Fluid Mechanics*, vol. 589, pp. 305–314, 2007. View at Publisher · View at Google Scholar · View at Scopus - J. H. Irving and J. G. Kirkwood, “The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics,”
*The Journal of Chemical Physics*, vol. 18, no. 6, pp. 817–829, 1950. View at Scopus - B. D. Todd, D. J. Evans, and P. J. Daivis, “Pressure tensor for inhomogeneous fluids,”
*Physical Review E*, vol. 52, no. 2, pp. 1627–1638, 1995. View at Publisher · View at Google Scholar · View at Scopus - D. Evans and G. Morriss,
*Statistical Mechanics of Nonequilibrium Liquids*, Cambridge University Press, Cambridge, UK, 2008. - W. G. Hoover, C. G. Hoover, and J. F. Lutsko, “Microscopic and macroscopic stress with gravitational and rotational forces,”
*Physical Review E*, vol. 79, no. 3, Article ID 036709, 2009. View at Publisher · View at Google Scholar · View at Scopus - D. H. Tsai, “The virial theorem and stress calculation in molecular dynamics,”
*The Journal of Chemical Physics*, vol. 70, no. 3, pp. 1375–1382, 1979. View at Scopus - M. P. Allen and D. J. Tildesley,
*Computer Simulation of Liquids*, Oxford University Press, Oxford, UK, 1989. - J. P. Hansen and I. R. McDonald,
*Theory of Simple Liquids*, Academic Press, New York, NY, USA, 2006. - D. C. Rapaport,
*The Art of Molecular Dynamics Simulation*, Cambridge University Press, Cambridge, UK, 2004. - D. K. Bhattacharya and G. C. Lie, “Nonequilibrium gas flow in the transition regime: a molecular-dynamics study,”
*Physical Review A*, vol. 43, no. 2, pp. 761–767, 1991. View at Publisher · View at Google Scholar · View at Scopus - D. Frenkel and B. Smit,
*Understanding Molecular Simulation: From Algorithms to Applications*, Academic Press, New York, NY, USA, 2002. - M. Han and J. S. Lee, “Method for calculating the heat and momentum fluxes of inhomogeneous fluids,”
*Physical Review E*, vol. 70, no. 6, Article ID 061205, 2004. View at Publisher · View at Google Scholar · View at Scopus - Y. Zheng, A. L. Garcia, and B. J. Alder, “Comparison of kinetic theory and hydrodynamics for Poiseuille flow,”
*Journal of Statistical Physics*, vol. 109, no. 3-4, pp. 495–505, 2002. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus - K. Xu, “Super-Burnett solutions for Poiseuille flow,”
*Physics of Fluids*, vol. 15, no. 7, pp. 2077–2080, 2003. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus