# 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.

/***** START OF DSMC CLASS *****/
var dsmcCell = (function() {
var AMU = 1.66053886e-27;	/*atomic mass*/
var K = 1.3806488e-23;		/*Boltzmann constant*/
var V_c=0.001*0.01*0.001;	/*cell volume*/
var sigma_cr_max=1e-18;		/*initial guess*/
var Fn = 5e8;				/*macroparticle weight*/

var delta_t=2e-5;			/*time step length*/
var it;						/*current iteration*/
var vdf;					/*public array to hold VDF data*/

var part_list = new Array(); /*particle list*/

/*particle class, creates particle with isotropic velocity,
see: http://www.particleincell.com/2012/isotropic-velocity/ */
function Part(vdrift,temp,mass)
{
this.v=[];

/*sample speed*/
var vth = sampleVth(temp,mass);

/*pick a random direction*/
var theta = 2*Math.PI*Math.random();
var R = -1.0+2*Math.random();
var a = Math.sqrt(1-R*R);

this.v[0] = vth*Math.cos(theta)*a;
this.v[1] = vth*Math.sin(theta)*a;
this.v[2] = vth*R + vdrift;

this.mass=mass;
}

/*creates initial particle populations*/
function init()
{
/*clear data, important if restarting*/
part_list = [];
it = 0;

/*first create slow particles*/
for (var i=0;i<1000;i++)
part_list.push(new Part(100,5,32*AMU));

/*and now the fast ones*/
for (var i=0;i<1000;i++)
part_list.push(new Part(300,5,32*AMU));

/*show initial plot*/
vdf = computeVDF(0,450,25);
plotXY(vdf);
}

/* performs a single iteration of DSMC*/
function performDSMC()
{
/*fetch number of particles*/
var np = part_list.length;

/*compute number of groups according to NTC*/
var ng_f = 0.5*np*np*Fn*sigma_cr_max*delta_t/V_c;
var ng = Math.floor(ng_f+0.5);	/*number of groups, round*/
var nc = 0;			/*number of collisions*/

/*for updating sigma_cr_maxx*/
var sigma_cr_max_temp = 0;
var cr_vec = [];

/*iterate over groups*/
for (var i=0;i<ng;i++)
{
var part1,part2;

/*select first particle*/
part1=part_list[Math.floor(Math.random()*np)];

/*select the second one, making sure the two particles are unique*/
do {part2=part_list[Math.floor(Math.random()*np)];}
while (part1==part2);

/*relative velocity*/
for (var j=0;j<3;j++)
cr_vec[j] = part1.v[j]-part2.v[j];
var cr = mag(cr_vec);

/*eval cross section*/
var sigma = evalSigma(cr);

/*eval sigma_cr*/
var sigma_cr=sigma*cr;

/*update sigma_cr_max*/
if (sigma_cr>sigma_cr_max_temp)
sigma_cr_max_temp=sigma_cr;

/*eval prob*/
var P=sigma_cr/sigma_cr_max;

/*did the collision occur?*/
if (P>Math.random())
{
nc++;
collide(part1,part2);
}
}

/*update sigma_cr_max if we had collisions*/
if (nc)
sigma_cr_max = sigma_cr_max_temp;

/*show diagnostics*/
vdf = computeVDF(0,450,25);
plotXY(vdf);

it++;

/*output more data to the console*/
console.log("it: "+it+"\tng: "+ng+"\tnc: "+nc+"\tscmax: "+sigma_cr_max);

/*run for 500 time steps*/
if (it<500)
setTimeout(performDSMC,10);
}

/*samples from 1D Maxwellian using the Birdsall method*/
function maxw1D()
{
return 0.5*(Math.random()+Math.random()+Math.random()-1.5);
}

/*returns random thermal velocity*/
function sampleVth(temp,mass)
{
var sum=0;
var vth = Math.sqrt(2*K*temp/mass);
for (var i=0;i<3;i++)
{
var  m1 = vth*maxw1D();
sum += m1*m1;
}

return Math.sqrt(sum);
}

/*evaluates cross-section using a simple (non-physical) inverse relationship*/
function evalSigma(rel_g)
{
return 1e-18*Math.pow(rel_g,-0.5);
}

/*returns magnitude of a 3 component vector*/
function mag(v)
{
return Math.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
}

/*performs momentum transfer collision between two particles*/
function collide (part1, part2)
{
/*center of mass velocity*/
var cm = [];
for (var i=0;i<3;i++)
cm[i] = (part1.mass*part1.v[i] + part2.mass*part2.v[i])/
(part1.mass+part2.mass);

/*relative velocity, magnitude remains constant through the collision*/
var cr=[];
for (var i=0;i<3;i++)
cr[i] = part1.v[i]-part2.v[i];

var cr_mag = mag(cr);

/*pick two random angles, per Bird's VHS method*/
var theta = Math.acos(2*Math.random()-1);
var phi = 2*Math.PI*Math.random();

/*perform rotation*/
cr[0] = cr_mag*Math.cos(theta);
cr[1] = cr_mag*Math.sin(theta)*Math.cos(phi);
cr[2] = cr_mag*Math.sin(theta)*Math.sin(phi);

/*post collision velocities*/
for (var i=0;i<3;i++)
{
part1.v[i] = cm[i]+part2.mass/(part1.mass+part2.mass)*cr[i];
part2.v[i] = cm[i]-part1.mass/(part1.mass+part2.mass)*cr[i];
}

}

/* bins velocities into nbins between min and max */
function computeVDF(min, max, nbins)
{
/*pointers for easier access*/
vel = [];
bin = [];

/*set delta range*/
var delta = (max-min)/(nbins-1);

/*shift left by half delta to center bins*/
min -= 0.5*delta;
max -= 0.5*delta;

/*set initial values*/
for (var i=0;i<nbins;i++)
{
vel[i] = i/(nbins-1);
bin[i] = 0;
}

/*compute histogram*/
for (var i=0;i<part_list.length;i++)
{
var part = part_list[i];
var vmag = mag(part.v);
var b = Math.floor((vmag-min)/(delta));
if (b<0 || b>=nbins-1) continue;
bin[b]++;
}

/*normalize values*/
var max_bin=0;
for (var b=0;b<nbins;b++)
if (bin[b]>max_bin) max_bin=bin[b];

for (var b=0;b<nbins;b++)
bin[b]/=max_bin;

return [vel, bin];
}

/*expose public members*/
return {init:init,performDSMC:performDSMC,vdf:vdf};
})();

/**** END OF DSMC CLASS *****/

/*grab canvas and related properties*/
var canvas=document.getElementById("cell");
var ctx = canvas.getContext("2d");
var height=canvas.height;
var width=canvas.width;

/*paint background before we scale the context*/
ctx.fillRect(0,0,width,height);

/*create axis*/
ctx.strokeStyle="black";
ctx.fillStyle="purple";
ctx.lineWidth=3;
setLineDashSafe([10,6]);
ctx.strokeRect(0.05*width,0.05*height,0.9*width,0.9*height);
setLineDashSafe([]);

/*save background*/
var saved_bk = ctx.getImageData(1, 1, width, height);

/*scale to [0,1][0,1] and put origin at left/bottom to simplify graphing*/
ctx.translate(0.05*width,0.95*height);
ctx.scale(0.9*width,-0.9*height);

/*show initial distribution*/
dsmcCell.init();

/*handler for the button click*/
function startSim()
{
dsmcCell.init();
dsmcCell.performDSMC();
}

/*generates XY plot from xy[2][]*/
function plotXY(xy)
{
var x=xy[0],y=xy[1];

/*restore background*/
ctx.putImageData(saved_bk, 1, 1);

/*start plotting*/
ctx.beginPath();
ctx.moveTo(0,0);
for (var i=0;i<x.length;i++)
ctx.lineTo(x[i],y[i]);

ctx.lineTo(1,0);
ctx.lineTo(0,0);
ctx.fill();
}

/*helper function since Firefox does not seem to support setLineDash*/
function setLineDashSafe(style)
{
if (ctx.setLineDash)
ctx.setLineDash(style);
}