Molecular-level study of compressible gaseous continua

For over 150 years, compressible gaseous continua have usually been formulated in terms of the Navier–Stokes-Fourier (NSF) equations for the macroscopic velocity, pressure and energy as functions of position and time. This study introduces molecular-level investigations of compressible gaseous continua using direct simulation in the Monte Carlo method and molecular dynamics. It is observed that the thermodynamic flux can be generated as long as molecular collisions occur, even without temperature drive and stress tensor work. This is contrary to the general opinion that temperature drive and stress tensor work are the only proximate causes driving thermodynamic flow.

The conventional approach to analysing continuum physics is to use a coarse-grained model in terms of macroscopic variables [1]. For example, fluid dynamics is usually formulated in terms of the Navier-Stoke-Fourier (NSF) equations for the macroscopic velocity, pressure and energy as functions of position x and time t. In addition, the NSF equations have been used for more than 150 years to study the fluid mechanics of gaseous continua. However, doubts persist regarding the completeness of the NSF equations for gas continua [2][3][4][5][6]. In particular, the existence of smooth solutions to the Navier-Stokes equations is still an open problem in physics. In fact, this problem is listed as one of the Millennium Prize Problems in mathematics [7].
Molecular-level study of compressible gaseous continua has received little attention because it requires unacceptably large memory demands in continuum states. However, it has been of great interest to researchers because it can provide fundamental knowledge on flow physics [8].
Our focus in this paper is on the molecular-level study of compressible gaseous continua using direct simulation in the Monte Carlo (DSMC) [9] method and molecular dynamics (MD) [10]. The simulations were performed under conditions sufficient for a small Knudsen number (Kn<0.001) to ensure that the gas behaved as a true continuum. We chose DSMC and MD as simulation methodologies owing to their deterministic characteristics, thereby allowing for realistic molecular behaviour such as intermolecular attractions, repulsions, movements and scatterings, including realistic interactions with solid surfaces [5].

Direct simulation of Monte Carlo
In DSMC, gas is represented by the velocity components and positions of many simulated molecules. The position coordinates, velocity components, and internal states of the simulated molecule are stored. These parameters are updated over time because the molecules are concurrently followed by representative collisions and boundary interactions in the simulated domain [8].
The primary approximation of DSMC is that molecular motions and intermolecular collisions are decoupled over short time intervals. The DSMC algorithm consists of four primary processes: moving the particles, indexing and cross referencing the particles, simulating collisions, and sampling the flow field [11]. The hard sphere (HS) model and diffuse reflection with full accommodation to the surface temperature is implemented in our simulation. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
For validation of the DSMC algorithm in continuum gas flow, we consider the numerical simulations for Ma K n 10, 1.45 10 4 = =´of Edney's type IV shock interaction. The results of the grid sensitivity investigation were presented in our previous work [11]. The number of cells/subcells used for our simulation was 97 060/958 490, proven in [11] to be finest. The computational time step was 1.2 ns. The simulation consisted of 1.56 10 5 time steps to ensure steady-state conditions. We compared the pressure coefficient and Stanton number with the experiments.
where p is the flow-field pressure, p ¥ is the far-field pressure of the flow field, r is the flow-field density, V ¥ is the far-field velocity of the flow field, and q is the thermal flux of the solid surface. The DSMC data for the pressure coefficient and Stanton number are shown in figure 1. The agreement between the DSMC calculation and experiment is generally good. The error between HS-DSMC and experimental data is controlled within 5%.

Molecular dynamics
In molecular dynamics simulations, the time evolution of a set of interacting particles is followed by the solution of Newton's equations of motion denotes the position vector of the ith particle, F i is the force acting upon the ith particle at time t and m i represents the mass of the particle. In addition, 'particles' usually correspond to atoms. To integrate this second order differential equation, the instantaneous forces acting on the particles and their initial positions and velocities need to be specified. Owing to the many-body nature of this problem, the equations of motion are discretized and then solved numerically. The position and velocity vectors define the MD trajectories describing the time evolution of the system. Therefore, they are propagated with a finite time interval using numerical integrators, and in the present study, the Verlet algorithm [12] is used. The position of each particle in space is defined by t r i ( ) and changing in time, whereas the velocities t v i ( ) determine the kinetic energy and temperature of the system. As the particles 'move' at the computer code, their trajectories may be displayed and analysed, providing averaged properties.
The basic formula of the Verlet algorithm can be derived from the Taylor expansions for the positions t r :  the system distributed in accord with a certain statistical ensemble. Our MD simulations were performed using the open source MD simulation package openFOAM [13], which has been validated in many cases.

Results
For Poiseuille flow [14,15] and the gas flow within a circular cylinder [16][17][18], viscous effects dominate the system, and the gas does not compress much; therefore, we consider them to test the applicability of NSF.

Poiseuille gas flow
As illustrated in figure 2, the force-driven compressible Poiseuille flow is defined as a stationary argon flow in a channel in the present DSMC and MD simulations. In addition, it is generated by the action of a uniform external force parallel to the walls in the x direction. The flow between two infinitely separated reservoirs is maintained at the same pressure and temperature by adjusting the density. The force parameter is defined as Here, a denotes the external force that drives the gas flow, h represents the channel height, R is the gas constant for monatomic argon gas, and T w is the wall temperature, which is set to 300 K in this study. The Knudsen number is defined as in two dimensions, Here, n and d represent the number area density and diameter of the molecule, respectively. The computational grid in DSMC consists of 1000 cells along the flow direction and 60 50 cells normal to the flow direction. The base time step for the simulation was 10 −9 s, although this was reduced in very highdensity regions. The solutions were obtained using approximately 6.4 million particles and averaged over a sampling interval of 0.5 ms, at which time the solution had attained steady state. In particular, the body force was constant.
The present MD system is composed of 2 10 6 molecules in a channel of length l d 1000 = and height h d 10 000 . = In addition, we also performed a simulation with a smaller system of 2 10 4 molecules in a channel of length l d 100 = and height h d 1000 .

=
In both cases, the number area density of molecules is n d 0.2 . well within the range of laminar continuum flow. The flow is generated by an external differential force F applied uniformly on each atom in the length-direction. The heat flux due to fluid friction is dissipated through the walls by maintaining the walls at constant temperature (300 K) using a Berendsen thermostat [19,20]. The MD calculations were performed for 280×10 6 steps with an increment of 1 femtosecond. In addition, we performed the integrations by the velocity-Verlet algorithm [21]. The average results across the channel were obtained by dividing the channel into many divided bins in the post-processing of the MD simulation data. The binned values were then time-averaged over many time steps. In each bin, the properties were averaged over all the atoms that caused a representative value for that bin. Spatial averaging was performed across the channel using the binning method [22] with 160 bins. It was observed that the simulation reaches the steady state after approximately 200×10 6 steps. Thus, the properties were averaged for the last 80×10 6 steps after the system had reached the steady state. The density in each bin was subsequently obtained by summing the masses of all atoms. Finally, the pressure and stress tensor and heat flux vector were calculated using the Irving-Kirkwood procedure. The calculation of the pressure and non-conserved variables (in particular, the normal viscous stress and tangential heat flux) directly from the molecular data using the statistical method was extremely important because it is not possible to validate the constitutive laws without information of the viscous stress and heat flux.
After checking the grid convergence, 100 30 30´meshes are implemented in NSF simulations with nonslip and slip wall conditions at constant temperature (300 K). The convergence error of the wall force is set at 10 . 7 - Figure 3 shows the tangential heat flux measured from the values at the centre. Here, the properties are nondimensionalized with respect to the values at the centreline, Surprisingly, the non-classical results of DSMC and MD can describe the occurrence of non-zero tangential heat flux, whereas classical NSF theory fails to do this in all three cases.

Gas flow within a circular cylinder
The second test considers steady gas rotation in a cylindrical container at an angular velocity i . )denotes an orthonormal right-handed set of unit vectors along each of the three coordinate axes. In contrast to the Brenner study [11], this test focused on measuring the radial heat flux distribution in a monatomic gaseous continuum. The temperature is predicted by the NSF equations to be uniform, and the heat flux is zero throughout the rotating gas, independent of radial position (although the pressure and hence density will vary with position due to the action of the centrifugal forces). The radius of the cylinder is 2000d.
The computational grid in DSMC consists of 50 cells along the circumference by 50 cells in the radial direction. The base time step for the simulation was 10 −8 s, although this was reduced in very high-density regions. The solutions were obtained using approximately 25 000 particles and averaged over a sampling interval of 0.2 ms. At this time, the solution had attained steady state.
A domain of 2500 000 gas molecules is simulated in MD, insuring that the continuum assumption is satisfied. Each simulation was solved in parallel on an HPC facility. The calculations were performed for 370×10 6 steps in increments of 1 femtosecond. The integrations were also performed using the velocity-Verlet algorithm. As with the study of Poiseuille gas flow, average results were also obtained by dividing the system into 40 radial bins to measure macroscopic field properties. Each radial bin has equal volume, so the width of each bin is different. The simulation reached the steady state after 300×10 6 steps. Therefore, the properties were averaged over the last 70×10 6 steps, after the system reached a steady state.
After checking the grid convergence, 40 40 meshes are implemented in NSF simulations with nonslip and slip wall conditions at constant temperature. The convergence error of the wall force is set at 10 .

-
The gas was assumed to be the ideal-gas law. Thus, the Mach number is defined by Here, R 0 represents the radius of the cylinder and T 0 is the initial temperature of the wall and the gas. With the gas at rest prior to beginning the system's rotation, the pressure and temperature in the cylinder were set at the values of p atm 1 0 = and T K 300 . 0 = Figure 4 shows several DSMC-and MD-simulated radial heat flux distributions in the cylinder at each of the four Mach numbers indicated. All preceding data show that the heat flow is clearly nonzero within the cylinder, increasing from the centreline to the wall, but the predication of NSF is zero heat flow.

Discussions
We introduced molecular-level investigations of compressible gaseous continua by DSMC and MD. The results of DSMC and MD can describe the occurrence of non-zero thermodynamic flux, whereas classical NSF fails to do this in each case.
Though it is commonly assumed that the NSF equations govern continuum gas flow, this belief, in fact, lacks supporting data. NSF assumes that the temperature-drive and stress-tensor-work are the proximate causes driving thermodynamic flow. This study shows that the thermodynamic flux can be generated as long as molecular collisions occur even without temperature drive and stress tensor work. This is contrary to the general opinion of NSF for over 150 years.
As analyzed by Xu [23,24], all hydrodynamic equations such as NSF are based on the fluid element assumption. Due to the individual isolated element approximation, there is only pressure, viscous friction and heat exchange between the elements, but no mass or molecules penetration, such as the absence of mass diffusion term. In other words, the particle from free penetration between elements is prevented by the intensive particle collision in the framework of NSF. The energy exchange only occurs through the boundary by the heat diffusion in the Fourier's law. In summary, the assumption of the isolated fluid element is a lack of changing neighboring effect, which may have defects in modeling and capturing specific flow behaviors. This is why nonzero thermodynamic flux is found in the solutions of DSMC and MD but not NSF in the present study.
In addition, the results show that solutions of the Navier-Stokes-Fourier (NSF) equations are in conflict with DSMC and MD for compressible gaseous continua. This indicates that the NSF equations do not fully describe transport phenomena occurring in compressible gaseous continua. The present study indicated that both the particle free penetration and fluid element should be considered in the development of new governed equations for gas flows.