Simulation and Optimization of an Inertial Electrostatic Confinement Fusion Reactor
by
Jimmy Zhan
A thesis submitted to the Department of Physics, Engineering Physics and Astronomy in
conformity with the requirements of ENPH 455
Supervisor: Dr. Robert Knobel
Faculty of Applied Science and Engineering
Queen's University
Kingston, Ontario, Canada
April 2, 2013
Abstract
The Inertial Electrostatic Confinement Fusion Reactor (IEC Fusion) is a potential method for controlled fusion power. The IEC fusion reactor confines plasma via a strong electrostatic potential between two concentric grids as ions stream towards the center. The electrostatic potential between the outer grid and the inner grid is on the order of magnitude of tens of kilo volts, and the outer grid is usually grounded. Gas species such as Deuterium and Tritium are injected into a vacuum chamber housing the grids, usually at pressures on the order of a millitorr. Once ionized by the electric field, the ions accelerate radially toward the central grid, colliding with other ions of sufficient relative velocity to undergo fusion. Those ions that shoot through the center without fusing will decelerate and fallback towards the center, making another pass. This method of confinement therefore has two advantages: Plasma recirculates around the central grid making multiple passes and hence increasing the likelihood of fusion. Optimal fusion cross section is on the order of hundred million Kelvins, while difficult to achieve conventionally, it translates to only tens of kilo electron volts, which is fairly easy to achieve via electrostatics. There are also difficulties associated with IEC fusion. Several loss mechanisms limit its ultimate power output, the most significant being Coulomb interactions, causing bremsstrahlung radiation and thermalization of ion velocities, and ion-grid collisions. This thesis investigates the fusion performance of a cylindrical IEC device, this will be done via computational modelling, and iterative numerical simulations of performance parameters. A set of optimum design parameters will be derived from the simulations, and a preliminary design for an IEC device will be given. The results were not surprising, there is an optimum operating voltage of approximately 190 kV. The optimum pressure and chamber size is approximately 1 meter at 1 millitorr. The grid transparency has a dramatic effect on the performance as well, but this is material limited to 98%. A novel method of reducing grid collisions is being researched by the US Navy currently, in which electrons are trapped to form a virtual cathode, using magnetic fields. This is termed the Polywell.
List of Figures
Figure 1 The original Farnsworth fusor.
Figure 3 The correct model for fusion cross section [10].
Figure 4 Coulomb interaction and Bremsstrahlung radiation of an electron and ion [10]
Figure 5 Charge exchange process of various isotopes of Hydrogen [15]
Figure 6 Comparison of fusion performance of various anode-cathode geometries [29].
Figure 7 The simulation reached near steady state in approximately 10 μs, with the number of ions and electrons equal. [29]
Figure 8. Ion-ion reactions occur only inside the inner grid, ion-neutral collisions rate increases approaching the center, and neutral-neutral collisions occur uniformly throughout device. [29]
Figure 9 A total neutron flux of roughly 8.8∙10^5 n/s was achieved, at a grid voltage of 100 kV and current of 40 mA. [29]
Figure 10 The basic schematic of an IEC fusion reactor [16].
Figure 11 The simulation cell, containing one macro-particle [21].
Figure 12 The Leapfrog method, showing the time offset between position and velocity [21].
Figure 13 The dependence of fusion rate as a function of grid voltage. The fuel pressure is kept at 0.001 millitorr, and the chamber radius is 1 m.
Figure 14 The dependence of the fusion rate plotted against the cathode grid transparency. Note that this is a log plot. The voltage is at 60 kV, the fuel pressure is kept at 0.001 millitorr, and the chamber radius is 1 m.
Figure 15 The dependence of fusion rate on the input fuel gas pressure. The voltage is kept at 60 kV, and the chamber radius at 1 m.
Figure 16 The dependence of fusion rate on the chamber size. The voltage is at 60 kV, and the pressure is at 0.001 millitorr.
Figure 17 The Paschen Curve for various gases, illustrating the “Paschen Minimum”. The breakdown voltage of air on either side of the minimum value has inverse relationship with the product of pressure and distance. [32]
Figure 20 Fusion Rate of Optimum Design.
Figure 21 Definition of the hard-sphere cross section model [10].
Figure 22 The Coulomb force repels, the nuclear force attracts, but only acts at a short range [10].
Figure 23 The classical model of fusion cross section [10].
Figure 24 The correct model of fusion cross section, including the effects of tunneling, resonance, and high-speed decay [10].
Figure 25 Experimentally measured cross sections σ(v) for various fusion reactions of interest. Notice that DT has a peak of 5 barns at 120 keV, and is clearly favorable in comparison to DD [10].
Figure 26 <σv> for various fusion reactions of interest, as functions of temperature. Note that the peak cross section for DT is at 70 keV [10].
Definitions of Variables
Boltzmann Constant
Planck Constant
Permittivity of free space
Speed of light
Electron charge
Charge of ion species x
Electron mass
Mass of ion species x
Density of electrons
Density of ion species x
Temperature of electrons
Temperature of ion species x
Fusion cross section
Velocity-averaged fusion cross section
Velocity vector
Position vector
Time
Pressure
Electrode gap distance (≈ radius of IEC chamber)
Potential
Fusion rate density
Energy per fusion reaction
Power density of process x
Acknowledgements
I want to particularly thank my thesis supervisor, Dr. Robert Knobel for providing me with unwavering support and guidance over this past year. This thesis would not have been completed without him.
I would also like to thank Dr. Jordan Morelli for lending me his book “Plasma Physics and Fusion Energy” for one year. The book provided the foundation upon which I constructed my thesis.
I am also grateful for the team at Particle in Cell Consulting, and Tech-X Corporation, for providing support on the utilization of PIC codes for plasma simulations.
Lastly, I would like to thank all my friends here at Queen’s, for this wonderful past four years.
Table of Contents
Inertial Electrostatic Confinement
Mathematical Modelling Methods
Chamber Size and Pressure Limit
Dependence on Cathode Grid Transparency
Simulation of My Proposed Design
Conclusions and Recommendations
Improvements in IEC Simulation
Improvements in IEC Technology
Appendix A Derivation of the Fusion Cross Section
Velocity Averaged Cross Section
Appendix B Additional MATLAB Plots to support Optimization
Introduction
The two conventional means of fusion are the inertial confinement and magnetic confinement fusion. Both of these methods are thermonuclear fusion. Thermonuclear fusion operates by increasing the average energy of the plasma, which has a Maxwellian distribution, so that some of the ions will have enough energy to undergo fusion [10].
In magnetic confinement reactors, plasma is confined via a magnetic field and it is heated to hundreds of millions of Kelvins. In inertial confinement fusion, a high power laser hits and ignites a pellet of compressed fusion fuel to initiate fusion [10]. Since this thesis is not on thermonuclear fusion, further details pertaining to their operations will not be given. The conclusion drawn is that thermonuclear fusion reactors are inefficient (heating up a Maxwellian distribution but only the tip is capable of fusion), large, complex, and compatible with only a few fusion fuels [16].
The Inertial Electrostatic Confinement fusion reactor was originally conceived by television pioneer Philo T. Farnsworth while investigating various vacuum tube designs for televisions. He noticed that electrons moving inside the vacuum tube can be stopped and accumulate at a location by the application of high frequency magnetic fields [16]. This is later known as the multipactor effect. Farnsworth devised the idea of plasma confinement using a vacuum tube. The walls of the reactor consists of electrons or ions held in place by the multipactor, fuel species are injected through the wall. He termed this design a virtual electrode, and the system a fusor [16]. Farnsworth later patented his design in 1968 (Figure 1).
The original Farnsworth fusors used arrays of ion guns and did not consider the possibility of ion recirculation [16]. In the late 1960, Robert Hirsch modified the original Farnsworth fusor, and designed spherical grid-based fusors, which he later patented (Figure 2). This design is termed the Hirsch-Meeks fusor. The machine consists of two concentric spherical grids that are mostly transparent, and held at different potentials. Ions accelerate radially inward, and those that did not collide would decelerate through the other side, stop, and accelerate back towards the center for another pass, thus conserving energy [16].
Figure 1 The original Farnsworth fusor.
Figure 2 Hirsch-Meeks fusor.
Fusion reactors based on the Hirsch-Meeks fusor configuration are simple and theoretically would have higher efficiencies than thermonuclear fusion reactors, because they accelerate the ions monoenergetically, and the ions recirculates through the core reactor region. In reality, IEC reactors are not perfect. The major problem is Coulomb interactions. They cause bremsstrahlung radiations, which are radiations generated from the deceleration of ions or electrons when they collide. Coulomb interactions also cause ion velocity thermalization in a matter of milliseconds [25], the plasma (initially monoenergetic) approach a Maxwell-Boltzmann distribution. The other problem with IEC is the non-perfect transparency of the inner grid [24]. A typical grid can be made no more than 98% transparent and still support its own weight, due to material limitations [16]. Ions recirculating through the core region will bombard the grid, losing energy and causing grid corrosion. Given that the average time for an ion to undergo fusion in a typical IEC device is several minutes, thermalization and grid collisions will contribute significantly to power loss [25].
Objectives
Qualitative Objectives
Research background theories relevant to nuclear fusion, and develop mathematical models to characterize an IEC device, numerically simulate the operation of the device, and iterate various design parameters. Factors such as fusion cross sections, ion density, and power density will be considered. Design parameters such as grid geometry, chamber size, gas pressure, and operating voltage will be optimized. Energy loss mechanisms will also be explored, including bremsstrahlung radiation, thermalization of ion velocities, and grid collisions.
Quantitative Objectives
Using the above results, the original objective was to design an IEC device capable of continuous operation, generating a flux of 10 million neutrons per second. This quantity is not arbitrary: Directly quantifying the power output of a fusion reactor is difficult without a means of extracting and converting that power into electricity. Thus neutron flux is a common measurement of fusion performance. IEC fusion reactors have been built by many amateurs in the recent past, and the best achieved a neutron flux of 1 million per second [16]. For a university level thesis, I decided to increase the performance by a factor of ten.
Revisions
It was noted that the scope of the original objective was too broad and optimistic, and thus it was decided to narrow down the objectives. In particular, emphasis was placed on the modelling and simulation of theoretical IEC fusion reactors rather than the design of a particular device. The simulation results will be used to optimize several design parameters, but no physical design will be given.
Theory
Fusion Basics
Fusion is a nuclear reaction in which two or more nuclei collide at high speed and form a new type of nucleus. During the process, some of the mass of the nuclei are converted into energy, via 
[1]. Fusion of nuclei with masses lower than iron (which has the largest binding energy per nucleon) gives off energy, while fusion of nuclei with masses greater than iron absorbs energy [1].
The fusion rate of a nuclear reactor is given by
Eq. 1 [10]
Cross Sections
The fusion cross section, denoted by σ, characterizes the probability that a pair of nuclei will collide and undergo fusion. I encourage you to read the full treatment and the derivation of the correct model of fusion cross section, in Appendix A. The result is that the cross section must be explained via quantum mechanical effects such as barrier tunneling, resonance effects, and wave-particle duality (Figure 3) [10]. Equally important is that the cross section must be integrated over our Maxwellian distribution in order to get the velocity-averaged cross section that applies to fusion reactors [10].
Figure 3 The correct model for fusion cross section [10].
Inertial Electrostatic Confinement
IEC fusion shares the same “fusion mechanics” as any other fusion, but differs in the method of plasma confinement. The reactor operates by accelerating ions towards a cathode grid, placed concentrically inside an anode grid of the same geometry (which is usually the metal wall of the vacuum chamber), via a potential difference. Various geometries for the grids can be employed, spherical and cylindrical being the most common [29]. Numerous studies over the past four decades demonstrated that IEC fusion yields very high power density and specific power [7].
Advantages of IEC Fusion
The major advantage is that since the ions are accelerated, they will reach the central cathode region monoenergetically, and the resulting plasma is non-Maxwellian. Thus the bulk of the plasma will undergo fusion instead of only the tip of the distribution [16].
Any given ion will likely recirculate many times through the central region before hitting the grid, and thus its probability of undergoing fusion is significantly increased [19].
Another advantage is that its design is inherently simple and light weight, thus making the IEC well suited for space travel technologies [18]. It does not need large arrays of lasers used for inertial confinement, and there is no large magnets and lithium sheaths as in magnetic confinement. A bonus due to lack of magnets is that IEC will not experience cyclotron radiation.
Key Issues
The main issue preventing IEC device from achieving its theoretical performance is Coulomb interactions. Coulomb interactions cause ion thermalization [16], and bremsstrahlung radiations [25]. This means that the IEC reactor will be forced to spend energy maintaining a non-thermal distribution.
Bremsstrahlung (braking) radiation is caused by the Coulomb interactions of charged particles and their subsequent deflections and radiations, emitting photons (mostly X-rays). This is shown schematically in Figure 4.
Figure 4 Coulomb interaction and Bremsstrahlung radiation of an electron and ion [10]
The full derivation of Bremsstrahlung radiation in a fusion reactor is long, the results are shown.
Eq. 2
where 
is the effective charge of an ion, and is equal to 1 for DT reactions [10].
One solution to reduce Coulomb interactions is to increase the potential well depth, and thus increasing the ion energy. This is because Coulomb collision cross sections are smaller at higher energies [10].
Another process adversely affecting the IEC reactor performance is charge exchange (Figure 5). When a moving ion collides with a neutral atom it will capture one of its electrons, thereby becoming neutral. The transfer of kinetic energy will be small (several eV) [6]. However a neutral atom will no longer be accelerated by the electric field, nor will it be trapped by the potential well. The neutral atom will draw additional energy for re-ionization.
Figure 5 Charge exchange process of various isotopes of Hydrogen [15]
The last issue is ion-grid collisions, leading to further thermalizations and grid corrosion. Currently, material science is not advanced enough to create a transparent enough grid that can support its own weight [16].
Plasma
Plasma is a state of matter consisting of ions, neutrals, and electrons. An important plasma parameter is known as the Debye length, as it determines the overall charge of a chunk of plasma. Neutrality of plasma holds for volumes of radius much greater than the Debye length, while on the scale of the Debye length there can be charge imbalance [21]. The Debye length is given as follows,
Eq. 3
Plasma contains charged particles, and charged particles interact with each other by the Coulomb force, given as,
Eq. 4
Previous Experimental Results
The University of Wisconsin-Madison Fusion Technology Institute carried out various IEC experiments and simulations. They designed a spherical geometry IEC reactor named UW-IEC [29]. They compared the fusion performance of spherical and cylindrical geometries, and concluded that spherical devices generally perform better with the same operating conditions, and scales better with increasing grid voltages (Figure 6) [29].
Figure 6 Comparison of fusion performance of various anode-cathode geometries [29].
The team at UW-M numerically simulated their UW-IEC device, using the Particle-In-Cell code (OOPIC) developed at UC-Berkeley. The simulator is 1D in space and 2D in velocity space, and contains Monte Carlo algorithms for electron-neutral and ion-neutral collisions (Figure 7, Figure 8, Figure 9).
Figure 7 The simulation reached near steady state in approximately 10 μs, with the number of ions and electrons equal. [29]
Figure 8. Ion-ion reactions occur only inside the inner grid, ion-neutral collisions rate increases approaching the center, and neutral-neutral collisions occur uniformly throughout device. [29]
Figure 9 A total neutron flux of roughly 8.8∙10^5 n/s was achieved, at a grid voltage of 100 kV and current of 40 mA. [29]
Mathematical Modelling Methods
Modelling of IEC device generally comes down to plasma modelling. This involves solving kinematic equations coupled with either Maxwell’s Equations (for electromagnetic fields), or Poisson’s Equation (for electrostatic fields) [3].
Single Particle Description
This model treats the plasma as individual ions and electrons moving in imposed (rather than self-consistent) electric and magnetic fields. In the case of an IEC device, the field would be a constant electrostatic field. The force on the particle is calculated from Lorentz Force Law, and the particles’ motions follow Newtonian mechanics.
Kinetic Description
This is the most popular plasma model. The kinetic model produces a distribution function, representing the density of plasma in 7d phase space (space, velocity, and time) [3]. The kinetic model can be solved either by the Boltzmann Equation, when collisions are needed, or by the Vlasov Equation, when long range Coulomb interactions are necessary (but forgoes treating collisions) [22]. The kinetic model can also be solved by the Fokker-Planck Equation, which is similar to the Boltzmann Equation with the collision terms simplified. The advantage of the kinetic model is that particles’ distribution function is defined everywhere, and therefore no assumption needs to be made about the type of distribution (ideal for non-Maxwellian plasmas). There are generally two types of kinetic models. One is based on representing a smoothed distribution function on a grid in phase space. The other one keeps track of the kinetic information by following a large number of individual particles (thus is more accurate, albeit at the expense of computational intensity). This second method is often achieved by a numerical method called Particle-In-Cell [3]. This thesis will use the second type of the kinetic model as the basis for IEC plasma simulation.
Fluid Description
The fluid description reduces some of the complexities in the kinetic description. The simplest fluid model treats the plasma as a single fluid, as in magnetohydrodrynamics, which is governed by Maxwell’s Equations coupled with Navier-Stokes Equations [3]. More general fluid models treat the plasma as two fluids (ions and electrons). The fluid model describes the plasma based on macroscopic quantities such as density, mean velocity, and mean energy []. These macroscopic fluid equations are obtained by taking velocity moments of the above kinetic equations (Boltzmann and Vlasov Equations). The fluid model of plasma must be accompanied by transport coefficients such as mobility, diffusion coefficient, and average collision frequencies, which are determined once a velocity distribution is assumed. Thus some inaccuracies will be introduced when assuming a particular distribution. The fluid model is often accurate enough when there are many collisions (or when the plasma density is high), causing the plasma velocity to become Maxwellian [27]. This is not ideal for IEC modelling because theoretically the distribution of an IEC plasma is monoenergetic, and low density fuel scenarios should also be investigated.
Numerical Simulation Methods
The most oft-used numerical method for solving partial differential equations for plasma simulation is the Particle-In-Cell (PIC) method [21]. In this method, individual particle’s kinetic information is tracked in a Lagrangian continuous phase space, and the moments of the distribution are computed simultaneously on a stationary Eulerian mesh grid [30]. PIC methods have been in use since the 1950s. The origin of the PIC method used in the simulation of plasma can be traced to the works of Buneman (1959) and Dawson (1960), as well as Hockney, Birdsall, and Morse [3].
The PIC method is very versatile, applicable to both high density plasmas with collisions (Maxwellian distributions), and low plasma density collisionless systems, such as solar wind and ion thruster simulations [21]. Plasma modelling can be complicated due to the external and induced electromagnetic fields, Coulomb interactions, presence of solid objects, and the different time scales on which electrons and ions are characterized. For an IEC device, I will assume there is no external magnetic field, and that the electric field is static. The current generated by the streaming plasma will be low enough (due to low density of plasma) so that self-induced magnetic field can be ignored [21]. In addition, because electrons have very small masses, their velocities will much larger than those of the ions. In the frame of reference of the ions, the electrons seem to move instantaneously. Thus, electrons will be modelled as a fluid and its density is the Boltzmann distribution, given below.
Eq. 5
This will allow integration time to be much larger, and thus reduce the number of computations required [21].
Calculating the Coulomb interactions between all individual particles lead to a 
computation number, the process would be very computationally intensive. PIC method solves this problem in two ways. Firstly, PIC codes use computational “macro-particles”, where each macro-particle can represent hundreds to millions of real particles. This is a valid approximation as long as the charge to mass ratios remain the same [21]. The number of real particles per macro-particle is called the specific weight. Secondly, instead of computing the Coulomb forces between individual particles directly, the PIC method computes a mesh grid of electric field from contributions of all charged particles in the simulation domain. Via this electric field, the force acting on each particle is calculated. This reduces the computation number from to 
[21].
Design
Figure 10 The basic schematic of an IEC fusion reactor [16].
The physical design of an IEC fusion reactor is shown in Figure 10. My “design” is to simulate such device, varying and iterating several performance parameters so that an optimum design can be found. The overall fusion performance of an IEC device is dependent on many variables. Most of those variables ultimately depend on the operating voltage, pressure, chamber size, and grid transparency [16]. This is shown in the following function.
Eq. 6
To find the optimum design, one therefore needs to take the partial derivative of the above function with respect to each parameter. Since the above function is not defined explicitly, “taking the partial derivative” will be done via numerical simulation. One parameter will be iterated during each simulation while keeping the others constant, then the value for that parameter which gives the maximum fusion performance is the optimum value. These optimum parameters will provide the optimum design for such an IEC device.
Simulation Algorithms
A MATLAB program was written to facilitate the numerical simulation of IEC device using the PIC method. The program was based on a simple MATLAB script developed by Particle in Cell Consulting, originally intended for use in solar wind plasma research. I heavily modified, enhanced, and adapted the template for IEC plasma simulation. This is a two-dimensional PIC algorithm, however with some modifications it will be used to approximate a cylindrical geometry IEC device. A full three-dimensional treatment of the PIC method for a common reactor size would be too computationally intensive. I explored the option of using a commercial 3D PIC plasma simulation software, named VSim. Plasma and fusion experts from VSim Corporation were consulted, and their opinion was that it would take a supercomputer to simulate an IEC device in three-dimension at the sizes and timescales that we are interested in. Therefore I opted for the simpler 2D MATLAB PIC algorithm.
The actual PIC computations occur inside a loop, which is outlined below.
Compute Charge Density
The charge density is a spatially dependent quantity representing the number of unit charges per unit volume. It is computed by distributing the total charge of all particles within a computational cell onto the four nodes of that cell, each dividing by the cell volume. Each cell has a dimension of 
by , this ensures that there will not be any intracellular Coulomb interactions [21].
Figure 11 The simulation cell, containing one macro-particle [21].
A macro-particle is inside a cell, its total charge is the specific weight times the charge per real particle. Charge of the particle is distributed among the four nodes with varying weight factors according to the node’s proximity to the particle. The yellow node receives the most weight while the green node receives the least (Figure 11). The weight factors are given below.
Eq. 7
Eq. 8
Eq. 9
Eq. 10
Where 1, 2, 3, 4 correspond to the blue, yellow, green, and pink nodes, respectively. 
is the fractional distance from the cell origin in the 
direction, and 
is the fractional distance in the 
direction. Exceptions are the walls of the domain, where only half or quarter of the charge is contributing.
Compute Electric Potential
The electric potential is then calculated on the nodes via an elliptic equation solver, which solves the discrete Poisson Equation, given below.
Eq. 11
Finite Difference Method (FDM) with central differencing is used. An elliptic partial differential equation requires boundary values [21]. In the case of an IEC reactor, there are only Dirichlet boundaries (no Neumann boundaries), which specify the values of the potential on the walls and grids.
Compute Electric Field
The electric field is simply the derivative of the potential, as shown.
Eq. 12
For the electric field along boundaries, the distance is halved.
Eq. 13
The component of the electric field is computed in a similar fashion.
Move Particles
Next, the particle motion is integrated through a time step of 
, via the Leapfrog algorithm. The Leapfrog is so named because times at which velocity and position are known are offset by half a time step, and thus they leap over each other (Figure 12).
Figure 12 The Leapfrog method, showing the time offset between position and velocity [21].
First, the electric field at the position of the particle is gathered, then the force and acceleration on the particle is computed. The integration of acceleration yields the velocity, and the integration of velocity yields position. The discrete version of this process is given below.
Eq. 14
Eq. 15
After the particles are moved, a check is performed to determine whether all particles are still in the computational domain. If any particle hits the wall of the chamber or the grids, they will be removed. In the process, they will either emit a secondary electron or simply cause an additional current draw on the power supply [7].
Generate Particles
This step is designed to emulate the non-instantaneous ionization of fuel species. After each iteration, more gas will be ionized. Newly generated ions will be given a thermal velocity, and a Maxwellian position.
Output and Iterate
Output of the simulation includes contour plots of electric potential, electric field, charge density, and ion density. The kinetic information in phase space of all particles are also presented. The different performance parameters of the IEC device are varied for each simulation. These parameters include grid voltage, fuel gas pressure, chamber size, and cathode grid transparency. Plots of the fusion performance and operating conditions are generated for each alteration of the above parameters, some of which are included in the thesis.
Limitations of Simulator
Single Species Model
The simulator uses one fuel species, as opposed to the proposed two species fusion (Deuterium-Tritium). This approximation should not cause a significant impact, as the only difference between a Deuterium and a Tritium ion is mass, which is 1 atomic mass unit [10]. They both behave the same ways to an electric field, and I have used the fusion cross section for DT fusion.
Ionization Approximation
Ionization is a complex process dependent on many variable [18]. In this simulation, ionization is approximated as a gradual process occurring at a specific rate that would eventually cause complete ionization of input fuel.
Geometry Approximation
The simulator uses an extrapolated 3D cylindrical geometry. The base of the simulation domain is a square in the Cartesian coordinates rather than the polar coordinates. (The four sides of a square is charged rather than the circumference of a circle). This is due to convenience of the Cartesian coordinates. I have noticed that this arrangement approximates a polar potential quite well. The third dimension is implemented by scaling the 2D plane. Thus the ion position and velocity distribution is not truly accurate in the third dimension by standards of spherical symmetry. However, this arrangement is quite accurate for a cylindrical IEC chamber. The cell size, which is based on the Debye length, and the number of nodes, specifies the total simulation domain size. Due to difficulties of simulating a large domain (see paragraph below), the cell size is scaled by a number to increase the domain size. This introduces further approximations in terms of ion distributions, but is necessary in order to simulate a meaningful chamber size and pressure without using a supercomputer.
Chamber Size and Pressure Limit
The size of the chamber and fuel pressure (number of particles) is limited in the simulation due to processing speed of the MATLAB script. As noted before, I am using the full kinetic model to simulate the plasma, and thus the kinetic information of every particle is tracked in phase space. This places a processing limit on the number of particles one can simulate. It takes approximately 30 hours to simulate a domain with a 
product of 0.0002 
. Therefore it is unfeasible to satisfy simultaneous requirements of high pressure and large chamber size. As noted before, this type of plasma simulation is best for studying low density plasmas. I had choose this because I wanted to study the thermalization process of Coulomb collisions, and not simply assume a thermalized plasma.
Material Choice
This simulator does not consider any material interactions in a vacuum environment subjected to high intensity ion bombardments. Effects such as surface erosion, sputtering, and contaminations are not considered. However, secondary electron emission is taken into account. Therefore, the optimum choice of material for the cathode grid cannot be determined from this simulator.
Loss Mechanisms
Of the three main loss mechanisms introduced, this simulator takes into account thermalization of ion-velocity and ion-grid collision processes. Bremsstrahlung radiation is not considered here due to the added complexity of the code. I considered this acceptable because Bremsstrahlung radiation is more of an issue for highly quasineutral plasmas with Maxwellian velocity distributions [36], whereas the IEC device has lower degrees of both of the above [36].
Analysis
In each of the following analysis, one parameter is varied in the simulation while keeping all the other parameters constant.
Dependence on Voltage
Figure 13 The dependence of fusion rate as a function of grid voltage. The fuel pressure is kept at 0.001 millitorr, and the chamber radius is 1 m.
Figure 13 is a log-log plot that shows as voltage is increased, fusion rate increases, but only to a point. This agrees with intuition because the fusion cross section as a function of energy also has a peak [10]. The voltage where the peak lies is higher in this plot, which I suspect is due to thermalization and grid collisions causing the bulk of the plasma to not receive the full energy of the potential well as kinetic energy. The optimum voltage is at approximately 190 kV.
Dependence on Cathode Grid Transparency
Figure 14 The dependence of the fusion rate plotted against the cathode grid transparency. Note that this is a log plot. The voltage is at 60 kV, the fuel pressure is kept at 0.001 millitorr, and the chamber radius is 1 m.
Increasing the transparency causes more than a linear performance increase (Figure 14). The cathode grid, if not perfectly transparent, causes ion-grid collisions. This directly reduces the fusion performance by reducing the number of high energy ions, and giving secondary electron emissions [17]. In effect, this causes further thermalization of the ion velocity distribution. Another issue with ions hitting the grid too often is that they cause premature grid corrosion and destruction [16].
Dependence on Pressure
Figure 15 The dependence of fusion rate on the input fuel gas pressure. The voltage is kept at 60 kV, and the chamber radius at 1 m.
Figure 15 shows the effect of varying fuel pressure. As stated earlier, the simulator has limited functionality for higher pressures (above 0.1 millitorr). Therefore the last point in the above plot is probably not accurate. This log-log plot shows that the fusion rate is proportional to the pressure. This result makes sense because the fusion rate equation is dependent on ion densities [10].
Dependence on Size
Figure 16 The dependence of fusion rate on the chamber size. The voltage is at 60 kV, and the pressure is at 0.001 millitorr.
This log-log plot shows that the fusion rate is roughly proportional to the chamber size (Figure 16). This makes sense because if the density of ion (pressure) is held fixed, the larger the chamber the greater the total number of reactions.
Paschen Minimum
From above results, it would seem that to optimize fusion rate, the pressure and size should be as large as possible. In reality, this is difficult to achieve due to Paschen’s Law (Figure 17), which states that for any given voltage, there is a range of that will cause dielectric breakdown [32]. The IEC device (which uses H2) must operate at values of above a point on the left side of the valley or below a point on the right side of the valley. Notice that the graph is asymmetrical, meaning it is easier to satisfy the dielectric requirement with small rather than large . Assuming 1 meter radius, one would need tremendous pressure if operating on the right side of the valley, and the glow discharge becomes an arc discharge.
Figure 17 The Paschen Curve for various gases, illustrating the “Paschen Minimum”. The breakdown voltage of air on either side of the minimum value has inverse relationship with the product of pressure and distance. [32]
The other issue with having high pressures is that Coulomb interactions would prevent any chance of achieving a monoenergetic ion distribution. The ion velocity would be Maxwellian, and the IEC device would in effect become a thermonuclear fusion reactor [10]. Evidently, the pressure needs to be kept low.
Proposed Design
In the regime of our simulations, the optimum and practical IEC fusion reactor would have a chamber radius of 1 m, operating with a pressure of 1 millitorr and a voltage of 190 kV. The cathode grid should be made of tungsten, with a transparency of 98%.
Simulation of My Proposed Design
Voltage: 190 kV
Pressure: 1 milliltorr
Size: 1 m (radius)
Transparency: 98%
Fields and Kinematics
Figure 18 Optimum Design.
Figure 19 Optimum Design
Fusion Rates and Current draw
Figure 20 Fusion Rate of Optimum Design.
Thus, the goal of achieving an IEC reactor capable of producing a neutron flux of 
is actually met.
Conclusions and Recommendations
The original objectives were obviously too broad. I would consider the revised objectives of this thesis complete. A research on fusion mechanics and IEC fusion device in particular was conducted. Various loss mechanisms were investigated, particularly the thermalization of ion velocities and grid collisions. The fusion performance of an IEC device is a function of many variables, and the design was essentially numerical simulations via PIC method of an IEC device, and optimizing several important performance parameters. A preliminary optimum design based on the findings of my simulations was given above.
Improvements in IEC Simulation
Several improvements can be made on the PIC simulator. These include the determination of Bremsstrahlung radiation from Coulomb interactions, charge exchange, and material interactions (grids corrosion). Another possible improvement is to optimize the code and perhaps porting it to C, so that it runs faster.
Improvements in IEC Technology
Ion Microchannels
It was noted in previous research on IEC devices that ions can be made to recirculate in what are termed microchannels between the grid wires, thereby largely avoiding losses due to grid collisions [4, 7, 18].
Multiple Grids
Research has shown that increasing ion confinement time (by reducing grid collisions) can be achieved by focusing the ion paths away from grid wires via insertion of additional concentric grids [7, 17]. These can be thought of as directed ion microchannels.
Polywell
The ultimate goal of directing ions into channels can be achieved via hybrid electric-magnetic confinement, known as a polywell. In a polywell configuration, magnetic mirrors are used to direct and confine electrons to the central region where they act as a virtual cathode. Ions then accelerate towards the virtual cathode, avoiding the issue of grid collisions. Currently, this research is conducted by the US Navy, with high hope of achieving net power.
Word Count: 4852 (±50)
References
- Atzeni, S., & Meyer-ter-Vehn, J. (2004). The Physics of Inertial Fusion: BeamPlasma Interaction, Hydrodynamics, Hot Dense Matter. Clarendon Press.
- Barnes, D., & Nebel, R. A. (1998). Physics of Plasmas 5.
- Birdsall, C. (1991). Particle-in-Cell Charged-Particle Simulations, Plus Monte Carlo Collisions With Neutral Atoms, PIC-MCC. IEEE Transactions on Plasma Science.
- Blair P. Bromley, L. C. (n.d.). Approximate Modeling of Cylindrical Inertial Electrostatic Confinement (IEC) Fusion Neutron Generator. Urbana, IL: Fusion Studies Laboratory, Department of Nuclear Engineering, University of Illinois at Urbana-Champaign.
- Chacón, L., Bromley, B., & Miley, G. (n.d.). Prospects of the Cylindrical IEC Fusion Device as a Neutron Source. Urbana: Fusion Studies Laboratory.
- Charge Exchange (and the creation of energetic neutral atoms). (n.d.). Retrieved from Southwest Research Institute: http://pluto.space.swri.edu/IMAGE/glossary/charge_exchange.html
- Dietrich, C. (2007). Improving Particle Confinement in Inertial Electrostatic Fusion for Spacecraft Power and Propulsion. Massachusetts: Carl C. Dietrich.
- Egle, B. (2006, November). Comparison of Cylindrical and Spherical Geometry in a New Inertial Electrostatic Confinement Device. The University of Wisconsin.
- Farnsworth, P. (1962). U.S. Patent No. 3358402.
- Freidberg, J. (2007). Plasma Physics and Fusion Energy. New York: Cambridge University Press.
- Gu, Y. (1997). Thesis. Department of Electrical and Computer Engineering, University of Illinois.
- Hatzky, R., & Bottino, A. (2010). Particle-in-Cell methods in Plasma Physics. Max-Plack Institut für Plasmaphysik, Euratom Association.
- Hirsch, R. (1967). Applied Physics.
- Keutelian, P., Krishnamurthy, A., Chen, G., Ulmen, B., & Miley, G. (n.d.). Progress in Numerical Simulation of HIIPER Space Propulsion Device.
- Kipritidis, J. (n.d.). Absolute densities of fast H ions in an IEC discharge. The University of Sydney.
- Ligon, T. (1998). The World’s Simplest Fusion Reactor. Analog.
- McGuire, T. (2007). Improved Lifetimes and Synchronization Behavior in Multi-grid Inertial Electrostatic Confinement Fusion Devices. Massachutsetts: MIT.
- Miley, G., Gu, Y., DeMora, J., Stubbers, R., Hochberg, T., Nadler, J., & Anderl, R. (1997). Discharge Characteristics of the Spherical Inertial Electrostatic Confinement (IEC) Device. IEEE TRANSACTIONS ON PLASMA SCIENCE.
- Nadler, J. (1992). Rev. Sci. Instruments.
- Nebel, R., & Barnes, D. (1998). Fusion Technology.
- Particle in Cell Consulting LLC. (2012). The Electrostatic Particle In Cell (ES-PIC) Method. Retrieved from Particle in cell: http://www.particleincell.com/2010/es-pic-method/
- Petkow, D. (2011). On the kinetic modeling of fusion processes in IEC devices. Wiesbaden, Germany: International Electric Propulsion Conference.
- Ren, C. (2011). Introduction to Particle-in-Cell Methods in Plasma Simulations. University of Rochester.
- Rider, T. (1991). A General Critique of Inertial-Electrostatic Confinement Fusion Systems. Massachusetts: MIT.
- Rider, T. (1994). Fundamental Limitations on Plasma Fusion Systems Not in Thermodynamic Equilibrium. Massachusetts: MIT.
- Sarri, G., Murphy, G., Dieckmann, M., Bret, A., Quinn, K., Kourakis, I., . . . Ynnerman, A. (2011). Two-Dimensioanl Particle-In-Cell Simulation of the Expansion of a Plasma into a Rarefied Medium. New Journal of Physics.
- Spirkin, A. (2006). A Three-dimensional Particle-in-Cell Methodology on Unstructured Voronoi Grids with Applications to Plasma Microdevice. Worcester.
- Spirkin, A. (2006). A Three-Dimensional Particle-in-Cell Methodology on Unstructured Voronoi Grids with Applications to Plasma Microdevices. Worcester.
- Tomiyasu, K., Santarius, J., & Kulcinski, G. (2003). Numerical Simulation for UW-IEC Device. Tokyo: Departmen of Energy Sciences.
- Tskhakaya, D., Holger, F., Schneider, Ralf, & Weiße. (n.d.). Chapter 6: The Particle-in-Cell Method, Computational Many Particle Physics. Springer, Berlin Heidelberg.
- Tskhakaya, D., Matyash, K., Schneider, R., & Taccogna, F. (2007). The Particle-In-Cell Method, Contributions to Plasma Physics.
- Paschen, F. (1889). Ueber die zum Funkenübergang in Luft, Wasserstoff und Kohlensäure bei verschiedenen Drucken erforderliche Potentialdifferenz. Annalen der Physik, 69-75.
- Townsend, J. (1910). The Theory of Ionization of Gases by Collision. Constable.
- Townsend, J. (1915). Electricity in Gases. Clarendon Press.
- Bussard, R. (2006). The Advent of Clean nuclear Fusion: Super-performance Space Power and Propulsion. International Astronautical Congress.
- Bussard, R., & King, K. (1991). Bremsstrahlung Radiation Losses in Polywell Systems. EMC2 Technical Report.
Appendix A Derivation of the Fusion Cross Section
A naïve description of the cross section is as follows. Assuming one nucleus is stationary (target) while the other nucleus (incident) is moving towards it. Imagine a spherical force field surrounds the target nucleus. The projection of that sphere perpendicular to the motion of the incident nucleus can be considered the cross section. If the incident particle pass through the cross section area, then the forces exerted between the particles is sufficiently strong to undergo fusion. If the incident particle does not pass through the cross section, then the target force is sufficiently weak that there is no collision [10].
For fusion reactions, which is dominated by short range strong nuclear force, fusion cross section is on the order of a nuclear diameter [10]. However, the actual behavior of the cross section is more complicated than given above.
Hard-Sphere Model
The above description of the fusion cross section is in essence the hard-sphere model. It assumes that each of the colliding particle acts as a billiard ball. The cross section is simply the diameter of the cross section of the nucleus (Figure 2121). This model then gives:
Eq. 16
Figure 21 Definition of the hard-sphere cross section model [10].
This model fails to take into account the fact that slow particles will not undergo fusion because of the repulsive Coulomb force.
Classical Model
The classical model attempts to determine the velocity dependence of by considering the effect of the repulsive Coulomb force and the attractive nuclear force. In the center of mass frame, conservation of energy requires that the sum of the initial kinetic energy of the colliding particles must exceed the Coulomb potential energy when they first come into contact [10]. If this is not satisfied, the particles (D and T) will be repelled rather than fuse (Figure 2222).
Figure 22 The Coulomb force repels, the nuclear force attracts, but only acts at a short range [10].
In mathematical terms, collision occurs when
Eq. 17
which simplifies to
Eq. 18
where
Eq. 19
is the reduced mass, and v is the relative velocity. Substituting the relevant parameters into the above equation yields the minimum center of mass kinetic energy 
, which is quite large, and thus this is a pessimistic model of the fusion cross section (Figure 233).
Figure 23 The classical model of fusion cross section [10].
Quantum Model
The correct model of the cross section involves quantum mechanical effects such as wave-particle duality, resonance effects, and low energy tunneling [10]. Tunneling provides a finite probability for low energy particles to penetrate the potential barrier. Resonance effect occurs under certain conditions of geometry, relative velocity, and potential energy, and results in an enhanced probability of nuclear reaction. Finally, wave-like effect of nuclei sometimes enable them to pass through each other without colliding (two closely coupled waves). The faster the particles’ relative velocity, the greater this effect [10]. These three quantum mechanical modifications on the classical model of fusion cross section give the following results (Figure 24).
Figure 24 The correct model of fusion cross section, including the effects of tunneling, resonance, and high-speed decay [10].
The analytical calculations are usually tedious, therefore actual cross sections are ultimately determined from experimental data, by directing a beam of monoenergetic particles at a stationary target (Figure 2524) [10].
Figure 25 Experimentally measured cross sections σ(v) for various fusion reactions of interest. Notice that DT has a peak of 5 barns at 120 keV, and is clearly favorable in comparison to DD [10].
Velocity Averaged Cross Section
Since the IEC device will at least partially thermalize due to Coulomb interactions, it is necessary to use the Maxwellian velocity distribution to average the fusion cross section [10]. The Maxwellian distribution is given as
Eq. 20
where 
.
With the above distribution function, one will integrate to find the velocity averaged fusion cross section . The following plot is a numerically determined velocity averaged fusion cross section for various fusion fuels (Figure 26 25).
Figure 26 <σv> for various fusion reactions of interest, as functions of temperature. Note that the peak cross section for DT is at 70 keV [10].
Knowing , one can then calculate the fusion output power given as [10],
Eq. 21
Appendix B Additional MATLAB Plots to support Optimization
Voltage
10 kV: Note that the velocity distribution is quite Maxwellian.
90 kV: Note that the ions are much more confined in the central region, and that the distribution is much more monoenergetic (see peak in graph below). This is because Coulomb cross section is relatively smaller at higher energies.
Pressure
0.001 millitorr: Note that the ions are quite sparse and rarefied, leading to a smaller fusion rate (if you can read the contour colour bar).
1 millitorr: The ion velocity distribution seems slightly more Maxwellian, however the overall ion concentration is much higher.
Size
1 mm: The chamber is only 1 mm in radius. Does not affect the ion distribution. But overall fusion rate is lower.
10 m: The chamber size is 10 m, thus the overall fusion rate is much higher. Nothing else changes.
Grid Transparency
85%: You can see that due to grid collisions, the ions are not confined at all to the core, they spread all over the chamber. The velocity distribution also looks bad, as it does not have any peaks (very Maxwellian).
100%: If the cathode grid is perfectly transparent. You can see that the ions are confined to the core perfectly. Also evidently is that the ion velocity is very monoenergetic, (you can see the good looking peak here). For these two reasons, the fusion rate is affected a lot.
