This thesis describes development of the DRACO plasma simulation code. DRACO is an electro-static (ES) code which uses the particle-in-cell (PIC) formulation to track plasma particles through a computational domain, and operates within the Air Force COLISEUM framework. The particles are tracked on a non-standard mesh, which combines the benefits of a Cartesian mesh with the surface-resolving power of an unstructured mesh. DRACO contains its own mesher, called VOLCAR, which is also described in this work.

DRACO was applied to a fully kinetic simulation of an ion-beam neutralization. The thruster configuration and running parameters were based on the NASA's 40cm NEXT ion thruster. The neutralization process was divided into three steps. Electron dynamics was studied by assuming an initial beam neutralization, which was accomplished by injecting both electrons and ions from the optics. Performing the simulation on a full-sized domain with cell size much greater than the Debye length resulted in a formation of a virtual anode. Decrease of the cell size to match the Debye length was not feasible, since it would require a million-fold increase in the number of simulation nodes. Instead, a scaling scheme was devised. Simulations were performed on thruster scaled down by a factor of 100, but its operating parameters were also adjusted such that the produced plasma environment did not change.

Loss of electrons at the boundary of the finite simulation domain induced a numerical instability. The instability resulted in a strong axial electric field which sucked out electrons from the beam. It was removed by introducing an energy based particle boundary condition. Combination of surface scaling and energy boundary resulted in physically sound simulation results. Comparison were made between the Maxwellian and polytropic temperatures, as well as between simulation electron density and one predicted by the Boltzmann relationship.

The cathode was modeled individually from the beam by introducing a positively charged collector plate at a distance corresponding to the beam edge. The local Debye length at the cathode tip was too small to be resolved by the mesh, even if mesh-refinement was incorporated. Since the simulation was not concerned with the near-tip region, two modifications were performed. First, the a limiting value of charge density at the tip was imposed. Second, the cathode potential was allowed to float. These two modifications were necessary to prevent development of a strong potential gradient at the cathode tip.

The modified cathode model was combined with ion injection from the optics to model the actual beam neutralization. Three configurations were tested: a single thruster, a 2x2 cluster with individual cathodes and a similar cluster with a single large neutralizer. Neither of the cases achieved neutralization comparable to one in the base-line pre-neutralized case. The reason for the discrepancy is not known, but it does not seem to be due a loss of electrons at the walls. The difference could be due to limited extent of the modeled physics. An additional work is required to answer this question.