Simple Particle In Cell Code in Matlab
This post includes the sample Particle-In-Cell (PIC) code that goes with our previous article on the electrostatic particle-in-cell method. It simulates the flow of uniform plasma past an obstruction – a charged plate in our case. This simple code is implemented in Matlab. If you don’t have Matlab, you can run the code using Octave, which is an open-source Matlab emulator. If you do end up getting Octave, as of the date of this article, the MinGW build available on Octave-Forge is out of date and contains a buggy contouring function. In order to use the latest version – and see the plots – you need to install Cygwin and add Octave as a Cygwin package…
The sample PIC code is contained in two files: flow_around_plate.m and eval_2dpot_GS.m. The first file is the actual PIC code, while the second file holds the potential solver. The solver computes electric potential using the Finite Difference method. The code simulates the flow of solar wind (low density plasma) over a charged plate. The code is not particularly fast (it’s in Matlab after all!), but shouldn’t take more than few minutes to run on a standard PC. If you run the code through Octave, make sure to first execute more off to avoid pagination of the output.
The simulation plasma environment is set on lines 24-28:
n0 = 1e12; %density in #/m^3 phi0 = 0; %reference potential Te = 1; %electron temperature in eV v_drift = 7000; %ion injection velocity, 7km/s phi_p = -7; %plate potential
By changing the plate potential, you can alter the size of the wake that forms behind the plate. At sufficiently large values, ions are attracted back to the plate. This behavior is analogous to the plasma environment that forms around satellites utilizing electric propulsion. Slow moving ions are extracted out of the thruster plumes and can impact surfaces without a direct line of sight to the thruster.
The simulation domain size and the number of time steps is set starting with line 35:
nx = 16; %number of nodes in x direction ny = 10; %number of nodes in y direction ts = 200; %number of time steps
At the default settings, 200 steps is sufficient to reach steady state. The simulation time step is such that a particle traveling at the drift velocity will cross 10% of the cell length per time step. Since particles speed up as they approach the plate, this conservative time step keeps particles from traveling more than a cell per time step. This is a necessary constraint in all PIC codes, but in general you want to keep the particles traveling a shorter distance (I generally use 25% in my simulations). The number of particles generated at each time step is specified on line 39,
np_insert = (ny-1)*15; %insert 15 particles per cell
Plate dimension are set on lines 48 and 49:
box(1,:) = [floor(nx/3) floor(nx/3)+2]; %x range box(2,:) = [1 floor(ny/2)]; %y range
For an explanation of the physics used in the code, please see the Particle In Cell page. Feel free to leave a comment if you have any questions.
You may also be interested in a Finite Element Particle In Cell example.
Leapfrog Particle Push (Velocity Integration)
Loading an isotropic velocity distribution
Direct Simulation Monte Carlo (DSMC) Method