Current Density Limit

Posted on May 27th, 2014
Previous Article :: Next Article

One difference between ion thrusters and Hall Effect thrusters is that ion thrusters are space charge limited devices. The region between the grids of the ion thruster is occupied solely by ions, and is thus subject to the Child Langmuir law. This law, given by
$$ J_{CL}=\frac{4}{9}\epsilon_0 \sqrt{\frac{2q}{m}} \frac{|\Delta \phi|^{1.5}}{(\Delta x)^2} $$
gives the maximum current density that can be extracted between two planar electrodes $latex Delta x$ apart and having potential difference of $latex Delta phi$. To illustrate the concept, I put together a small Matlab program for my propulsion class at the George Washington University. The program (current.m) is given below. It simulates flow of singly-charged Xenon ions, emitted from an anode at x=0, and moving towards a cathode located 1cm away, with potential difference of 100V. For this configuration, the Child Langmuir law predicts maximum current density of 0.0477 A/m2.

% Simulation to study Child Langmuir current density limit
% For more see:
more off;			%turn off pagination in Octave
%physical constants
AMU = 1.6605e-27;	%atomic mass in kg
E = 1.6022e-19;		%elementary charge in Coulomb
EPS0 = 8.8542e-12;	%permittivity of free space F/m
ME = 9.1094e-31;	%electron mass
%set boundaries
phi_left = 0;
phi_right = -100;
x_left = 0;
x_right = 0.01;
%compute mesh parameters
nx = 100;				%number of mesh nodes
delta_x = (x_right-x_left)/(nx-1);	%cell spacing, there are nx-1 cells
%set matrix data
A = zeros(nx);
A(1,1) = 1;		%fixed potential at left
A(nx,nx) = 1;		%fixed potential at right
for i=2:nx-1
	%finite difference stencil for second derivative
	A(i,i-1) = 1;
	A(i,i) = -2;
	A(i,i+1) = 1;
%ion data
ion_mass = 131.3*AMU;
ion_charge = 1*E;
%time data
delta_t=0.5*delta_x/sqrt(2*ion_charge/ion_mass*(abs(phi_right-phi_left))) %time step
it_max = 3000;		%number of time steps
%injection current
j_CL = (4/9)*EPS0*sqrt(2*ion_charge/ion_mass)*(abs(phi_right-phi_left))^1.5/(x_right-x_left)^2;
I = 1*j_CL		%this is also current density due to unit length in y and z
np_inject = 10;	%number of particles to inject per time step
ion_spwt = I*delta_t/(ion_charge*np_inject);	%number of actual particles per macroparticle
part_count = 0;		%number of particles
efield = zeros(nx,1);
%main loop
for it=1:it_max
	%*** 1) generate new particles ***    
	for p=1:np_inject		%inject particles
		%increment particle counter
		part_x(part_count) = 0;	%initial position
		part_v(part_count) = 0;	%initial velocity
		part_new(part_count) = 1;	%mark new particles
		%compute electric field at particle position
		fi = (part_x(part_count)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		part_ef = efield(i)*(1-di)+efield(i+1)*(di);	
		%rewind the velocity back by 0.5*delta_t, this part of the leapfrog scheme
		part_v(part_count) = part_v(part_count) - (ion_charge/ion_mass)*part_ef*0.5*delta_t; 
	%*** 2) compute field quantities***
	u = zeros(nx,1);
	count = zeros(nx,1);		%this field holds number of macroparticles
	for p=1:part_count
		%get particle mesh index
		fi = (part_x(p)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		count(i) = count(i)+(1-di);
		count(i+1) = count(i+1)+ (di);
		u(i) = u(i)+(1-di)*part_v(p);
		u(i+1) = u(i+1)+ (di)*part_v(p);
	den = ion_spwt*count/delta_x;	%divide by "volume"
	den(1) = den(1)* 2.0;			%only half area at boundaries
	den(nx) = den(nx)*2.0;
	%average velocity, make sure we don't divide by zero
	u(count~=0) = u(count~=0)./count(count~=0);
	u(count==0) = 0;
	j = den.*u*ion_charge;	%current density
	%*** 3) set charge density
	rho = den*ion_charge;		
	%*** 4) solve potential ***
	b = -rho/EPS0*delta_x*delta_x;
	b(1) = phi_left;
	b(nx) = phi_right;
	phi = Ab;
	%*** 5) update electric field ***
	efield = -gradient(phi,delta_x);
	%*** 6) move particles ***
	for p=1:part_count
		%obtain electric field at the particle position
		fi = (part_x(p)-x_left)/delta_x + 1;
		i = floor(fi);
		di = fi-i;
		part_ef = efield(i)*(1-di)+efield(i+1)*(di);	
		% set random delta_t for newly born particles
		if (part_new(p))
			dt = (1.0+1.0*rand())*delta_t;
			part_new(p) = 0;	%clear new particle flag			
			dt = delta_t;
		%update particle velocity from F=m(dv/delta_t)=qE
		part_v(p) = part_v(p)+(ion_charge/ion_mass)*part_ef*dt;
		%update particle position
		part_x(p) = part_x(p)+part_v(p)*dt;
		%remove particles that crossed boundaries
		if (part_x(p)<x_left || part_x(p)>=x_right)
			%remove particle by replacing the current array data with the last
			part_x(p) = part_x(part_count);
			part_v(p) = part_v(part_count);
			part_new(p) = part_new(part_count);
			part_count=part_count-1;		%decrement particle counter
			p=p-1;		%also decrement current index to rescan this position
	%*** 7) output some diagnostics
	%plot results
	if (mod(it,25)==0 || it==it_max)
	disp(sprintf('it = %dt j_ave = %4.3g',it,mean(j(j~=0))));
		x_pos = [1:nx]*delta_x+x_left;
 		title('Electric Field');
 		title('Current Density');
 		axis([x_left,x_right, 0, 1.2*I]);
		axis([x_left,x_right, -2000, 12000]);
		%for saving animation

Code Testing

First, we can make some simple sanity checks, by considering what happens if there is no electric field. This was done by commenting out the line to calculate electric field:

%*** 5) update electric field ***
%	efield = -gradient(phi,delta_x);

and changing the injection velocity to a non-zero value:

	part_v(part_count) = 1000;	%initial velocity

Figure 1 shows results for this case. As expected, the velocity remains constant and so does the injection current.

noef_j noef_v
Figure 1. Results with no electric field


Next, we can consider an actual simulation with the electric field enabled. For the first case, we set the injection current density to $latex 0.2 j_{CL}$ by using the following line

I = 0.2*j_CL		%this is also current density due to unit length in y and z

As shown in Figure 2, for this case the current density remains constant at the prescribed value of 0.0095 A/m2. The electric field and potential are also plotted.

01phi 01ef 01j 01v
Figure 2. Results with injection current density at 0.2 JCL

The results get more interesting once the injection current density exceeds the Child Langmuir law. Figure 3 shows these results for injection current of 2 J_CL. Once the discharge becomes space charge saturated, a negative electric field forms at the anode, which then prevents additional ions to enter the simulation. Only as many ions are allowed to pass as is needed to maintain the current density at the Child Langmuir limit.

2j 2v
Figure 3. Results with current density set to 2.0 JCL. Once the discharge becomes space charge saturated, a negative electric field forms at the anode to prevent injection of additional ions.

We can demonstrate this behavior further by giving ions a small initial velocity. This will cause them to travel a short distance from the injector electrode before being turned back. We can also visualize the dynamic behavior by creating an animation. To do this, I used the process described here to export png files for individual frames and then used ImageMagick to create the animated gif with the following commands:

$ find -name 'v*.png' -maxdepth 1 -exec convert {} -resize 400x400 small/{} ;
$ cd small
$ convert -delay 10 -loop 0 v*.png v_anim_small.gif

The figure below shows the average velocity. Negative value indicates that the majority of ions are moving in the opposite direction, towards the emitter electrode.

animation of current density limit
Figure 4. Animation of the dynamic behavior for a space charge saturated discharge.
Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.

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.