One of the governing equations for electrostatic plasma simulations is the Poisson’s equation,
In the types of discharges we typically consider here at PIC-C, the ion density is obtained from kinetic particles (the particle-in-cell method) while electrons are approximated as a fluid. The electron thermal velocity assuming 5eV electrons is over one million meters per second. On the other hand, typical ion velocities in plumes are around 15 km/s. This is a speed ratio of three orders of magnitude. The PIC method requires that particles travel less than a cell length per time step. As such, the simulation time step will be determined by the fastest species. Hence, if want to include electrons, we need to run the simulation at 1000 times smaller time step than what would be satisfactory for the ions. In other words, the simulation will take 1000 times longer to complete the same physical time span.
Luckily, we can often get away with not simulating electrons as kinetic particles. This is true especially when magnetic fields are not present or we are dealing with motion in direction along the magnetic field lines. Consider the electron momentum equation,
In the frame of reference of ions, electrons are seen to respond instantaneously to a disturbance, and thus the time dependent term can be neglected. The convective term can also be generally neglected since its magnitude is small (or zero if we consider spatially uniform velocity). This then leaves us with the force balance, . We can next make three substitutions, , , and (assuming isothermal electrons), which gives us . This equation can be easily integrated along some arbitrary path to obtain
. This very important relationship is known as the Boltzmann relationship. It is encountered frequently in plasma physics. The importance of the Boltzmann relationship lies in the fact that it relates local potential and the local electron density. Given a set of reference parameters, electron potential can be computed everywhere in the domain (where the governing assumptions hold) from simply knowing the electron density.
Nonlinear Poisson’s Equations
The Boltzmann relationship can also be inverted to obtain an expression for electron density,
This is the missing piece needed for our electrostatic potential solver! Instead of simulating electrons as particles, we can use this relationship directly in the Poisson’s equation,
Unfortunately, this relationship results in a non-linear equation, since now the expression we are solving for (potential) appears on the right hand side as well. In matrix form, we have . Luckily, solving non-linear equations is no rocket science – we just need little help from Mr. Newton.
Newton’s Method for Systems
Newton’s method is a simple algorithm for finding function roots. Consider some arbitrary equation . We can use Taylor’s expansion to compute the value of this function in the vicinity of a known point,
. For small differences the second order terms drop out and we can write
Numerically, the above equation is used to get consecutively closer approximations to the function root, x_new = x_old - f(x_old)/f_prime(x_old). We can use a similar approach to treat a system of equations. The matrix analogue to the single equation is , where
is the Jacobian. This algorithm is generally implemented in two steps two avoid the need to compute a matrix inverse. First, the system is solved for using your favorite linear solver. The solution is then updated as . The process continues until , some prescribed tolerance.
The non-linear Poisson solver algorithm
So how does this method look in practice when applied to the Poisson’s equation? Scroll down for the entire Matlab source code (note, the coding was actually done in Octave since we don’t have access to a Matlab license). To keep things simple, we apply the algorithm to a 1D system bounded on both ends with fixed potential nodes. The potential is -5V on the left boundary and +5V on the right boundary. Uniform plasma density of and temperature is used.
Using the standard finite difference algorithm, we can write our governing equations as
Here the constant . In other words, we have the system (we are using for generality). The two vectors on the right hand side correspond to the constant and the non-linear term. As the first step in our algorithm, we compute the non-linear vector (the linear doesn’t change between iterations).
%1) compute bx for n=1:nn if (fixed_node(n)) bx(n)=0; else bx(n) = -C*den0*exp((x(n)-phi0)/kTe); end end
The example is made more general by using a fixed_node array to indicate which nodes are fixed and which are fluid. We compute electron density only on the fluid nodes. Next we add the two right hand side vectors and subtract to obtain .
%2) b=b0+bx b = b0 + bx; %3) F=Ax-b F = A*x-b;
We next need to compute the Jacobian. Our function consists of three terms: the linear system given by the matrix multiplication, a vector which is not a function of x, and another vector which is. The derivative of the non-linear RHS vector is . The Jacobian is thus given by . We first compute this P vector on fluid nodes,
%4) compute P=dbx(i)/dx(j), used in computing J for n=1:nn if (fixed_node(n)) P(n)=0; else P(n) = -C*den0*exp((x(n)-phi0)/kTe)/kTe; end end
and set the Jacobian,
%5) J=df(i)/dx(j), J=A-diag(dbx(i)/dx(j)) J = A - diag(P);
The rest is easy. We first solve for y and then use it to update the solution. The final profile is plotted in Figure 1.
%6) Solve Jy=F y = J\F; %7) update x x = x - y;
The iterations continue until convergence
%8) compute norm; l2 = norm(y); if (l2<1e-6) disp(sprintf("Converged in %d iterations with norm %g\n",it,l2)); break; end
Complete Source Listing
Here is the complete source code in Matlab. The code was actually tested using Octave, so let us know if it doesn’t work in Matlab proper. Also, feel free to contact us if you need a Java implementation.
%Demo non-linear Poisson solver % we are solving: % d^2phi/(dh^2) = -rho/eps0*(ni-n0*exp((phi-phi0)/kTe)) % with phi=-5, 5 on left and right boundary % See http://www.particleincell.com/2012/nonlinear-poisson-solver/ % for additional info %clear screen clc %constants EPS0 = 8.85418782e-12; Q = 1.60217646e-19; %setup coefficients den0=1e16; phi0=0; kTe=5; %precomputed values lambda_D = sqrt(EPS0*kTe/(den0*Q)); %for kTe in eV dh = lambda_D; %cell spacing C = -Q/EPS0*dh*dh; %setup matrixes nn=15; %number of nodes A=zeros(nn,nn); fixed_node = zeros(nn,1); b0=zeros(nn,1); %left boundary A(1,1)=1; b0(1)=-5; fixed_node(1)=1; %right boundary A(nn,nn)=1; b0(nn)=5; fixed_node(nn)=1; %internal nodes for n=2:nn-1 %FD stencil A(n,n-1)=1; A(n,n)=-2; A(n,n+1)=1; fixed_node(n)=false; b0(n)=C*den0; end %initial values bx = zeros(nn,1); P = zeros(nn,1); J = zeros(nn,nn); x = zeros(nn,1); y = zeros(nn,1); %--- Newton Solver ---- for it=1:20 %1) compute bx for n=1:nn if (fixed_node(n)) bx(n)=0; else bx(n) = -C*den0*exp((x(n)-phi0)/kTe); end end %2) b=b0+bx b = b0 + bx; %3) F=Ax-b F = A*x-b; %4) compute P=dbx(i)/dx(j), used in computing J for n=1:nn if (fixed_node(n)) P(n)=0; else P(n) = -C*den0*exp((x(n)-phi0)/kTe)/kTe; end end %5) J=df(i)/dx(j), J=A-diag(dbx(i)/dx(j)) J = A - diag(P); %6) Solve Jy=F y = J\F; %7) update x x = x - y; %8) compute norm; l2 = norm(y); if (l2<1e-6) disp(sprintf("Converged in %d iterations with norm %g\n",it,l2)); break; end end disp(x'); plot(x);