Loading an isotropic velocity distribution

Posted on May 24th, 2012
Previous Article :: Next Article

Particle methods such as the Particle In Cell Method (PIC) or Direct Simulation Monte Carlo (DSMC) simulate a gas discharge by advancing particles through small steps until some final steady state is reached. But how do we actually get the particles into the simulation in the first place? This article deals with this very important topic, known as particle loading.

The velocity with which particles are born essentially specifies one of the simulation initial conditions. The PIC / DSMC methods simply model the evolution of this initial distribution. Hence, it is crucial that the initial distribution accurately represents the physical situation. Regardless how accurate your solver is, you can’t expect anything but garbage if you feed the solver garbage inputs. Our simulation codes such as Starfish or CTSP allow the user to load particles from both surface and volume sources using a wide range of velocity distribution models. Surface sources are used to model gas injection. These sources inject particles in the direction of the surface normal according to a prescribed distribution. The simplest example would be a uniform beam that simply injects particles with the same velocity and direction. A more realistic model is a drifting Maxwellian. In this source, the particle velocity varies according to the temperature of the gas – the hotter the gas, the greater the variation. The Maxwellian distribution is also used to load background species – those of you with a physics background will right away know that this distribution describes the distribution of velocities in a thermalized gas (one that has undergone a sufficient number of collisions).

Method 1 (Incorrect)

Here we consider the simplest form of a volume source, a source that injects particles with uniform speed and a uniformly-distributed velocity direction. If such a source is a point source, it can be imagined to be a infinitely small light source, with light beams emitting uniformly in all directions. Hence at time 0 we will have only an infinitely small point. At some later time t1 the leading front of the particles will describe the surface of a sphere. The sphere will get bigger as the time goes on, with the density decreasing accordingly.

So how do we inject such particles? First let’s consider the most obvious, but incorrect approach.

/*pick random direction*/
n[0] = -1.0 + 2*Math.random();
n[1] = -1.0 + 2*Math.random();
n[2] = -1.0 + 2*Math.random();
 
/*normalize to a unit vector*/
double n_mag = Math.sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
n[0] /= n_mag;
n[1] /= n_mag;
n[2] /= n_mag;
 
/*assign velocity*/		
part.v[0] = n[0]*v_mag;
part.v[1] = n[1]*v_mag;
part.v[2] = n[2]*v_mag;

This code picks a random direction along each axis, scales the vector to unity, and multiplies it by the velocity magnitude. Simple, but incorrect. If you run this code for a while, you will get a density contour such as the one shown in Figure 1 (the example code for this example is provided in the post on particle data structures). We were expecting a sphere, but instead the particles are expanding as a cube.

cube gas density
Figure 1. The first method gives an incorrect density – we were expecting a sphere but instead got a cube

Method 2 (Isotropic Direction)

The reason why the first method doesn’t work is that it is actually behaving correctly, however, the algorithm is incorrect. Think about what a random number generator does. Its goal is to uniformly fill the space between the lower and the upper limit ([0:1) for the normalized case). In 1D, a random number generator will uniformly fill a line. In 2D, it will uniformly fill a square. And in 3D, it will fill a cube, which is exactly what we saw.

Instead, we need to get a bit more clever. We can divide the problem into two steps. First, let’s consider the x-y plane. We can obtain a uniform circular distribution by first selecting a random angle, and then using trigonometry to obtain the x and y coordinates. In other words, \(\theta=[0:2\pi)\), \(n_x=\cos(\theta)\) and \(n_y=\sin(\theta)\). This loading is shown below in Figure 2.

/*pick a random angle*/
double theta = 2*Math.PI*Math.random();
 
/*compute vector*/
n[0] = Math.cos(theta);
n[1] = Math.sin(theta);
n[2] = 0;
cube gas density
Figure 2. Getting there – a 2D uniform distribution

Almost there, we just need to add the third direction. Let’s think about what we want to accomplish. We want the vector to be uniformly distributed along the z axis. At the same time, we want to assure that the vector has a unit magnitude. This will require scaling the x-y component by some factor a. Hence we write \(n_z = [-1:1]\), \(n_x=a\cos(\theta)\) and \(n_y=a\sin(\theta)\). We need \(n_x^2+n_y^2+n_z^2=1\), or \(a^2\cos^2(\theta)+a^2\sin^2(\theta)+n_z^2=1\). This gives us \(a=\sqrt{1-n_z^2}\) by the trig identity. We modify our example code accordingly,

/*pick a random angle*/
double theta = 2*Math.PI*Math.random();
 
/*pick a random direction for n[2]*/
double R = -1.0+2*Math.random();
double a = Math.sqrt(1-R*R);
 
n[0] = Math.cos(theta)*a;
n[1] = Math.sin(theta)*a;
n[2] = R;

and voila, we get the right distribution, as shown in Figure 3! As a little aside, the implementation here has a little bug. Java’s Math.random returns values in the range of [0:1) (zero inclusive but one exclusive). Hence, there will be a slight preferential minus z direction, since \(n_z=[-1:1)\). This could be corrected by using a modified RNG that is one-inclusive.

isotropic density
Figure 3. The expected spherical distribution

2 comments to “Loading an isotropic velocity distribution”

  1. Amitavo
    May 24, 2012 at 9:01 am

    THIS ARTICLE IS VERY USEFUL FOR ME AS I AM A BEGINER FOR PIC AND THANKS FOR POSTING. I NEED SOME SUPPORT.I WANT TO DISTRIBUTE MY 30 PARTICLE HAVING BEAM RADIUS OF 10.1MM IN R-Z PLANE OF A CAVITY.I AM ALSO USING BORIS SCHEME FOR PARTICLE ADVANCEMENT.I AM FACING PROBLEM HOW I CAN DISTRIBUTE THEM IN A CORRECT WAY.IS THERE ANYBODY WHO CAN GUIDE ME OR SUPPORT ME IN DOING THIS??PLEASE HELP ME.

  2. May 24, 2012 at 9:41 am

    Thanks Amitavo. What exactly are you trying to accomplish? Uniform density or uniform spacing between particles? If it’s the latter, then that’s quite easy, simply load the particles at offsets of 10.1mm/29. However, if you want uniform density then you need to take into account the radial position. In RZ, the cell volume is given by something like dz*dr*r*dtheta. In that case what you want to do is calculate the average mass density (total mass of the 30 particles divided by the injection volume) and then compute the mass to inject at each cell by multiplying by the cell volume. You can also get the volume from the source surface area multiplied by the injection velocity.

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre lang="" line="" escaped="" cssfile=""> In addition, you can use \( ...\) to include equations.