This article first appeared on the website of Micropropulsion and Nanotechnology Laboratory (MpNL), my research group at the George Washington University. It’s reprinted here for convenience.
Particle-In-Cell (PIC) is a technique commonly used to simulate motion of charged particles, or plasma. Here I show you the basics of the method, as applicable to studies of phenomena such as solar wind propagation, or analysis of electric thruster plumes. These discharges are characterized by low plasma densities. At low density, plasma behaves more like a collection of discrete particles, than a single continuous fluid. High-density plasmas are simulated using the extension of computational fluid dynamics into electromagnetics, magnetohydrodynamics.
Modeling of plasmas is complicated by the presence of external and self-induced electromagnetic fields, inter-particle interactions, presence of solid objects, and the different characteristic time scales at which ions and electrons propagate. In order to speed up the calculation, simplifying assumptions are generally taken to match the problem at hand. Here, we are assuming that the current generated by the plasma medium is low enough so that the self-induced magnetic field can be ignored. This is a valid assumption for the class of problems we are dealing with here. As such, the set of underlying Maxwell’s equations is reduced and we obtain an Electro-Static PIC code. In addition, we assume that electrons follow the Boltzmann relationship. The simulation then consists of only the heavy particles, ions and neutrals. This simplification has a tremendous impact on the computational speed, since the time integration can be performed on the much larger ion time scale. Finally, we assume that gas densities are sufficiently low so that particle collisions are negligible.
So what is plasma? It is a collection of charged particles (ions and electrons) as well as non-ionized atoms that is globally charge neutral. In other word, assume that you could grab a chunk of plasma. As long as the radius of this volume is larger than a certain characteristic length (called the Debye lenght), the collection will contain approximately equal amount of positive (from ions) and negative (from electrons and negative ions) charge. The ratio of ionized particles to neutral atoms is the ionization fraction. In
electric thrusters, about 90% of the working gas is ionized.
Figure 1. Representation of plasma as a collection of ions (+), electrons (-) and neutral atoms (o). Charged particles interact with each other through self-induced electric fields.
Charged particles interact with each other by attracting particles of opposite charge and repelling those with the same charge. This is the Coulomb force and is given by: . Conceptually, we could simulate plasma by taking a collection of particles representing the real physical ions and electrons and directly computing this force. However, plasma simulations generally require at least 1 million particles in order to reduce numerical errors. Since Coulomb force leads to an n2 problem, computation of a single time step would require at least 1 trillion operations. Quite a daunting task, especially considering that a typical simulation requires several thousand of such steps!
Introducing the Particle-In-Cell (PIC) method. This method also uses computational particles (called macro-particles) to represent the real ions, electrons and neutrals. However, instead of computing the Coulomb force directly, we obtain the force acting on the particles from the electric field. Force acting on the particles is given by the Lorentz force, , which is a size-n problem. Here q is the particle charge, E is the electric field, and v and B are the particle velocity and magnetic field, respectively. In this example we will assume that plasma is not subjected to any magnetic field (self-induced or applied) and thus the second term in the Lorentz force equation can be ignored. As was eluded to previously, it is not computationally feasible to simulate every physical ion, electron or neutral. As such, each computational particle represents a large number (several tens of million in thruster simulations) of real particles. The ratio of real particles per macro-particle is called the specific weight, sp_wt.
Motion of each macroparticle is governed by Newton’s Second Law,
Here q is the particle charge and m is the particle mass. Electric field is computed from , where is the electric potential. It is given by , the Poisson equation. Here ρ is the charge density and ε0 is the permittivity of free space. In terms of ion and electron number densities, charge density is written as . The subscripts i and e denote ions and electrons, respectively, and Zi is the average ion charge number.
Due to their tiny mass, electrons move at very high velocities (on the order of 5×105 m/s in the types of plasmas considered by MpNL). Resolving electron motion without introducing numerical errors requires the use of extremely small computation time steps (around 10-12 s). With such a small time step, it is not feasible to study any large-scale problems. For instance, it takes couple tens of µs for an ion beam to propagate several thruster lengths from the spacecraft. At the electron time scale, such a simulation would require 10 million time steps. In addition, from the frame of reference of the ions, electrons move instantaneously. Hence it makes sense to simplify the analysis by considering electrons as a fluid. The electron density is then given by the Boltzmann relationship,
The Particle In Cell method consists of an initial setup, the main loop, and a final clean up / results output. All computation happens in the loop. The loop consists of the following steps:
- Compute Charge Density: particle positions are scattered to the grid
- Compute Electric Potential: performed by solving the Poisson equation
- Compute Electric Field: from the gradient of potential
- Move Particles: update velocity and position from Newton’s second law.
- Generate Particles: sample sources to add new particles
- Output: optional, save information on the state of simulation
- Repeat: loop iterates until maximum number of time steps is achieved or until simulation reaches steady state
Charge density is a scalar variable spatially varying in space. It indicates the number of charge units per unit volume. We compute it by distributing charge of all particles onto the nodes of computational cells, and then dividing by the corresponding node volume.
Figure 2. Schematic of the scatter operation. Charge of the particle (in gray) is distributed among the surrounding nodes. Charge contributed to each node is based on the proximity of the particle to that node.
The 1st order (linear) scatter operation is schematically shown in Figure 2. This figure shows a snippet of the computational mesh, and the simulation particle, shown in gray. Charge of the particle is distributed among the nodes of the cell the particle is in. Hence then name for this method, particle in cell. Please note that the total charge carried by each macro-particle is sp_wt*q. The yellow node receives the largest fraction, due to the particle’s close proximity to this node. The smallest amount is contributed to the green node. The weight factors are given by the four area fractions,
where 1, 2, 3 and 4 correspond to the blue, yellow, green and pink nodes, respectively. hx is the fractional distance of the particle from the cell origin in the x direction. Node volumes equal the cell volumes, ΔxΔy. Exceptions exist along domain boundaries, since here only a half (or quarter for corners) of cells is contributing. More complicated boundary conditions require additional treatment. For instance, periodic boundaries require summation of nodes on the plus (+) and the negative (-) side of the domain. This interpolation technique can also be used to scatter data with non-rectangular cells by mapping the physical coordinates to a regular logical domain.
Next, the electric potential is calculated by calling an elliptic equation solver. The Poisson equation can be discretized as
where we used the finite difference with central differencing. This method is modified along mesh boundaries. An elliptic equation describes a boundary value problem and as such boundary conditions must be prescribed along all boundaries. Two types of boundaries exist: Dirichlet and Neumann. The first, Dirichlet, specifies the value of potential. The second boundary specifies the value of the derivative, or electric field. Plasma problems will typically contain a combination of both boundary types. For instance, Figure 3 shows the boundary conditions for the flow of solar wind around a charged plate problem simulated by the sample MATLAB PIC code.
Figure 3. Mix of boundary conditions in the sample PIC code.
Electric field is computed by differencing electric potential, or along boundaries . Only the x direction derivative is shown here, and the one sided boundary equation is applicable to the x- face. The remaining derivatives are obtain in an analogous fashion.
We next integrate particle motion through a time step Δt. The Leapfrog method is commonly used in PIC codes. It is fast and numerically stable. First, we integrate velocity through the time step. Next, position is updated. The name comes from the fact that times at which velocity and positions are known are offset from each other by half a time step. As such, the two quantities leap over each other. This idea is sketched in Figure 4.
Figure 4. The Leap-Frog method.
Numerically, the particle push is given by
The electric field at the position of each particle is computed by a gather operation, which is inversely analogous to the process used to compute charge density.
After particles are moved to new positions, it is necessary to verify that all particles are still in the computational domain. Two boundary interactions are possible. The particles can either exit the domain, or can collide with solid objects. Computational boundaries are either open (or absorbing, allowing particles to leave), reflective (elastically returning particles into the domain) or periodic
(particles are transported to the opposite side of the domain). The reflective boundary is used to identify planes of symmetry. Details of particle-surface interaction, which can result in erosion of native material, is not considered here.
New particles are generated by sampling sources. Sources include spacecraft thrusters or solar wind. This step is problem specific.
In some cases, all particles will be loaded initially, and the simulation will only computes their final distribution state. Typically we are interested in sampling the Maxwellian distribution. As given in Birdsall, this distribution can be approximated as:
where M is some number and Ri is an i-th random number in range [0:1). The value chosen for M controls the accuracy of this method. Birdsall used 12 to prevent entries larger than 6 times the thermal velocity. In the sample MATLAB code, we use M=3.
Finally, since the Leap-Frog method requires velocity and position to be offset from each other, it is necessary to update the velocity of sampled particles by -Δt. This step is frequently omitted. The size of the numerical error introduced by this omission will depend on several factors: the strength of the electric field near the source and size of the time step. This step is also omitted in the sample code, for clarity as it would require duplication of the velocity update code.
Output from PIC codes includes the spatial distribution of plasma parameters such as potential, charge density, electron temperature, as well as particle data such as velocities and current densities. In addition, time-dependant output of global aggregate diagnostic data, such as total kinetic and potential energy of the simulation, is helpful in diagnosing code performance.
The loop repeats until some condition is satisfied. Simulations with continuous sources are run until a steady state is achieved. The steady state is characterized by the net number of particles in the simulation domain remaining constant between time steps. In other words, the number of new particles generated by sources is balanced by the number of particles leaving the domain through the boundaries.
- Birdsall, C.K. and Langdon, A.B.,Plasma Physics via computer simulation, 1985
- Chen, F.F., Introduction to plasma physics and controlled fuison, 1984
- Lieberman, M.A., and Lichtenberg, A.J., Principles of plasma discharges and material processing, 2005
That’s it. Next see the simple Particle In Cell code in Matlab post for implementation. Feel free to leave a comment if you have any questions.Tweet