This note summarizes the kinetic theory for the Direct Simulation Monte Carlo (DSMC) method. The main contents are from the classic book of I.D. Boyd, Nonequilibrium Gas Dynamics and Molecular Simulation.
The model input for three typical methods:
- Interatomic potential energy surface (PES) $\to$ molecular dynamics (MD) calculations.
- The collision cross section $\to$ direct simulation Monte Carlo (DSMC).
- Transport property coefficients $\to$ computational fluid dynamics (CFD) calculations.
Our goal here is establish a relation between microscopic molecular and macroscopic gas dynamics. The procedure is as follows:
- molecular description of a nonequilibirum gas $\to$ Boltzmann equation
- Take moments of Boltzmann equations $\to$ a set of conservation equations
- In the limit of near-equilibrium flow $\to$ conservation equations reduce to the same form as Navier-Stokes equations.
- Compare the two sets of the equations $\to$ macroscopic state properties and transport properties.
1 Chapman-Enskog analysis
It determined the precise velocity distribution function that reduces the Boltzmann equation to the Navier-Stokes equations. Two important features:
- rigorously verifies the Newtonian and Fourier transport models that relate the shear stress tensor and that flux vector to linear functions of velocity and temperature gradients, respectively.
- establishes the quantitative limits on the accuracy of these transport models and therefore on the accuracy of the Navier-Stokes equations for nonequilibrium flows.
1.1 Analysis for BGK
By performing Chapman–Enskog analysis, there are no longer unknown coefficients. Rather, the transport terms and transport coefficients are fully determined according to the collision operator in the Boltzmann equation. however, the Prandtl number predicted by the BGK model is unity, which is a well-known deficiency of the BGK collision operator.
1.2 Analysis for Boltzmann equation
When the full Boltzmann collision integral is included the coefficients become integral expressions with no closed-form solution. The mathematical procedure used to derive these coefficient expressions can be found in Chapman and Cowling (1952) and is also presented in Chapter 10 of Vincenti and Kruger (1967).
The Prandtl number now is 2/3 (for monotonic atoms), which is physically accurate.
Transport coefficients can be expressed as a function of cross section. [see Eq.(5.103 ~ 106)]. The viscosity cross section and momentum cross section are similar to the total collision cross section.
2 Collision cross section
2.1 Definition
Collision cross section is the integral result of the scattering angle.
The scattering angle is a function of the impact parameter, relative speed, and interatomic potential.
2.2 Comparison for the MD and DSMC
The potential energy surface (PES) that governs interatomic forces is what determines the scattering angle of individual collisions. PES is the model input for molecular dynamics calculations. The collision cross sections (Eqs. 5.103 and 5.104) are simply the integrated result of all impact parameters that lead to a finite scattering angle. Since in a dilute gas impact parameters and the pre-collision orientations of molecules are completely random, it is not necessary to model every possible collision arrangement and only the integrated result (i.e., the cross section) is required. For this reason, the collision cross section is a physically meaningful and convenient quantity in many molecular models of dilute gases. It appears in the collision integral in the Boltzmann equation and serves as the main model input for the DSMC method.
3 Hard-Sphere Interactions
The Hard-Sphere model is employed to find the function form of scattering angle.
In summary, using a hard-sphere interatomic potential results in a collision cross section that is independent of relative velocity (g), and therefore a collision integral and viscosity coefficient inversely proportional to the hard-sphere diameter squared (d^2) and proportional to sqrt(T). This temperature dependence is not accurate for real gases and is an artifact of the hard-sphere assumption.
4 DSMC method
The DSMC method simulates the Boltzmann equation, including all relevant physics, and thus makes none of the assumptions and simplifications discussed above, such Hard-Sphere assumption.
The DSMC method leverages three physical characteristics of a dilute gas:
- Molecules move in free flight without interaction for time scales on the order of the local mean collision time.
- The impact parameters and initial orientations of colliding molecules are random.
- There are an enormous number of molecules per cubic mean free path and only a small fraction need be simulated to obtain an accurate molecular description of the flow.
The corresponding three technique/assumptions of DSMC:
- simulated molecules can be moved in straight lines for a fraction of their mean collision time without any loss in accuracy.
- the impact parameters of colliding molecules in a dilute gas are random and do not need to be deterministically simulated as they are with MD.
- highly precise distribution functions, and thus a full molecular description of the flow, can be obtained by modeling only a small fraction of the real number of molecules.
Each technique is related to the other two, and the combination of these three techniques enables the DSMC method to accurately simulate macroscopic nonequilibrium flow fields.
A key concept of the DSMC method is that a simulated particle does not represent a distribution of real molecules. Rather, a simulated particle represents a large number of identical real molecules.
A DSMC simulation cannot determine precisely which real molecules actually collide and what impact parameters characterize the initial conditions of each collision. This loss of determinism is also a result of moving simulated particles in straight lines for a fraction of their mean collision time. In effect, although the distribution functions are accurately resolved at spatial and temporal scales of the mean free path and mean collision time, the precise locations of the real molecules comprising those distributions are no longer known below these scales.
In this manner, collisions between simulated particles within each DSMC cell are performed stochastically. During each time step, simulated particles are randomly paired up within each cell, and subsequently tested for a collision. Thus all possible collision pairs, consisting of all classes of molecular velocities, internal energies, and chemical species, are sampled in a Monte Carlo fashion from the actual distributions of particles within each cell.
4.1 The cell size and time step required for DSMC
Since simulated particles are randomly paired for collisions within each cell, DSMC collision cells must be sized at or below the local mean free path. If cells were significantly larger than the mean free path, molecules separated by large distances could be randomly selected to collide. Such collisions would transfer mass, momentum, and energy over unphysically large distances. The same error would be produced if particles were moved in straight lines for time steps much larger than the mean collision time.
The resulting errors in the computed flow field would be analogous to the numerical dissipation error that results from using a coarse mesh or too-large time step for a CFD calculation. Essentially, sharp flow gradients would not be accurately resolved and heat and momentum transfer to surfaces would also be incorrect.
4.2 Central DSMC algorithm
-
Generate particles at inflow boundaries.
- Move all particles in straight lines along their molecular velocity vectors for a timestep less than the local mean collision time.
- Apply boundary conditions to particles that collide with a surface.
- Remove particles that exit the simulation domain.
-
Perform collisions stochastically within each cell, assuming impact parameters for collision pairs to be random.
-
Sample particle properties in each cell.
- Return to (1).
Molecular properties can be sampled in each cell to compute velocity and internal energy distribution functions, as well as average quantities (macroscopic quantities) such as density, bulk velocity, temperature, and pressure. If there are many particles in each cell then such average quantities and even distribution functions can be resolved with high precision at each time step of simulation. This may be desired for the simulation of unsteady flows using DSMC.
In summary, the DSMC method takes advantage of the inherent properties of dilute gases by using simulated particles that each represent a large number of identical real molecules, moving particles with time steps on the order of the mean collision time, and stochastically selecting collision pairs and initial orientations within volumes (computational cells) on the order of the mean free path. These are rigorous simplifications based on sound physical principles.
4.3 Particle Movement and Sorting
After all particles are moved and sorted locally into cells, all remaining physical models operate only on the list of particles contained within each cell. Thus, the rest of the algorithms are to be applied to the particles within a single cell, independently of any other cell in the computation.
4.4 Collision Rate
The number of collisions that occur within a DSMC cell during a single time step (the simulated collision rate) is directly linked to the local transport properties of the gas.
In the limit of equilibrium flow, a DSMC simulation should reproduce the equilibrium collision rate of a real gas and therefore reproduce the real gas viscosity.
Assume there are $N$ molecules in the volume $V$, they formed $N(N-1)/2$ possible molecule pairs. Since each DSMC particle represents a large number ($W_p$) of identical real molecules, then the following probability should be applied all $N_p(N_p-1)/2$ particle pairs within the cell:
\[P_{DSMC} = (\pi d^2 W_p)g\Delta t_{DSMC}/V_{DSMC}\]In this manner, the total number of collisions to be performed within a DSMC cell during a single time step for a hard-sphere gas is
\(N_{coll} = \frac{1}{2}N_p(N_p -1)\frac{\pi d^2W_p<g>\Delta t_{DSMC}}{V_{DSMC}}\)
where $
For a simple gas, in the limit as $W_p \to 1$ (and thus $N_p-1 \to N-1 \approx N$, for large N), this is identical to the analytical expression for the hard-sphere collision rate in an equilibrium gas (given by kinetic theory).
Actually, in G.A. Bird’s classic DSMC book, Bird first “guess” the expression of $N_{coll}$ and then prove that it can correct reproduce the equilibrium collision rate.
However, it is important to note that since the probability of collision $P_{DSMC}$ is applied to each particle pair individually, the DSMC method correctly predicts different collision rates for different particle-pair classes. For example, for hard-sphere particles, pairs with higher relative velocities collide more frequently compared to pairs with lower relative velocities. Thus, in addition to reproducing the correct equilibrium collision rate, DSMC also predicts the correct nonequilibrium collision rate and is accurate for an arbitrary (non-Maxwellian) velocity distribution function, found within the interior of a shock wave, for example.
- Modification $\to$ No-Time-Counter (NTC)
The cross section of a molecule pair is not constant ($\sigma \neq \pi d^2$), but a generally a strong function of the relative velocity ($\sigma$ $=\sigma(d,g)$). Since the average relative velocity of molecule pairs in a gas is proportional to the gas temperature, this significantly alters the temperature dependence of the local collision rate (and therefore transport properties) compared to the hard-sphere result, and must accounted for.
The hard-sphere procedure is extended to more realistic interactions, by allowing the cross section for a molecule pair to include a dependence on relative velocity. The probability applied to each molecule pair now becomes,
\[P_{DSMC} = \sigma(d,g)gW_p\Delta t_{DSMC}/V_{DSMC}\]and the number of collisions performed within a DSMC cell becomes, \(N_{coll} = \frac{1}{2}N_p(N_p-1) \frac{<\sigma(d,g)g> W_p \Delta t_{DSMC}}{V_{DSMC}}\) where $<\sigma(d,g)g>$ represents the average of all molecule pairs in the cell.
Further considerations
- If $N_p$ is small, the fluctuation of $N_p$ leads to inaccurate result of collision rate. Typical requirement for steady flow is $N_p>10$, while for unsteady flow, it should be much larger.
- The $P_{DSMC}$ applied to $N_p(N_p-1)$ pairs is typically small. Instead of applying a small probability to O(N_p^2) particle pairs, one can obtain an identical result(the required number of collisions, or collision rate) by applying a larger probability to a smaller number of pairs.
A maximum for the number of particle pairs expected to collide is determined as \(N_{coll-max} = \frac{1}{2}N_p(N_p-1)\frac{[<\sigma(d,g)g>]_{max} W_p \Delta t_{DSMC}}{V_{DSMC}}\)
As emphasized by G.A. Bird, if $N_p$ doubled by halving $W_p$, the $N_{coll-max}$ only doubles, which is linear as $N_p$ (not $N_p^2!$).
The round off value of $N_{coll-max}$ is the actual number of particle pairs to be tested pairs, or attempted pairs: \(N_{attempted} = floor(N_{coll-max} + 0.5)\)
Assuming that each particle can collide only once per time step, this number must be limited as \(1 \leq N_{attempted} \leq floor(N_p/2)\)
Note: $N_{coll-max}$ is the number of collision pairs that exactly corresponds to a constant value of $[<\sigma(d,g)g>]{max}$ applied to all particle pairs, wheres $N{attempted}$ is the collision pairs that will actually be tested in DSMC cell. Owing to the discrete nature of simulation (the two value may not equal) , the differences \(F_{correction} = N_{coll-max}/N_{attempted}\) must be corrected at each time step.
The true collision rate is then simulated by randomly selecting $N_{attempted}$ particle pairs, and accepting each pair for an actual collision with probability \(P_{DSMC}= F_{correction} \times \frac{\sigma(d,g)g}{<\sigma(d,g)g>]_{max}}\) By substituting $N_{coll-max}$ into above equation, the $<\sigma(d,g)g>]{max$ disappears. The product of $P{DSMC} \times N_{attempted}$ reduced exactly to the number of collisions that should be performed within a DSMC cell.
4.5 Variable Hard-Sphere model
The VHS model allows the total cross section of an interaction to have a dependence on the relative velocity of the collision pair. Thus, VHS is able to predict the correct temperature dependence of the viscosity coefficient.
The cross section corresponding to the VHS model uses a variable diameter: \(d = d_{ref}(\frac{g_{ref}}{g})^\nu\)
and therefore \(\sigma_T^{VHS} = \pi d_{ref}^2 (\frac{g_{ref}}{g})^{2\nu}\)
The value of the power-law exponent ($\nu$) can be set to match the slope of a real cross section. On substituting this VHS viscosity cross section into the equation for the viscosity coefficient given by the kinetic theory, we have \(\mu^{VHS} = f(d_{ref},g_{ref}) T^{1/2+\nu}\)
where $f$ is a function of reference diameter and relative velocity. We can define a reference viscosity, then the final expression for the viscosity simulated by the VHS is: \(\mu^{VHS} = \mu_{ref}^{VHS}(\frac{T}{T_{ref}})^w\) where \(\mu_{ref}^{VHS} = \frac{15\sqrt{2\pi m_rkT_{ref}}}{2(5-2w)(7-2w)\pi d^2_{ref}}\) and \(w \equiv \nu + 1/2 = 2/\alpha + 1/2,\) with the parameters list in the table:
| w | $\nu$ | $\alpha$ |
|---|---|---|
| viscosity exponent | power-law exponent for VHS | exponent for inverse power-law model |
Therefore, the input species parameters for VHS model are:
- molar mass, $m$;
- reference temperature, $T_{ref}$;
- reference diameter, $d_{ref}$;
- viscosity exponent, $\omega$.
It is important to note, however, that although the correct equilibrium viscosity is simulated, the nonequilibrium collision rate is not based on any underlying data, and therefore is not strictly validated.
4.6 DSMC - how it works?
The majority of DSMC collision models are specific to the collision pair, for example, species $i$ colliding with species $j$.
All the macroscopic flow quantities can be sampled from the particles straightforward.
- Kinetic theory establishes the relationship between transport coefficients (e.g., viscosity) with the cross section.
- VHS model is developed to relate the total cross section $\sigma_T^{VHS}$ to the relative velocity of the collision pair， which can be used to determine the collision probability.
- Experimental data for the viscosity is employed to determine reference diameter.
- The total collisions within a DSMC cell in a time step, or collision rate, is calculated by simple kinetic theory.
- The TCE model is employed to reproduce the required collision rate.
- The particle collisions and move are decoupled in DSMC. Therefore, after the collision, the particles move in straight lines.
- Apply the boundary conditions.
- Sampling procedure, conduct the statistical average to obtain the steady results and fluid properties.
