Franz J. Vesely > CompPhys Tutorial > Selected Applications > Hydrodynamics

 < >
 Part III: Ch. 8

8.4 Direct Simulation Monte Carlo / Bird method

Originally developed for dilute gas flow in engineering and in space science [Bird 94], [Nanbu 83].
Recent extensions for dense gases: [Alexander 95].

Basic idea of the DSMC method for a dilute gas of hard spheres:
1. Setup:
• Divide the sample into cells of volume $V_{c}$, each with $N_{c} \approx 20-40$ particles with given positions and velocities.
• The side length of the cells should be smaller than but of the order of the mean free path.
• Boundary conditions appropriate to the problem are defined, usually specular (reflecting wall) and/or periodic boundaries.
• A time step $\Delta t$ smaller than the typical intercollision time is assumed.

2. Translate all particles according to $\vec{r}_{i} \rightarrow \vec{r}_{i}+\vec{v}_{i}\, \Delta t$, applying the given boundary conditions.

3. Within each cell, draw $M_{c}$ pairs of particles $(i,j)$ that are candidates for a collision:

1. Let the probability of a pair $(i,j)$ to collide depend only on their relative speed $v_{ij} \equiv \left| \vec{v}_{j}-\vec{v}_{i}\right|$ and not to their positions. The argument for this is that all particles in one cell are within free path range of each other. The probability for the pair $(i,j)$ to collide is thus simply proportional to the relative speed: $p_{c}(i,j) \propto v_{ij}$.

Recalling the rejection method of Section 3.2 we draw pairs $(i,j)$ in accordance with this probability density: assuming the maximum of $v_{ij}$ for all pairs in the cell to be known, draw a random number $\xi$ from a uniform distribution in $[0,1]$ and compare it to $v_{ij}/v_{max}$.

Calculating $v_{max}$ would amount to the expensive scanning of all pairs of particles in the cell. Therefore, use an estimated value of $v_{max}'$. If that value is larger than the actual $v_{max}$, the density $p_{ij}$ is still sampled correctly but with a slightly lower efficiency.

2. The total number of collision pairs to be sampled in a cell during one time step is determined as follows. For a gas of hard spheres with diameter $d$ the average number of pair collisions within the cell is $M_{c} \equiv Z \, V_{c} \, \Delta t = \frac{\textstyle \rho^{2} \pi d^{2} \langle v_{rel}\rangle}{\textstyle 2} \, V_{c} \, \Delta t$ where $Z$ is the kinetic collision rate per unit volume, $\rho = N_{c}/V_{c}$ is the number density, and $\langle v_{rel}\rangle$ is the average relative speed. In order to have $M_{c}$ trial pairs survive the rejection procedure of step (a) we have to sample

$M_{trial} \equiv M_{c}\, \frac{\textstyle v_{max}'}{\textstyle \langle v_{rel}\rangle}= \frac{\textstyle \rho^{2} \pi d^{2} v_{max} '}{\textstyle 2} \, V_{c} \, \Delta t$

collisions.

4. Now perform the collision for the pair $(i,j)$. Since the post-collision velocities are determined by the impact parameter which is unknown, they must be sampled in a physically consistent way. In the hard sphere case this is most easily done by assuming an isotropic distribution of the relative velocity $\vec{v}_{ij}$ after the collision. Since the relative speed $\left| \vec{v}_{ij}\right|$ remains unchanged, the problem is reduced to sampling a uniformly distributed unit vector. Marsaglia's recipe may be used for this (see Figure 3.8).