Direct Simulation Monte Carlo (DSMC) Method

Introduction

In this article, we’ll introduce Direct Simulation Monte Carlo (DSMC), which is a method for simulating interparticle collisions. It is frequently used to simulate neutral gas flows and to add collisions to plasma simulations. This method has been popularized by G.A. Bird, and to this day, his book Molecular Gas Dynamics and the Direct Simulation of Gas Flows serves as the bible for the DSMC method. Prof. Bird also runs a website containing various DSMC resources. Please make sure to check out these references for the details of DSMC, the write up here is definitely not exhaustive!

DSMC is the analogue to the particle in cell method for neutral gases. It is a method for computing the motion of neutral flows directly using simulation particles instead of solving fluid conservation equations, as is done in the CFD approach. The main advantage over CFD is that DSMC does not assume a priori any distribution function. The CFD method is generally applicable only to thermalized flows – ones in which gas molecule velocities follow the Maxwellian distribution function. The DSMC method works by moving the simulation particles through a finite time step through a discretized computational domain. After each push, particles are paired up with other particles located in the same computational mesh and collision probability is computed for these pairs. Those pairs that are colliding then have their velocities modified in a way that conserves energy and momentum.

Since by itself, the PIC method does not account for particle collisions, the particle in cell method must be coupled with a collisional algorithm to account for momentum transfer and charge exchange collisions. Previously we have introduced one way of adding collisions: the Monte Carlo Collisions (MCC) method. The problem with MCC is that it is not momentum conserving, and is thus physically suitable only when collisions have only a negligible effect on the target population. This significantly restricts what we can do with it. But no worries, DSMC is here to the rescue!

Quick DSMC Theory Crash Course

Let’s consider a single simulation cell containing a number of simulation macroparticles. Bird gives the probability of a collision for a macroparticle in this cell as
$$P=F_N\sigma_Tc_r\Delta t/V_c$$
where \( F_N\) is the macroparticle weight, \(\sigma_T\) is the total collision cross-section, \(c_r\) is the relative velocity, \(\Delta t\) is the simulation time step, and \(V_c\) is the cell volume. You can think of this expression as
$$P=F_N\left(\dfrac{\sigma_Tc_r\Delta t}{V_c}\right)$$
or the ratio of the volume swept by the cross-section to the total volume of the cell scaled by the macroparticle weight.

No Time Counter Method (NTC)

One order to check for collisions would be to iterate over all \(N\) particles, and compute probability with all remaining particles. This would result in \(N(N-1)/2 \sim (N^2)/2\) pairs. For a large number of particles, this method clearly becomes very computationally inefficient. Bird’s No Time Counter (NTC) method was designed to help with this issue. It allows us to estimate ahead of time the maximum number of pairs that need to be checked. The maximum collision probability is
$$P_{max}=F_N\left(\sigma_T c_r\right)_{max} \Delta t / V_c$$
where \(\left(\sigma_T c_r\right)_{max}\) is a parameter chosen ahead of time using approximate predictions of cross-section and velocity. As Bird points out in his book, the actual value is not all that important since it ends up getting cancelled out. The number of pairs to check is then
$$(1/2)N\bar{N} F_N\left(\sigma_T c_r\right)_{max} \Delta t / V_c$$
Here the second \(\bar{N}\) is the average particle count. The average is used to reduce the statistical time-step to time-step oscillations. For each pair, we then compute the probability as follows
\(P = \dfrac{\sigma_T c_r}{\left(\sigma_T c_r\right)_{max}}\)
this value is compared to a random number and the collision occurs if \(P>R\).

Variable Hard Sphere

The collision cross-section for a hard-sphere, billiard ball collision is \(\sigma=\pi d^2\). This model is not physically realistic for real molecules. Bird proposed a new model which takes into account the decrease in cross-section as relative velocity between molecules increases. The diameter of the sphere is allowed to vary according to
$$d=d_{ref}(c_{r,ref}/c_r)^\nu$$

All that’s needed now is the code for actually performing the collision. Bird goes into the details of deriving the equations in Chapter 2 of his book. In short, momentum is conserved in the collision such that
\(m_1\mathbf{c}_1 + m_2\mathbf{c}_2 = m_1\mathbf{c}^*_1+ m_2\mathbf{c}^*_2 = (m_1+m_2)\mathbf{c}_m\)
Bird also shows that energy conservation dictates that the magnitude of relative velocity \(\mathbf{c}_r =\mathbf{c}_1-\mathbf{c}_2\) remains constant during the collision. We can also write the post-collision velocity as
$$ \mathbf{c}_1^* = \mathbf{c}_m +\dfrac{m_2}{m_1+m_2}\mathbf{c}_r^*$$
and
$$ \mathbf{c}_2^* = \mathbf{c}_m -\dfrac{m_1}{m_1+m_2}\mathbf{c}_r^*$$
Since we know the magnitude of the relative velocity from the pre-collision velocities, we just need to select a new direction for the post-collision velocity vector. In general, the post-collision scattering is anisotropic – the direction will depend on the pre-impact velocities and differential cross-sections are then used to select the in-plane angle. However, here for simplicity we assume that the scattering is isotropic. We then let
$$\cos \theta = 2R_1-1$$
and
$$\phi = 2\pi R_2$$
where the two R1 and R2 are two random numbers. The relative velocity is then given by
$$\mathbf{c}_r^*=<\!\cos \theta,\sin\theta\cos\phi,\sin\theta\sin\phi\!>|\mathbf{c}_r|$$

Zero D Example

We now have all the theory needed to implement a simple DSMC example. The example is implemented using Javascript and HTML5 canvas element. Simply click the “Start the Simulation” button below the canvas to run 500 iterations of the algorithm from above. This example illustrates thermalization of two gas populations. The plot shows the histogram of the velocity distribution function. Initially, it shows two distinct populations, but after a while, due to collisions, the two populations merge into a single one.

Your browser does not support the <canvas> element. Please use a different browser or click here if you are viewing this in an email.

Continue reading…

This article continues in Part 2: details of the Javascript DSMC code