Direct simulation Monte Carlo (DSMC) method uses probabilistic Monte Carlo simulation to solve the Boltzmann equation for finite Knudsen number fluid flows.

The DSMC method was proposed by Graeme Bird,^[1][2][3] emeritus professor of aeronautics, University of Sydney. DSMC is a numerical method for modeling rarefied gas flows, in which the mean free path of a molecule is of the same order (or greater) than a representative physical length scale (i.e. the Knudsen number Kn is greater than 1). In supersonic and hypersonic flows rarefaction is characterized by Tsien's parameter, which is equivalent to the product of Knudsen number and Mach number (KnM) or M${\displaystyle ^{2))$/Re, where Re is the Reynolds number.^[4][5] In these rarefied flows, the Navier-Stokes equations can be inaccurate. The DSMC method has been extended to model continuum flows (Kn < 1) and the results can be compared with Navier Stokes solutions.

The DSMC method models fluid flows using probabilistic simulation molecules to solve the Boltzmann equation. Molecules are moved through a simulation of physical space in a realistic manner that is directly coupled to physical time such that unsteady flow characteristics can be modeled. Intermolecular collisions and molecule-surface collisions are calculated using probabilistic, phenomenological models. Common molecular models include the hard sphere model, the variable hard sphere (VHS) model, and the variable soft sphere (VSS) model. Various collision models are presented in.^[6]

Currently, the DSMC method has been applied to the solution of flows ranging from estimation of the Space Shuttle re-entry aerodynamics to the modeling of microelectromechanical systems (MEMS).

DSMC Algorithm

The direct simulation Monte Carlo algorithm is like molecular dynamics in that the state of the system is given by the positions and velocities of the particles, ${\displaystyle \{\mathbf {r} _{i},{\textbf {v))_{i}\))$, for ${\displaystyle i=1,\ldots ,N}$. Unlike molecular dynamics, each particle in a DSMC simulation represents ${\displaystyle F_{N))$ molecules in the physical system that have roughly at the same position and velocity. This allows DSMC to rescale length and time for the modeling of macroscopic systems (e.g., atmospheric entry). Specifically, the system volume is ${\displaystyle V=(NF_{N})/n}$, where ${\displaystyle n}$ is the number density and each collision between simulation particles represents ${\displaystyle F_{N))$ collisions among molecules in the physical system. As a rule of thumb there should be 20 or more particles per cubic mean free path for accurate results.^{[citation needed]}

The evolution of the system is integrated in time steps, ${\displaystyle \tau }$, which are typically on the order of the mean collision time for a particle. At each time step all the particles are moved and then a random set of pairs collide. In the absence of external fields (e.g., gravity) the particles move ballistically as ${\displaystyle \mathbf {r} _{i}(t+\tau )=\mathbf {r} _{i}(t)+\mathbf {v} _{i}(t)\tau }$. Any particle that reaches a boundary or a surface has its position and velocity reset accordingly (e.g., periodic boundary conditions). After all the particles have moved, they are sorted into cells and some are randomly selected to collide. based on probabilities and collision rates obtained from the kinetic theory of gases. After the velocities of all colliding particles have been reset, statistical sampling is performed and then the process is repeated for the next time step.

Collisions

On each timestep the particles are sorted into spatial cells and only particles in the same cell are allowed to collide. Typically the dimension of a cell is no larger than a mean free path. All pairs of particles in a cell are candidate collision partners, regardless of their actual trajectories.

The details of how collisions are calculated in DSMC depend on the molecular interaction model; here we take the hard spheres model, which is the simplest. In the hard spheres model, the collision probability for the pair of particles, ${\displaystyle i}$ and ${\displaystyle j}$, is proportional to their relative speed,

$${\displaystyle P_{\mathrm {coll} }[i,j]=((|\mathbf {v} _{i}-\mathbf {v} _{j}|} \over {\sum _{m=1}^{N_{\mathrm {c} ))\sum _{n=1}^{m-1}|\mathbf {v} _{m}-\mathbf {v} _{n}|))}$$

${\displaystyle N_{\mathrm {c} ))$

1. A pair of candidate particles, ${\displaystyle i}$ and ${\displaystyle j}$, is chosen at random and their relative speed, ${\displaystyle v_{\mathrm {r} }=|\mathbf {v} _{i}-\mathbf {v} _{j}|}$, is computed.
2. The pair is accepted as collision partners if ${\displaystyle v_{\mathrm {r} }>v_{\mathrm {r} }^{\mathrm {max} }\Re }$, where ${\displaystyle v_{\mathrm {r} }^{\mathrm {max} ))$ is the maximum relative speed in the cell and ${\displaystyle \Re }$ is a uniform deviate in [0, 1).
3. If the pair is accepted, the collision is processed; the velocities of the particles are reset but positions are unchanged.
4. After the collision is processed or if the pair is rejected, return to step 1.

This procedure is correct even if the value of ${\displaystyle v_{\mathrm {r} }^{\mathrm {max} ))$ is overestimated, although it is less efficient in the sense that more candidates are rejected.

After the collision pair is chosen, their post-collision velocities, ${\displaystyle \mathbf {v} _{i}^{*))$ and ${\displaystyle \mathbf {v} _{j}^{*))$, are evaluated. Writing the relative velocity in terms of spherical angles, ${\displaystyle \theta }$ and ${\displaystyle \phi }$

$${\displaystyle \mathbf {v} _{\mathrm {r} }^{*}=v_{\mathrm {r} }[(\sin \theta \cos \phi ){\hat {\mathbf {x} ))+(\sin \theta \sin \phi ){\hat {\mathbf {y} ))+\cos \theta \,{\hat {\mathbf {z} ))]}$$

${\displaystyle 2\pi }$

${\displaystyle \phi =2\pi \Re _{1))$

${\displaystyle \Re _{1))$

$${\displaystyle P_{\theta }(\theta )\,d\theta ={\textstyle {\frac {1}{2))}\sin \theta \,d\theta }$$

${\displaystyle q=\sin \theta }$

${\displaystyle P_{q}(q)\,dq=({\textstyle {\frac {1}{2))})\,dq}$

$${\displaystyle \cos \theta =q~\mathrm {and} ~\sin \theta ={\sqrt {1-q^{2))}~\mathrm {where} ~q=2\Re _{2}-1}$$

$${\displaystyle \mathbf {v} _{i}^{*}=\mathbf {v} _{\mathrm {cm} }^{*}+{1 \over 2}\mathbf {v} _{\mathrm {r} }^{*}\qquad \mathbf {v} _{j}^{*}=\mathbf {v} _{\mathrm {cm} }^{*}-{1 \over 2}\mathbf {v} _{\mathrm {r} }^{*))$$

$${\displaystyle \mathbf {v} _{\mathrm {cm} }={1 \over 2}(\mathbf {v} _{i}+\mathbf {v} _{j})={1 \over 2}(\mathbf {v} _{i}^{*}+\mathbf {v} _{j}^{*})=\mathbf {v} _{\mathrm {cm} }^{*))$$

$${\displaystyle v_{\mathrm {r} }=|\mathbf {v} _{i}-\mathbf {v} _{j}|=|\mathbf {v} _{i}^{*}-\mathbf {v} _{j}^{*}|=v_{\mathrm {r} }^{*))$$

From the collision frequency, ${\displaystyle f_{\mathrm {coll} ))$, given by kinetic theory the total number of hard sphere collisions in a cell during a time ${\displaystyle \tau }$ is

$${\displaystyle M_{\mathrm {coll} }={1 \over 2}(N_{\mathrm {c} }-1)F_{N}f_{\mathrm {coll} }\tau =((N_{\mathrm {c} }(N_{\mathrm {c} }-1)F_{N}\pi d^{2}\langle v_{\mathrm {r} }\rangle \tau } \over {2V_{\mathrm {c} ))))$$

${\displaystyle d}$

${\displaystyle V_{\mathrm {c} ))$

$${\displaystyle ((M_{\mathrm {coll} )) \over {M_{\mathrm {cand} ))}=((\langle v_{\mathrm {r} }\rangle } \over {v_{\mathrm {r} }^{\max ))))$$

${\displaystyle \tau }$

$${\displaystyle M_{\mathrm {cand} }=((N_{\mathrm {c} }(N_{\mathrm {c} }-1)F_{N}\pi d^{2}v_{\mathrm {r} }^{\max }\tau } \over {2V_{\mathrm {c} ))))$$

${\displaystyle v_{\mathrm {r} }^{\max ))$