Simple Particle In Cell Code in Matlab

Introduction

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…

File Organization

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.

Input Parameters

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.

Figure 1. Simulation results, plasma potential and ion density

Download

Download the files here: flow_around_plate.m, eval_2dpot_GS.m

Did you find this article useful? If so, please share the link with your friends, coworkers, and colleagues. And don't forget to subscribe to the newsletter and follow us on Google+ and Twitter. To find out more about how we can assist you on your project, please see the request for proposal page.