Numerical simulations of the plumes of electric propulsion thrusters, DB VanGilder

Tags: simulations, electron temperature, experimental data, ion, Hall thruster, ion density, plume, collisions, charged particles, particles, Thruster, particle, ions, plasma behavior, c Douglas Bryan VanGilder, physical models, Ion thruster, rectangular cells, node position, Computational techniques, NUMERICAL SIMULATIONS, cell boundaries, distribution function, ion velocity, Current Density, electron number density, collision, Electric Propulsion, Cornell University, Zero dimensional simulations
Content: NUMERICAL SIMULATIONS OF THE PLUMES OF ELECTRIC PROPULSION THRUSTERS A Dissertation Presented to the Faculty of the Graduate School of Cornell University in Partial Ful llment of the Requirements for the Degree of Doctor of Philosophy by Douglas Bryan VanGilder January 2000
c Douglas Bryan VanGilder 2000 ALL RIGHTS RESERVED
NUMERICAL SIMULATIONS OF THE PLUMES OF ELECTRIC PROPULSION THRUSTERS Douglas Bryan VanGilder, Ph.D. Cornell University 2000 The plumes of two electric propulsion thrusters are simulated using a computer code which combines the direct simulation Monte Carlo (DSMC) and the Particlein-Cell (PIC) techniques. The DSMC method used for rare ed gas ow problems is described. The PIC method used to simulate the plasma behavior in the plume is also described. A discussion of the assumptions and the physical models used for the plasma dynamics is included. Details of the collision models are also discussed. Computational techniques permit the physical models and collision models in the simulations to be examined. These techniques include a particle weighting scheme and a routine for simulating expansion into a nite background pressure. Both techniques address resolution problems which occur when species have large di erences in density. The simulations capture the plume behavior qualitatively. The sensitivity of the plume pro le to various physical assumptions and ow conditions is examined. Comparisons with available experimental data reveal which assumptions are valid and which quantities are insensitive to the uncertainties. The electron temperature
and its variation in the plume are important features which a ect the ion pro le. Many of the simulations show good agreement with much of the data. Simulations which include the full chamber of the experimental facility in which the data were obtained lead to good agreement with the data from the centerline to large angles in the thruster back ow region. This level of agreement is not found if the background pressure of the vacuum chamber is not included in the simulations. Discrepancies between the simulations and the experimental data indicate areas where the modeling needs improvement.
Biographical Sketch The author was raised in the northern panhandle of West Virginia. After high school, he went to nearby Bethany College for his Bachelor of Science degree in mathematics and physics. While there, he joined the Phi Kappa Tau Fraternity in 1991, met his future wife, and still found some time to study. He graduated in May of 1994. Afterwards, he ventured to central New York to earn his Doctorate in Aerospace Engineering. In the process, he met some great people and also got married. iii
To my parents, for teaching me the value of an education. And to my wife Jennifer, for sharing in my dreams. iv
Acknowledgements The funding for the work presented in this dissertation was kindly provided by NASA Lewis Research Center through grant NAG3-1451 and by the Air Force O ce of Scienti c Research through Grant F49620-96-1-0091. Computational resources were provided on an IBM SP-2 at the Cornell Theory Center. Thanks to the members of my special committee: Professor Hammer for making PLASMA PHYSICS \mostly" understandable, Professor de Boer for challenging me both in and out of the classroom and for the experimental data, and especially my advisor Professor Boyd for his guidance and advice. He has been a wonderful source for knowledge and aid. He has also made my work at Cornell enjoyable. Also, Professor Caughey was kind enough to ll a vacancy at my defense. I would also like to thank Dr. Brad King and others at the University of Michigan for providing much of the experimental data presented in this work. Many people at Cornell have made my experience more enjoyable. These include Shankar, Tom, and Walt who made me feel welcome as the \new guy" in the o ce (which lasted 3 years), Wyatt and Mike whose stays were brief but quite memorable. Thanks, Mike, for the help in my Japan trip. I owe much gratitude to the following people as well: To Keith and Dan, for helping me make the transition from undergrad and v
answering ALL my questions. Dan for the computer guidance and Keith for whatever else, e.g. \Keith's grid program". Dan, your stories are priceless (even when heard again and again). Keith, thanks for the softball memories. To Jitendra, for his constant cheerfulness and the (last-minute) time we spent on the plasma hw. Why did you bail on the second semester course? To Bill, for being a great lifting partner who made me run EVERY day (not really) and for the screen. I am glad you joined our o ce. Thanks, Jitu and Bill, for reading numerous copies of this document. To Elizabeth, for always having a di ering opinion and for her enlightening sense of humor. Each of you have been wonderful friends and o cemates. To Sharon, for making sure I always knew what I needed to do. To the gang at Binghamton: Dana, Ken, Scott, Martina, Jon and Ti ani, for all the silliness, the cookouts, and the chipping game... Friends like you are appreciated. I'd especially like to say thanks to my family. To my parents, for all their love and support. You always have encouraged my endeavours and taught me to work towards my goals. To my brother Don, who plays the role of a big brother very well. Your accomplishments set a high bar to match. Although, I haven't emulated them, you have been a source of inspiration. Also to my grandparents and my in-laws, for the love and support as I strived to complete my thesis work. And to my wife, for everything. You are an excellent companion. Without your love, understanding, kindness, cooking, etc. this work never would have been completed. Thanks for enduring the stress and for participating in this journey with me while you worked on one yourself. I love you and look forward to spending our lives together. Without these people, this dissertation would never have become a reality. vi
Table of Contents
1 Introduction
1
1.1 Electric Propulsion : : : : : : : : : : : : : : : : : : : : : : : : : : : 1
1.1.1 Electrothermal : : : : : : : : : : : : : : : : : : : : : : : : : 2
1.1.2 Electrostatic : : : : : : : : : : : : : : : : : : : : : : : : : : : 3
1.1.3 Electromagnetic Thrusters : : : : : : : : : : : : : : : : : : : 5
1.2 Motivation : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 7
1.2.1 Contamination : : : : : : : : : : : : : : : : : : : : : : : : : 7
1.2.2 Experiment : : : : : : : : : : : : : : : : : : : : : : : : : : : 8
1.2.3 Computer Modeling : : : : : : : : : : : : : : : : : : : : : : 9
1.3 Approach : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 10
2 Numerical Techniques
11
2.1 DSMC : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 11
2.1.1 Uses : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 11
2.1.2 Collisions : : : : : : : : : : : : : : : : : : : : : : : : : : : : 14
2.1.3 Grids : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 16
2.2 PIC : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 18
2.2.1 Cells : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 18
2.2.2 Equations : : : : : : : : : : : : : : : : : : : : : : : : : : : : 20
2.3 Combining these Numerical Techniques : : : : : : : : : : : : : : : : 23
3 Physical Modeling
25
3.1 Species : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 25
3.1.1 Neutral Atoms : : : : : : : : : : : : : : : : : : : : : : : : : 26
3.1.2 Ions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 29
3.1.3 Electrons : : : : : : : : : : : : : : : : : : : : : : : : : : : : 31
3.2 Equations : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 31
3.3 Boundary Conditions : : : : : : : : : : : : : : : : : : : : : : : : : : 37
3.3.1 Ion Thruster : : : : : : : : : : : : : : : : : : : : : : : : : : : 38
3.3.2 Hall thruster : : : : : : : : : : : : : : : : : : : : : : : : : : 39
3.4 Cathode ow : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 40
3.5 Collision Phenomena : : : : : : : : : : : : : : : : : : : : : : : : : : 40
3.5.1 Neutral-Neutral : : : : : : : : : : : : : : : : : : : : : : : : : 41
3.5.2 Neutral-Ion : : : : : : : : : : : : : : : : : : : : : : : : : : : 41
vii
3.5.3 Ion-Ion : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 44
4 Computational Modeling
48
4.1 Grids : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 48
4.1.1 DSMC Grid : : : : : : : : : : : : : : : : : : : : : : : : : : : 49
4.1.2 PIC Grid : : : : : : : : : : : : : : : : : : : : : : : : : : : : 49
4.1.3 Implicit Scheme : : : : : : : : : : : : : : : : : : : : : : : : : 50
4.2 Particle Weighting Scheme : : : : : : : : : : : : : : : : : : : : : : : 52
4.3 Chamber Back Pressure : : : : : : : : : : : : : : : : : : : : : : : : 55
4.4 Computational Domain : : : : : : : : : : : : : : : : : : : : : : : : : 57
4.4.1 Ion thruster : : : : : : : : : : : : : : : : : : : : : : : : : : : 58
4.4.2 Hall thruster (small) : : : : : : : : : : : : : : : : : : : : : : 60
4.4.3 Hall thruster (full chamber) : : : : : : : : : : : : : : : : : : 62
4.5 Parallel Implementation : : : : : : : : : : : : : : : : : : : : : : : : 62
4.6 Thruster Boundary Conditions : : : : : : : : : : : : : : : : : : : : : 64
5 Ion Thruster Simulations
66
5.1 Inlet Pro le : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 66
5.2 General Pro le : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 71
5.3 Ion Flux Comparisons : : : : : : : : : : : : : : : : : : : : : : : : : 75
5.4 Ion Density : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 81
5.5 Charge-Exchange Ion Pro le : : : : : : : : : : : : : : : : : : : : : : 83
5.6 Electric Field Results : : : : : : : : : : : : : : : : : : : : : : : : : : 89
5.7 Conclusions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 92
6 Hall Thruster Simulations
96
6.1 Inlet Pro le : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 97
6.1.1 Variations in Ion Pro le : : : : : : : : : : : : : : : : : : : : 99
6.1.2 Ion Pro le for Electron Model Comparisons : : : : : : : : : 101
6.2 Results from Sensitivity Study : : : : : : : : : : : : : : : : : : : : : 104
6.2.1 Current Density : : : : : : : : : : : : : : : : : : : : : : : : : 105
6.2.2 Ion Velocity and Density : : : : : : : : : : : : : : : : : : : : 107
6.2.3 Heat Flux : : : : : : : : : : : : : : : : : : : : : : : : : : : : 108
6.3 Results from Modeling Study : : : : : : : : : : : : : : : : : : : : : 114
6.3.1 Electron Density : : : : : : : : : : : : : : : : : : : : : : : : 114
6.3.2 Near Field Current Density : : : : : : : : : : : : : : : : : : 115
6.3.3 CEX Ion Pro le : : : : : : : : : : : : : : : : : : : : : : : : : 119
6.3.4 Full Chamber Geometry : : : : : : : : : : : : : : : : : : : : 120
6.3.5 Potential : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 125
6.4 Parallel Study : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 125
6.5 Conclusions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 129
viii
7 Analysis of Collision Phenomena
131
7.1 Inlet Pro le : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 132
7.2 Simulations : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 133
7.3 Results : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 134
7.4 Conclusions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 144
8 Conclusions
148
8.1 Summary of Conclusions : : : : : : : : : : : : : : : : : : : : : : : : 149
8.1.1 Ion thruster : : : : : : : : : : : : : : : : : : : : : : : : : : : 149
8.1.2 Hall thruster : : : : : : : : : : : : : : : : : : : : : : : : : : 149
8.2 Future Work : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 152
Bibliography
155
ix
List of Tables 1.1 Typical Performance Parameters of Various Electric Propulsion Systems. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 2 2.1 Typical values at the thruster exit for the UK-10. : : : : : : : : : : 12 2.2 Mean free paths for various collisions in the UK-10 ion thruster plume. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 13 2.3 Typical values at the thruster exit for the SPT-100. : : : : : : : : : 13 2.4 Mean free paths for various collisions in the SPT-100 Hall thruster plume. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 14 3.1 Ratio of collision frequency to plasma frequency for each thruster. : 32 3.2 Mean free paths for various collisions in the UK-10 ion thruster plume. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 33 3.3 Mean free paths for various collisions in the SPT-100 Hall thruster plume. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 34 4.1 Typical plasma parameters for the simulations. : : : : : : : : : : : 65 5.1 Table of various UK-10 simulations. : : : : : : : : : : : : : : : : : 71 6.1 Values at the thruster exit for Case 1. : : : : : : : : : : : : : : : : 100 6.2 Ion conditions at thruster exit for sensitivity study. : : : : : : : : : 101 6.3 Values at the thruster exit for reference (Ref) case. : : : : : : : : : 102 6.4 Input conditions for modeling study. : : : : : : : : : : : : : : : : : 103 6.5 Integrated heat ux at both radial locations. : : : : : : : : : : : : 113 6.6 Integrated current (in Amps) at various axial locations. : : : : : : 118 7.1 Values at the thruster exit for reference (Ref) case. : : : : : : : : : 133 7.2 Input conditions for various simulations. : : : : : : : : : : : : : : : 134 x
List of Figures 1.1 Resistojet : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3 1.2 Ion Thruster : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 4 1.3 Hall Thruster : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 7 2.1 Subcells created by particle position along with the weighting factors. 22 2.2 One time step of the DSMC-PIC algorithm. : : : : : : : : : : : : : 24 3.1 Comparison of analytical electron temperature to measurements by Kim. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 35 3.2 Comparison of charge-exchange collision cross sections. : : : : : : : 44 3.3 Scattering angle versus relative velocity for various impact param- eters. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 46 4.1 Implicit grid scheme. : : : : : : : : : : : : : : : : : : : : : : : : : : 51 4.2 Test for background routine. : : : : : : : : : : : : : : : : : : : : : 57 4.3 Computational domain for ion thruster simulations. : : : : : : : : : 58 4.4 Computational grids for UK-10 simulations. : : : : : : : : : : : : : 59 4.5 Ratio of radial grid spacing to local Debye length. : : : : : : : : : 60 4.6 Computational domain. : : : : : : : : : : : : : : : : : : : : : : : : 61 5.1 Gaussian ion ux pro le with near eld measurement. : : : : : : : 67 5.2 Focusing of concave acceleration grid. : : : : : : : : : : : : : : : : 69 5.3 Radial pro les of ion ux at various axial locations for the base case. 72 5.4 Integral ux contours for both electron temperatures. : : : : : : : : 72 5.5 Contours of charge-exchange ions for uniform velocity case with and without back pressure. : : : : : : : : : : : : : : : : : : : : : : : : : 75 5.6 Contours of average axial velocity of thruster neutrals with and without ions. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 76 5.7 Radial pro les of ion ux at various axial locations from measure- ments by de Boer. 10] : : : : : : : : : : : : : : : : : : : : : : : : : 77 5.8 Comparisons of ion ux along axis for various cases with data. 10] : 77 5.9 Comparisons of spherical velocity test case with de Boer's ion ux data 10] at 10 cm and 15 cm from the exit. : : : : : : : : : : : : : 78 5.10 Comparisons of ion ux with de Boer's data 10] at 15 and 25 cm from the exit. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 79 xi
5.11 Comparisons of ion ux from various ion velocity pro les to experiment 5 cm from the exit. 10] : : : : : : : : : : : : : : : : : : : : : 80 5.12 Comparisons of ion ux 25 cm from the thruster exit for various uniform ion velocity cases with data. 10] : : : : : : : : : : : : : : : 81 5.13 Comparisons of ion density to experimental measurements by Pollard 12] taken at various angles 30 cm from the center of the exit. : 82 5.14 Comparisons of ion density to Pollard's data 12] at various angles 61 cm from the center of the exit. : : : : : : : : : : : : : : : : : : : 83 5.15 Comparisons of ion density to Pollard's data 12] at various angles 122 cm from the center of the exit. : : : : : : : : : : : : : : : : : : 84 5.16 Comparisons of charge-exchange ion ux at the thruster exit for various uniform velocity cases. : : : : : : : : : : : : : : : : : : : : 85 5.17 Comparisons of charge-exchange ion ux in the radial direction 20 cm from the thruster exit. : : : : : : : : : : : : : : : : : : : : : 86 5.18 Comparisons of charge-exchange ion radial velocity 20 cm from the thruster exit. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 86 5.19 Normalized axial velocity distribution functions of charge-exchange ions at the beam edge (r = 5cm), 10 cm away from the thruster exit plane. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 88 5.20 Normalized radial velocity distribution functions of charge-exchange ions at the beam edge (r = 5cm), 10 cm away from the thruster exit plane. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 88 5.21 Normalized radial velocity distribution functions of charge-exchange ions away from the beam (r 20cm), 45 cm away from the thruster exit plane. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 89 5.22 Comparisons of the magnitude of ion current density 10 cm from the thruster exit in the back ow region. : : : : : : : : : : : : : : : 90 5.23 Comparisons of plasma potential along axis for various cases with de Boer's data for oating potential 10] : : : : : : : : : : : : : : : 91 5.24 Comparisons of radial electric eld 30 cm away from the thruster exit. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 92 5.25 Comparisons of potential 10 cm above the thruster. : : : : : : : : : 93 5.26 Comparisons of radial electric eld 10 cm above the thruster. : : : 93 6.1 Normalized current density at the thruster exit used as input for simulations. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 103 6.2 Comparisons of ion radial velocity with Manzella's measurements 14] at an axial distance of 11 mm. : : : : : : : : : : : : : : : : : : : : 104 6.3 Comparisons of current density with King's data 17] at radial dis- tances of 50 cm and 1 m. : : : : : : : : : : : : : : : : : : : : : : : 105 6.4 Comparisons of current density with King's data 17] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : 106 6.5 Comparisons of current density with King's data 17] at a radial distance of 1 m. : : : : : : : : : : : : : : : : : : : : : : : : : : : : 106 xii
6.6 Comparisons of average ion velocity with King's data 16] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : 108 6.7 Comparisons of average ion velocity with King's data 16] at a radial distance of 1 m. : : : : : : : : : : : : : : : : : : : : : : : : : : : : 109 6.8 Comparisons of ion density at a radial distance of 50 cm with King's data. 16] : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 109 6.9 Comparisons of heat ux with King's data 16] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 111 6.10 Comparisons of heat ux with King's data 16] at a radial distance of 1 m. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 111 6.11 Comparisons of heat ux with King's data 16] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 112 6.12 Comparisons of heat ux with King's data 16] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 113 6.13 Comparisons of electron density from sensitivity study with measurements by Myers 20] at a radial distance of 31 cm. : : : : : : : : 115 6.14 Comparisons of electron density from modeling study with measurements by Myers 20] at a radial distance of 31 cm. : : : : : : : : 116 6.15 Comparisons of current density from reference case with measurements by Kim. 18] : : : : : : : : : : : : : : : : : : : : : : : : : : : 117 6.16 Comparisons of current density from variable electron temperature case with measurements by Kim. 18] : : : : : : : : : : : : : : : : : 117 6.17 Comparisons of current density from Case 1 with measurements by Kim. 18] : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 119 6.18 Comparisons of ion density for the various simulations at an axial distance 4 cm from the thruster exit in the back ow region. : : : : 120 6.19 Comparisons of axial current density for the various simulations at an axial distance 4 cm from the thruster exit in the back ow region.121 6.20 Comparisons of current density for the various simulations at a radial distance 0.7 m above the thruster. : : : : : : : : : : : : : : : 121 6.21 Comparisons of current density with measurements by King. 17] : : 122 6.22 Comparisons of axial current density 4 cm from the thruster exit (in the back ow region) for the full geometry cases and the smaller domain cases. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 123 6.23 Comparisons of current density 2.425 m above the thruster for the full geometry cases. : : : : : : : : : : : : : : : : : : : : : : : : : : 124 6.24 Comparisons of current density with and without back pressure. : : 124 6.25 Comparisons of potential with measurements by Marrese 21] at an axial distance of 48 cm. : : : : : : : : : : : : : : : : : : : : : : : : 126 6.26 Comparisons of plasma potential with measurements by King 17] at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : 126 6.27 Comparisons of instantaneous potential values at a radial distance of 50 cm. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 128 xiii
6.28 Comparisons of instantaneous axial electric eld values at a radial
6.29
distance of 1 Comparisons
m. of
:::::::: charge-exchange
:X:e+: +:
:::: current
::::: density
: in
:: the
::: back
:
128
ow region. : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 129
7.1 Comparison of reference case with experimental data at r = 50 cm and = 0 and 20 . : : : : : : : : : : : : : : : : : : : : : : : : : : 136 7.2 Comparison of reference case with one allowing ion-ion CEX collisions at r = 50 cm and = 0 . : : : : : : : : : : : : : : : : : : : : 137 7.3 Comparison of case with larger divergence angle with experimental data at r = 50 cm and = 0 and 20 . : : : : : : : : : : : : : : : : 138 7.4 Comparison of case with and without Coulomb collisions at r = 50 cm and = 20 . : : : : : : : : : : : : : : : : : : : : : : : : : : : 139 7.5 Comparison of case with and without Coulomb collisions at r = 50 cm and = 20 . : : : : : : : : : : : : : : : : : : : : : : : : : : : 139 7.6 Comparison of case with larger divergence angle (20 ) to 10 case at r = 50 cm and = 30 . : : : : : : : : : : : : : : : : : : : : : : : 140 7.7 Comparison of cases with larger divergence angle (20 ) to 10 case at r = 50 cm and = 40 . : : : : : : : : : : : : : : : : : : : : : : : 141 7.8 Comparisons of current density from both divergence angle simulations with measurements by Kim. : : : : : : : : : : : : : : : : : : : 142 7.9 Comparisons of current density from cases allowing ion-ion CEX collisions with measurements by King. : : : : : : : : : : : : : : : : 142 7.10 Comparison of various cases at r = 50 cm and = 30 . : : : : : : : 143 7.11 Peaks in distribution function vs. angle for various cases with ex- perimental data at r = 50 cm. : : : : : : : : : : : : : : : : : : : : : 144 7.12 Comparison of reference case with experimental data at r = 50 cm and 1 m and = 150 . : : : : : : : : : : : : : : : : : : : : : : : : : 145
xiv
Nomenclature A = Area At = Area of thruster exit B = Magnetic eld (Tesla) E = Electric eld (V=m) F = Force (N) Ft = Thrust (N) I = Current (A) R = Radius of curvature S = Weight to node for PIC weighting scheme T = Temperature (K unless speci ed as eV ) Wj = Particle weight of species j Xe = Xenon a = Acceleration (m=s2) b = Collision impact parameter dref = Reference collision diameter for VHS model e = Charge of an electron (C) f(V ) = (Normalized) Ion energy distribution function xv
j = current density (A=m2) k = Boltzmann constant (J=K) l = Distance of separation m_ = Mass ow rate m = mass (kg) n = number density (m;3) ne1 = Reference electron number density (m;3), de ned where potential is 0 p = Pressure (N=m2) pb = Background pressure qij = Charge at node (i,j) (C) q = Charge (C) rin = Inner radius of Hall thruster rout = Outer radius of Hall thruster rt = Thruster radius (m) ux = Axial velocity (m=s) v = velocity (m=s) vrel = relative velocity (m=s) vth = thermal velocity (m=s) z = Ratio of ion charge to electron charge dd~Vt = Total derivative (local plus convective rate of change) xvi
Greek Symbols = change in t = time step = Potential drop (V) log( ) = Coulomb logarithm = Neutral atom polarizabilty = cyclotron frequency = Divergence angle = the ratio of speci c heats = Dishing depth (m) 0 = permittivity of free space (F=m) = Perturbation to logarithmic term in potential = E ciency u = Utilization e ciency d = Debye length (m) = Reduced mass (kg) ei = Electron-ion collsion frequency = Mass density (kg=m;3) c = Charge density (C=m3) = Cross section & = Variable used to adjust width of density pro le = collision time (inverse of collision frequency) = Potential (V ) xvii
' = Interaction potential = Collision scattering angle = Double ion fraction ! = Viscosity coe cient exponent !p = Plasma frequency Subscripts 0 = initial a = Neutral xenon back = background cex = charge exchange e = electrons f = fast (high energy) i = ions j = species j m ; t = momentum transfer n = neutrals prop = propellent r = radial direction s = slow (low energy) x = axial direction = azimuthal direction Superscripts + = singly-charged ion ++ = doubly-charged ion xviii
Chapter 1 Introduction 1.1 Electric Propulsion With the increased demand for satellites for various commercial and defense applications, such as communications and global positioning, there is a need for onboard thrusters with long lifetimes. Increasing the mission lifetime of the satellites from 10 years 1] to about 15 years 2] would lead to substantial cost reductions by decreasing the frequency of launches. Compared to chemical propulsion systems, electric propulsion devices o er a higher speci c impulse. As indicated in Table 1.1, these thrusters can achieve an order of magnitude higher speci c impulse than chemical systems. This can lead to substantial savings in the amount of propellent needed and allow satellites to remain in orbit longer. Since these devices also produce considerably less thrust than chemical rockets, they are well-suited for space applications where gravitational forces are weaker than on Earth, and the exterior pressure is very low. Examples of space applications include satellite repositioning, north-south station-keeping, and orbit transfer. Electric propulsion thrusters also have the potential for deep space missions. In fact, an ion thruster has been used for the primary propulsion for NASA's Deep Space One mission. 3] There are several varieties of thrusters that are categorized as electric propulsion devices. The three main categories are described in the following subsections by 1
2 Table 1.1. Typical Performance Parameters of Various Electric Propulsion Systems.
Thruster Thrust Speci c Thruster Input
mN] Impulse s] E ciency %] Power kW]
Chemical rocket 4,000-400,000 200-400
50
100
Resistojet
2-100
200-300
65-90
1
Ion thruster 0.01-200 1500-5000
60-80
1-10
Hall thruster 0.01-2000 2000-5000
30-50
1-30
MPD arcjet 0.001-2000 1000-8000
30-50
0.1-200
way of an example of each type. The classi cation is based on how acceleration is produced. Various performance values for the di erent thrusters are presented in Table 1.1. 2] Approximate values for chemical rockets are included for comparison. Here, the e ciency is de ned as the ratio of thrust-power out to the power put into the system. 1.1.1 Electrothermal Electrothermal thrusters use electrical means to heat the propellent which is then expanded through a converging-diverging nozzle. Perhaps the simplest electric propulsion thruster is the resistojet. A diagram of this thruster is presented in Fig. 1.1. In this device, an electric current is supplied to metal coils. The propellent then ows over these coils and is heated by convection. It is then accelerated gas dynamically using a converging-diverging nozzle.
Gas feed
Input _+ power
Heater
3
Nozzle
Plenum
Expansion
Figure 1.1. Resistojet Several di erent gases have been used as propellent. Ammonia and hydrazine are popular, because when decomposed the molecular hydrogen can be accelerated to high velocities. The velocity is proportional to p1m and hydrogen has a low molecular mass. The speci c impulse is proportional to the velocity, thus light propellents are bene cial. 1.1.2 Electrostatic Electrostatic thrusters use electric elds to accelerate charged propellent that has been ionized inside the thruster. A simpli ed diagram of an ion thruster is presented in Fig. 1.2. Inside ion thrusters, neutral propellent is ionized in an ionization chamber. This can be accomplished either by electron bombardment or by passing the propellent through a hot porous metal (usually tungsten) contact ionizer. In each case positive ions are created, and a plasma is formed inside the chamber. A plasma is a high temperature gas comprised of ions, electrons, and neutrals in proportions such that it is electrically neutral. Downstream of this main discharge plasma, two or three perforated grids are electrically charged. A potential drop on the order of 1000 volts accelerates the ions to high velocities. The rst grid keeps the electrons in
4 Screen Grid
Accelerator Grid
Propellent Feed
Main Discharge Plasma
Ion Flow

e-
Neutralizing Cathode Figure 1.2. Ion Thruster the cavity, while the next grid accelerates the ion beam. Some devices use a third grid with a slightly higher potential than the second grid in order to decelerate and collimate the ion beam. Ion thrusters require an exterior cathode to supply electrons which neutralize the ion beam, because the grid system keeps electrons inside the thruster. If these electrons were not introduced to neutralize the beam, the spacecraft would become charged and space charge e ects would cause the ions to reverse direction or limit greatly the amount of current owing. 4] The hollow cathode which produces a neutral plasma for electron transport is the preferred neutralizer currently used. 5] The propellent most commonly used is xenon, because it is an inert gas with a high atomic weight (131:4 g=mol) and a relatively low ionization potential (12:1 eV ). The ionization potential indicates the amount of energy required to ionize the gas. Thus, a low value allows creation of ions by using less energy than needed for a gas
5 with a high ionization potential. The high atomic weight also improves the energy expenditure per ion. 6] A large mass to charge ratio leads to higher exhaust energy and thrust than a lighter atom for a xed exhaust velocity. The ion thruster considered in the present investigation is the UK-10 thruster which was developed in the United Kingdom and is based on the T5 device. 5] This cylindrical thruster is 10 cm in diameter and uses two (or three) grids which are dished inwards to diminish beam divergence. This thruster is of the Kaufman type, i.e. it uses electron bombardment to ionize the gas inside the ionization chamber. A dc discharge is established between an interior cathode and a cylindrical anode. A magnetic eld provided by solenoids serves to direct the primary electrons for ionization of the neutrals by collisions with these electrons. A more detailed description of this thruster as well as other ion propulsion systems can be found in Ref. 5]. 1.1.3 Electromagnetic Thrusters Electromagnetic thrusters use magnetic elds as well as the electric elds to accelerate highly ionized plasmas. The example considered in the present study is the Hall thruster. Xenon is most commonly used as the propellent for these thrusters as well. Samanta Roy describes these devices as essentially being gridless ion thrusters. 7] The diagram presented in Fig. 1.3, which is similar to one found in Ref. 8], shows the acceleration mechanism of a Hall thruster. These devices have an annular shaped acceleration region with an anode at the base of the thruster. An applied electric eld directed axially originates at this anode and terminates on an exterior cathode at the end of the thruster. This external cathode provides electrons which drift towards the anode to ionize the neutral gas which ows from
6 holes in the anode. Other electrons from the cathode neutralize the ion beam in terms of charge and current. A magnetic eld in the radial direction is applied throughout the annulus using an electromagnet. The electrons interact with this magnetic eld and generate the accelerating electric eld for the ionized propellent. The interaction of the applied electric and magnetic elds gives rise to an E B drift on the ions and electrons. Since the electrons are more mobile than the ions, an azimuthal current is generated, the Hall current. As the electrons drift towards the anode, their axial mobility is reduced by the interaction with the magnetic eld. This reduction in mobility allows the plasma to maintain a very large axial electric eld which accelerates the ions to nearly the full applied potential. In the closed drift region, the magnetic eld must be large enough so that e e >> 1, but i i << 1. 9] A magnetic eld of this magnitude gives the desired electron behavior while having a small in uence on the ions. Otherwise, a large ion azimuthal velocity would diminish the Hall current. A derivation by Shumlak et al. 10] relates the Hall parameter ( e e) to the propellent velocity and the length of the acceleration zone. This relationship can also be used to relate the length of the acceleration chamber to the thruster e ciency. The previous paragraph applies to all Hall current thrusters. The stationary plasma thruster (SPT) is considered in this study. Various SPT's have been own successfully for several years on Russian spacecraft. Perhaps the most widely studied is the SPT-100 where the number indicates the outer diameter of the annulus in millimeters. Typical discharge potentials are 100 ; 700 V . The magnetic eld is generally 100 ; 200 G. Another type of Hall thruster is the anode layer thruster (TAL) which has a much shorter acceleration region. Without wall exchange processes, the electron temperature is very high and ionization and acceleration of
7
Electrons
(Hall current) Gas
Distributor
Ions
Anode
E
B Plasma B Electrons
Cathode Figure 1.3. Hall Thruster these ions occurs in a thin region near the anode. Stronger magnetic elds are used for these thrusters. Another type of electromagnetic thruster is the magnetoplasmadynamic (MPD) thruster. These thrusters use a volumetric j B force to accelerate ionized gases to produce steady continuous thrust. This force results from the interaction of crossed electrical currents in the gas driven by magnetic elds that are external and/or self-induced. 1.2 Motivation 1.2.1 Contamination Although electric propulsion thrusters are well-suited for space applications, there is concern about the interaction of the plasma in the plume with the spacecraft.
8 Sputtering of solar arrays by high energy ions may occur, if the divergence of the ion beam allows ions to angles above 40 degrees from the centerline. These ions not only damage solar arrays, but the sputtered material may be deposited on spacecraft surfaces. The highly-energetic ions may sputter ceramic material from the interior of Hall thrusters and molybdenum from the grids of ion thrusters which may also contaminate the spacecraft. The plasma in the plume is not entirely the high energy propellent. A chargeexchange (CEX) plasma is created from collisions between un-ionized propellent and ions both in the acceleration region and in the plume. In these interactions, electrons are transferred from the neutrals to the ions, resulting in highly-energetic neutrals and slow ions. For ion thrusters, these charge-exchange ions can lead to additional sputtering of grid material into the plume. This not only corrodes the acceleration grids but the material is then available for deposition on spacecraft surfaces. The behavior of the entire plasma and the self-consistent electric elds it creates a ect where this sputtered material, which may also be ionized, is directed. Another problem is that the plasma may interfere with transmissions from the spacecraft. For these reasons, it is important to understand in detail the dynamics of the plume. 1.2.2 Experiment To examine the plumes of these thrusters, experiments are performed in vacuum chambers. These experiments are conducted in order to quantify various plasma properties in the plume. A wealth of experimental data is available for the plumes of the UK-10 ion thruster and the SPT-100 Hall thruster. 11]{ 22] These data include
9 ion and electron densities, current density, ion velocity, electron temperature, heat ux, and plasma potential. However, it is likely that the experimental facilities a ect the plume and thus the properties being measured. The environment the thruster would nd in space is unlikely to be represented. Many vacuum chambers can be maintained at fractions of milliPascals. However, the operating conditions of these thrusters give low plasma densities, and thus interactions between the propellent plasma and the background gas in uence the propellent pro le. This background pressure can lead to higher thrusts and a larger amount of back ow. Also, the walls of the chamber may in uence the plasma. This in uence may be due to the potential of the walls or material sputtered from them. 1.2.3 Computer Modeling Computational modeling allows the dynamics of the plume and interaction with its environment to be examined without any in uences of experimental facilities. Thus, it is possible to emulate the environment the thrusters would experience in space. The ability to simulate the plumes of these devices allows a wider variety of operating conditions to be considered than would be feasible in a laboratory. Computer codes may also allow the e ects of the test facility to be quanti ed, so that future experimental data can be adjusted accordingly. Other investigators have simulated the plumes of electric propulsion thrusters. 7], 23]{ 29] Their investigations are described in more detail in the next chapter. Computer codes must rely on experimental data for validation. The UK-10 ion thruster and the SPT-100 have been chosen for study because of the amount of experimental data available. The goals of the present study are to develop a capability to simulate the plumes of these devices, assess the accuracy of the model through comparison
10 with available experimental data, and determine the sensitivity of the computed results to various parameters. 1.3 Approach To accomplish these goals, a computer code which combines the direct simulation Monte Carlo 30](DSMC) and the Particle-in-Cell 31](PIC) techniques is developed to understand in detail the plasma behavior of the plumes of Hall thrusters and ion thrusters. These two established techniques are well-suited for the rare ed nature of the plumes of these devices. Both techniques employ simulated particles which represent a set number of physical molecules, atoms, ions, or electrons. The PIC method determines the trajectories of charged particles as predicted by imposed and self-consistent electric elds. The DSMC method is used to simulate the collisional e ects in the ow eld. Both charge-exchange and momentum transfer collisions are modeled. Ions and neutral atoms from the thruster and background atoms from the experimental facilities are simulated. The two numerical methods are described in some detail in the next chapter. A description of the physical models used and the associated assumptions are found in Chapter 3. This is followed in Chapter 4 by speci c computational modeling approaches used in this implementation. The inlet conditions employed for each thruster are discussed in Chapters 5 and 6. Also in Chapters 5 and 6, the results for the ion thruster and stationary plasma thruster simulations are presented and compared with data from various experimental measurements. Chapter 7 presents results from a study of the collisional behavior in the SPT-100 plume. Some conclusions and suggestions for future work are included in Chapter 8.
Chapter 2 Numerical Techniques 2.1 DSMC The direct simulation Monte Carlo (DSMC) method is a numerical technique which uses particles for simulating nonequilibrium gas ows. 30] Computational particles which represent real molecules or atoms are tracked through physical space while collisions between them are modeled. Each particle has its own position, mass, velocity, and internal energy. Collisions between these particles can lead to the exchange of momentum and energy, or they may result in chemical reactions. The velocity distribution function is approximated by using these discrete particles. Macroscopic ow properties are then obtained from sampling the microscopic particle properties. The models used to collide the particles and exchange energy are mostly phenomenological in order to represent the physics of the gas ow at the macroscopic level. 2.1.1 Uses The DSMC method is well-suited for nonequilibrium gas ow problems in which binary collisions dominate. Continuum uid mechanics assumes thermal equilibrium locally in the gas. By simulating the gas via discrete particles which represent a large number of real atoms or molecules, nonequilibrium e ects can be captured. 11
12 The Knudsen number is de ned as the ratio of the mean free path, the average distance traveled by a single molecule between collisions, to a characteristic length of the problem. For Knudsen numbers about 10;3 or lower, the continuum approach is valid. For large Knudsen numbers, approximately 10 or higher, the ow is free molecular, thus collisions are extremely rare. For Knudsen numbers in between these regimes, the DSMC technique is most applicable. Examples of ow problems in which DSMC is used include hypersonic shocks, expansion inside the nozzle of a thruster, micro-electrical mechanical systems, and the plumes of thrusters. 32]{ 37] This technique has also been used to simulate the ow and impingement onto a substrate for thin lm deposition. 38] The current study considers the plumes of electric propulsion thrusters. Table 2.1 includes typical densities and velocities for the simulations of the UK-10 ion thruster plume. These values are for the operating conditions and assumptions given in Chapter 5. The mean free paths for various collisions in the ion thruster plume are given in Table 2.2 using the collision cross sections presented in Section 3.5. Values for collisions with the background particles are based on a background density of 7:37 1016 m;3. Table 2.1. Typical values at the thruster exit for the UK-10.
Species Xe Xe+
Density Velocity Temperature
1015m;3] km/s]
K]
822.
0.199
375
6.54
40.2
500
13
Table 2.2. Mean free paths for various collisions in the UK-10 ion thruster plume.
Type/Species Xeprop(m) Xeback(m)
Xem;t Xe+cex Xe+m;t
1.4
16.1
2.6
29.
53.9
600.
Table 2.3 presents typical densities and velocities for the simulations of the SPT-100 Hall thruster plume. These values are for the operating conditions and assumptions given in Section 6.1. The mean free paths for the Hall thruster plume are given in Table 2.4 using the same expressions for the cross sections as used for Table 2.2. Values for collisions with the background particles are based on a background density of 1:5 1018 m;3. The values for both thrusters indicate the rare ed nature of the plumes of these devices. Table 2.3. Typical values at the thruster exit for the SPT-100.
Species Xe Xe+ Xe++
Density Velocity Temperature
1017m;3] km/s]
eV]
11.1
0.325
0.086
1.74
17.6
4
0.58
17.6
4
14
Table 2.4. Mean free paths for various collisions in the SPT-100 Hall thruster plume.
Type/Species Xeprop(m) Xeback(m)
Xem;t Xe+cex Xe+m;t Xe+ce+x Xe+m+;t
1.2
0.92
1.5
1.1
17.5
13.4
0.87
0.67
17.5
13.4
2.1.2 Collisions Particle properties are obtained from a speci ed distribution function at the inlet of the computational domain. Particles are paired for possible collisions. The probability of a collision is a function of the collision cross section and the densities of the species. Various techniques can be used to determine if a pair of particles collide. These include the Time Counter (TC) method, 30] the No Time Counter (NTC) method, 30] and a scheme developed by Bagano and McDonald. 39] The macroscopic collision rate must be reproduced by the instantaneous probabilities of collision. Each technique accomplishes this. The NTC method is given as a computational improvement over the TC method, because the number of collision pair selections is determined before entering the collision coding for a given time step. 30] The probability of collision is obtained by: ( vrvelr)eml ax , where the denominator is the maximum of the product of the collision cross section and the relative velocity. This value is updated if the probability exceeds 1. The Bagano and McDonald scheme is computationally e cient like the NTC method, but it avoids
15 the dependence on the maximum value by using the ratio of the expected number of collisions in a cell during a time step to the number of possible collisions of that type to obtain the probability. This scheme is chosen for the present implementation. The collision dynamics are performed in the center of mass reference frame. Particles may exchange translational as well as internal energy during a collision. Models for exchanging internal energy include rotational 40] and vibrational 41,42] energy. However, the present study only considers translational energy transfer, because the propellent gas is monatomic xenon. The cross sections are functions of collisional energy and depend on the mode of energy to be transferred. In order to update the particles' positions, a time t must be speci ed. Therefore, an algorithmic approach in which a time step is used for each iteration is employed. This time step should be a fraction of the mean time between collisions in order to capture the collisional e ects in the ow properly. The collisions serve to alter the distribution function. The density, momentum, and temperature of the ow are derived from the zero, rst, and second moments of the velocity distribution function respectively. Sampling the computational particles, which is essentially averaging, gives these macroscopic ow properties. For steady ows, the sampling can be extended for a large number of time steps in order to reduce statistical error. It has been shown that this error is a function of the number of particles in a cell and the number of sampling time steps. 43] For unsteady ows, care must be taken to obtain realistic averages. Collisions may also lead to chemical reactions. If the reaction rate is large enough, it is important to include these types of collisions. The charge-exchange reactions mentioned in Chapter 1 are a simpli ed version of chemistry. Other
16 reactions include combination, dissociation, or exchange reactions. Each of these may result in new species being formed. Modeling these types of ows requires additional information concerning possible reaction products and cross sections. The present implementation only considers the charge-exchange reactions because xenon is an inert gas. 2.1.3 Grids The computational domain is represented by a grid comprised of cells. These grids may be either structured or unstructured depending upon the shape of the cells. The cells are used to choose particles for collision and to sample the particles to obtain the macroscopic ow properties. The cells are generally scaled to a fraction of the mean free path, so that only particles within a mean free path are paired for collision. These pairs are chosen randomly in a cell regardless of actual position. The size of the cell is also related to the time step. A given particle should stay in a cell for only about 4-5 time steps to prevent the distribution from being unduly biased by sampling a particle an excessive number of times. The DSMC method can only resolve spatial gradients in the ow eld on the cell size scale. The geometry may be zero dimensional to fully three dimensional. Zero dimensional simulations allow physical modeling and collisional behavior to be examined while the movement is unnecessary. A three dimensional implementation of DSMC is described in Ref. 44]. Fully three dimensional modeling requires considerably more memory and execution time than two dimensional modeling. Not only is more cell information needed, but also more particles must be employed. These two requirements lead to substantially longer simulations. Therefore, the present study assumes symmetry in the azimuthal direction and uses cylindrical coordi-
17 nates. Axisymmetric ows can be represented by a two dimensional grid, but the movement routine uses all three velocity components. The geometry is assumed to be a single azimuthal plane. Particles are moved in all three directions and are then mapped back into this plane. The positions and velocities of the particles are stored in the particle data structures. The particles use their positions to move through the cells according to their velocities. In the present implementation, particles which reach cell boundaries are passed to the neighboring cell. The boundaries of each cell are de ned as either walls, symmetry lines, in ows, out ows, or cells which share an edge (called neighbors). Walls may be fully specular, fully di use, or a combination of the two. Fully specular walls simply re ect the particle back into the simulation without altering properties other than the velocity vector. Fully di use walls accommodate the particles to the wall temperature and sample properties from a distribution before returning them to the cell. It is often necessary to treat some of the collisions between particles and walls as specular and the rest as di use. 30] Symmetry lines act essentially as specular re ection walls. Particles reaching in ow or outow boundaries are removed from the simulation. It is often important to include the facility background pressure in simulations when comparing with experimental data. To represent this e ect, particles can be introduced at the boundaries at a speci ed density and temperature, but this may take a considerable amount of time to reach a steady state. Another approach is to impose ambient conditions uniformly in the cells at the start of a simulation and allow these particles to move through the domain. However, this may require substantially more particles than are desired for a simulation. A scheme which imposes a background density and temperature is described in Chapter 4. Temporary particles are created in the cells
18 for possible collisions, but they are not tracked by the code. 2.2 PIC Analogous to gas ows, plasmas can be modeled either by a uid description or by a kinetic description. Fluid simulations involve solving the magnetohydrodynamic (MHD) equations analogous to computational uid dynamics (CFD) codes which solve the uid equations. Kinetic simulations consider the plasma as particles interacting through electromagnetic elds. The latter can be accomplished by either numerically solving the plasma kinetic equations or by tracking the motion of collections of charged particles via particle simulations similar to the DSMC method. The Particle-in-Cell (PIC) technique uses this particle approach. The PIC method is well-developed, but only the portions of interest to the present implementation are discussed. A detailed description of the method can be found in Ref. 31]. This method has been applied to magnetically and inertially con ned fusion plasmas, electron and ion guns, microwave devices, and plasma propulsion. Work by Roy 7,23] employed the PIC method exclusively to model ion thruster plumes. 2.2.1 Cells Like the DSMC method which uses particles which represent a certain number of real molecules or atoms, the PIC method moves particles which represent neutrals, ions, and electrons through space. Again this domain is represented by cells in a computational grid. Unlike DSMC which de nes the cell properties at the center of the grid cells, PIC uses properties de ned at the nodes of the grid. The reasons for
19
this are made clear in the next subsection which describes the important equations
for the method. The present formulation assumes rectangular cells of various size.
Extension to unstructured grids would alter the relationship between the cells and
the particles.
In simulating plasmas, the cell size, time scale, and number of particles per cell
must be set in order to represent the essence of the plasma physics. Fundamental
properties of plasmas include the Debye length ( d) and the plasma frequency (!p). The Debye length relates the thermal velocity (vth) of the plasma to the plasma fre-
quency. It is the shielding distance around a test charge and the scale length inside
which particle-particle e ects are most signi cant. Beyond this length, collective
e ects dominate.
The plasma frequency is the characteristic frequency of the oscillations of the
electrons about their equilibrium positions. Mathematically, it is:
!p
=
s
ne2 0m
:
The Debye length is related to this parameter by:
d
=
!p vth
:
To update the particles' positions properly according to the physics, the time scale should correspond to the inverse of the plasma frequency. If particle-particle collisions are signi cant in the plasma ow, the cell size should be on the order of or less than the Debye length. Otherwise, in quasi-neutral plasmas in which collective behavior is more signi cant, larger cells can be used.
20
2.2.2 Equations
The PIC technique is a means of tracking charged particles in electric and magnetic
elds. These can be applied or self-consistent elds. The resulting force gives an
acceleration to the particles. The equations of motion are then:
d~v dt
=
F~ m
(2:1)
d~x dt
=
~v:
(2:2)
The following nite di erence equations are used in place of these:
; ~vnew
t
~vold
=
F~old m
(2:3)
; ~xnew
~xold t
=
~vnew
(2:4)
where t is a time step. It is apparent that the new position is incorrect based
on the acceleration. However, if the velocity of the particle \lags" its position by
half of a time step, this error is diminished. This common method for integration
is called the \leap-frog method". The force is then applied at time \t" when the position is ~xold. It is necessary to adjust the initial velocity at t = 0 to ~v(; t=2)
based on the force at t = 0.
The force due to a magnetic eld is a rotation of the velocity. In order to
capture this e ect, the electrostatic force is applied in half time steps. First the
particle is accelerated a half time step. Then the velocity is rotated. Finally, the
other half acceleration is applied to the velocity.
For electrostatic problems the force due to the electric eld, de ned at the
nodes of the grid, comes from:
E~ = ;r :
(2:5)
21
Using central di erencing, the nite di erence equation of this is:
; Ej =
j;1 2
x j+1 :
(2:6)
The potential can be obtained if the charge density is known. The charge density
at the nodes is determined by the particles in each cell. Several schemes can be
used. The simplest is to use the number of particles in a cell and its volume to
determine the density and then assume this value is centered at a particular node.
For a regular grid, this ignores the other three nodes comprising that cell. To better
represent the density, each particle can be split and mapped to the nodes of the cell
it is in. In one dimension, the particle could contribute a fraction of its weight (Wp) to each node based on its proximity to that node. In two dimensions this would
become a weight by area scheme, using the four subcells with the particle position
acting as the fth node. This approach can be extended to three dimensions as
well. For axisymmetric simulations, Ruyten derived a formulation that conserves
both charge and charge density when mapping particles to the nodes. 45] This
scheme is used in the present implementation. Using the weights \S" shown in
Fig. 2.1, the mapping can be described by the following equations:
Sj Sj+1 Si Si+1
= = = =
(xx(xrriixi++j+;+11;11r;;;x;j2)ixx2(x(rrii3()j2rr(+j2j2++1r11j;+;+1rj2r2+)j2r)j3r;j
; r)
r)
(2.7) (2.8) (2.9) (2.10)
The particle is then mapped to each node based on the normalized \area" ratios
of the subcells presented in Fig. 2.1. For example, the product of Si and Sj is the fraction of the particle weight that is assigned to node i,j. Each of these schemes
22 (xi+1,rj+1)
Sj
Sj+1
(x,r) particle
(xi,rj) Si+1
Si
Figure 2.1. Subcells created by particle position along with the weighting factors.
allows for cells of various sizes and assumes a node-centered volume.
Generally, the potential is then obtained by nite di erencing of Poisson's equa-
tion:
; ; j;1
2 j+ ( x)2
j+1 =
c 0
(2:11)
where x is constant in the formulation. However, the present implementation
uses the electron momentum equation to determine the electric eld assuming
quasi-neutrality. This implementation is described and justi ed in Chapter 3.
Once the electric eld is obtained at the nodes of the grid, the e ects on each
particle from its surrounding nodes need to be determined. The resulting acceler-
ation is simply given by mapping the electric eld to the particle using the same
formulation as is used for determining the charge density from the particles. In
fact, this is necessary to prevent a self-force on a particle.
There is an analog for electromagnetic simulations. The current density can be
obtained from particle velocity and position. Determining the electric and magnetic
23 eld then involves solving Maxwell's equations. To make the equations complete, imposed boundary conditions for the potential must be given. Either the potential or the electric eld can be speci ed at the grid boundaries. These are called Dirichlet and Neumann conditions, respectively. The condition applied depends on the computational domain used for the simulations. 2.3 Combining these Numerical Techniques Combining these two techniques seems a natural approach to simulating the plasma plumes of electric propulsion thrusters. Both use computational particles representative of physical particles and move them through space in a computational grid. Each also uses distribution functions to introduce particles at inlet boundaries. Work by Oh and Hastings 24,25] combined PIC and DSMC to simulate Hall thruster plumes. Gatsonis et al. 26] simulated pulsed plasma thrusters with a similar algorithm. Wang et al. 27] modeled ion thrusters using a PIC-MCC (Monte Carlo collision) algorithm. Their work and Roy's 7,23] used an analytical expression to specify the spatial distribution of neutral atoms. A source term was then used for the chargeexchange ions based on an analytical production rate. The present implementation models the propellent neutrals as particles. Charge-exchange ions are generated directly from collisions between these neutral atoms and ions, and the background neutrals with the ions. The algorithm is illustrated in Fig. 2.2 for one time step. The time step is on the order of the inverse of the ion plasma frequency, which is used for ion PIC simulations. This is considerably smaller than the collision time scale generally used
24 for DSMC simulations. Therefore, very few collisions occur in a given time step. The collisions and movement are easily separable. Separate grids are maintained to allow cells to be scaled depending on the physics of interest. The data structures are de ned such that the particles are associated with the cells of the DSMC grid, and movement is performed within these constructs. Particle positions can then be used to map charge density to the correct PIC nodes. This approach has advantages that are discussed in Chapter 4.
PIC Weight particle charges to nodes qij = qion ; nij = qij /e Solve for Electric field as derived from electron momentum equation Weight E-field to particles Eion= Eij Calculate particle acceleration aion= (qion/mion) Eion
DSMC Select collision pairs randomly for every particle Collide particles Change in velocity V
Update particle velocity, V V(t) = Vo + aiont + V Update particle position, X x(t) = xo + Vt
Figure 2.2. One time step of the DSMC-PIC algorithm.
Chapter 3
Physical Modeling There are several issues that need to be addressed to simulate the plumes of plasma thrusters. Computational grids must be generated based on disparate length scales that describe the plasma and collision phenomena. Appropriate boundary conditions for the electric eld must be imposed. The plasma behavior in the plume is determined by the propellent ions and neutrals, the background gas, the electric elds, and the collision dynamics. The following sections describe the important physical aspects of the present implementation. Aspects speci c to computational modeling are discussed in the next chapter.
3.1 Species
The charge-exchange plasma in the plumes of these thrusters is of critical concern
for spacecraft contamination. This plasma is created from collisions between highly
energetic ions and neutral atoms with velocities considerably lower than the ions.
These collisions lead to the transfer of an electron from the neutral to the ion.
These reactions can be described by the following equation:
Xe+f + Xes ! Xef + Xe+s :
(3:1)
Here, the superscript \+" indicates ions of any charge and the subscripts indicate highly-energetic (fast \f") or slow (\s") ions. From this equation, it is apparent
25
26 that not only are the propellent ions important, but charge-exchange ions may also be signi cant. Their signi cance is determined by the neutral xenon ow which a ects the production rate of charge-exchange ions. The other species of interest is the electrons, since this is a plasma ow. The following subsections describe how these species are modeled.
3.1.1 Neutral Atoms
The neutral atoms arise from two sources. Some originate from the thruster as
propellent that is not ionized, while the others constitute the background gas of
the experimental test facility.
The xenon ow is introduced inside the ionization regions of the thrusters as a
neutral gas. Through collisions with electrons, the gas becomes ionized. However,
not all of the propellent is ionized. This residual neutral gas is not accelerated by
the electric eld and, thus, it does not reach a high velocity. Although utilization
e ciencies for these thrusters exceed 90%, the neutral number density is larger than
the ion number density, because the ions are quickly accelerated to high velocities
and extracted. This utilization e ciency is de ned as:
Atmi iX ons(nj (vx)j )
u=
j m_
(3:2)
Neutral velocities are about two orders of magnitude lower than that of the ions, because they are una ected by the accelerating electric eld in the thrusters. Therefore, the neutral propellent density is signi cant, unless this e ciency is nearly 100%. Because these neutrals are found in the ow, charge-exchange reactions can occur inside the acceleration region and in the near plume.
27 Previous computational work on ion thrusters by Wang et al. 27] and by Roy 7] used an analytical expression to specify the neutral number density. Using this density, an analytic source term was then used to create the charge-exchange ions. The present implementation tracks the neutral atoms in the simulations and creates the charge-exchange ions directly from collisions with propellent ions. It is important to represent the neutral density in the plume as accurately as possible, because the production rate depends linearly on this density. Collisions between neutral atoms in which momentum and energy are exchanged are also modeled. These include propellent neutrals colliding with other propellent neutrals or with background neutrals. The cross section for these collisions is a function of relative velocity, and the variable hard sphere (VHS) model is used for the de ection angle. In the VHS model, the de ection angle is chosen randomly from an isotropic distribution. 30] In the center of mass reference frame, the relative velocity is unchanged with collision. Thus, the postcollisional velocities can be obtained from the relative velocity and the de ection angle. The probability of collision is a function of the cross section which depends on the collision type. The collision algorithm is based on the scheme developed by Bagano and McDonald. 39] The functions used for the cross sections are discussed in the nal section of this chapter. Sonic conditions for the neutrals are assumed at the thruster exit based on a stagnation temperature of 500 K for the ion thruster and 1000 K for the SPT. A study by Boyd et al. 37] demonstrates that this assumption of sonic ow agrees well with experiments of neutral xenon ows. It also compares the resulting neutral pro le with that given by the analytical model assumed by Roy. 7] Lower neutral densities are predicted by these simulations. Comparison with experimental mea-
28 surements is not included, because neutral properties are di cult to obtain when the thruster is operating. The background density is large enough that it is di cult to separate thruster neutrals from background particles. In the present study, the neutral properties are assumed to be uniform across the thruster exit plane. The wall e ects on this species are less signi cant than the ions, which are neutralized upon collision with the channel walls. Ground-based experiments have a nite background pressure determined by the capacity of the pumping system. Although these systems can give low densities, there is signi cantly more gas in the chambers when the thrusters are operating. The density of the vacuum chambers is considerably higher than is found in orbit where the thrusters operate on satellites. In low earth orbit (LEO), about 300800 km altitude, the number density varies from 1012 to 1014m;3. The density is orders of magnitude lower for geosynchronous orbit (GEO), about 35,800 km. The densities of the vacuum chambers where the experiments used in the present study were performed are considerably higher. 11,17] The values given in Section 2.1.2 are obtained from using pressure measurements and assuming a room temperature ideal gas. The background densities are 7:37 1016 m;3 for the UK-10 ion thruster and 1:5 1018 m;3 for the SPT-100 Hall thruster. At these operating conditions for the thrusters, the propellent densities are comparable in magnitude to that of the background gas. In regions beyond the near eld, these neutrals are more pronounced than the propellent neutrals. Therefore, the background density must be included in the simulations. It is assumed to be composed entirely of xenon neutrals for three reasons. First, the background pressure is considerably higher when the thruster is operating, suggesting that the gas from the thruster is the main species. Also, the propellent ions are moving at much higher velocities and
29 thus leave the domain more easily. Finally, the ions likely recombine at the walls and become neutrals. The neutrals in the facility collide with ions and neutrals which originate from the thruster. Typical mean free paths of these collisions are given in Tables 2.2 and 2.4. The background particles are not simulated directly, because it is not necessary to update their positions or know their exact properties. Instead, in each computational cell, temporary particles are created each time step with velocities sampled from a Maxwellian distribution at 295 K. The background density is assumed to be uniform. This is reasonable away from the thruster exit where the e ects of the background gas are more signi cant than those of the thruster neutrals. It is not necessary to model collisions between pairs of background atoms. It is assumed that the background distribution is unchanged with collisions. This assumption is reasonable, since the largest change would be due to the fast atoms created by charge-exchange reactions. The density of these fast atoms is about three orders of magnitude below the total facility background density for the operating conditions considered in this study. 3.1.2 Ions The ions are of two varieties: propellent ions from the thruster and charge-exchange ions created by collisions with neutrals. These are treated as di erent species in the computer code. The density of the charge-exchange ions is initially set to zero. The propellent ions may be in various charge states, because electron bombardment of xenon can produce multiply-charged ions. However, these ions are predominantly singly-charged. The fraction of doubly-charged ions ranges from about 10% for ion thrusters 13] to as much as 25% for Hall thrusters. 16] These ions are included in the
30
SPT simulations but not in the ion thruster simulations. Ions with higher charges
have been assumed to be negligible, because triply-charged ions are estimated to
constitute less than 1% of the ions in the plume. 18] The doubly-charged ions are
also of two varieties, propellent ions and the charge-exchange ions created in the
plume.
Speci ed operating conditions of these devices include current, voltage drop,
ow rate, and thrust. Once the fraction of doubly-charged ions is chosen, using
three of these constraints leads to the values for the average ion velocity and
density at the exit plane. Near eld experimental measurements in ion thruster
plumes suggest a Gaussian distribution of current density which varies with radial
position. 11,46] The ion density is lowest at the walls of these thrusters, because
collisions with the walls lead to neutralization of the atoms. Hall thrusters have
longer acceleration regions, therefore this problem is signi cant for them as well.
Therefore, a modi ed Gaussian pro le similar to the one assumed by Wang et
al. 27] for the current density is used at the thruster exit plane for both devices.
This expression terminates at the thruster radius, but is not zero there. The current
density must satisfy:
Z
I = jx(r)dA
(3:3)
where jx(r) = iX ons nj(r)vx j(r)qj]. When a divergence angle is assumed at the
thruster exit, the axial velocity can be given by:
vx(r) = jvj cos (r)]
(3:4)
where (r) is the divergence angle as a function of radial position. The ion number
density can be de ned using a Gaussian function as:
n(r)
=
nmax
exp
";
(r
; &
r)2
#
(3:5)
31 where r is zero for the ion thruster and the radial position halfway between the inner and outer radius for the Hall thruster. The value for & determines the width of the pro le. Various values are used for the simulations. For the ion thruster simulations, & is set to the square of the thruster radius, as is done in Ref. 27]. The details of these assumptions are presented in the chapters which include each thruster's results. 3.1.3 Electrons The di erence in magnitude of the ion and electron thermal velocities makes it di cult to include the electrons as particles in the simulations. Moreover, it is the ion behavior that is of interest in these plumes. If the electrons were modeled as particles, it would require using the electron time scale as the time step. Since this is much smaller than the ion time scale, the simulations would take an inordinate amount of time to reach steady state. Therefore, the electrons are not treated as particles in the simulations. Instead, the plasma is assumed to be quasi-neutral. Hybrid methods which combine uid equations for the electrons with particle methods for the ions are often used to simulate plasmas. 7,31] This approach is simpli ed in this study. The electron uid momentum equation is used to model the electrons. The next section describes this implementation. 3.2 Equations The predominant e ect related to the charge of the ions is the self-consistent electric eld created in the plume of these thrusters. The electrons, being more mobile than the ions, react quickly to this eld and are followed by the ions. Therefore, the
32
governing equation is the electron momentum equation:
mene
dV~e dt
=
;nee(E~
+ ~ve
B~ ) ; rp ; neme ei(~ve ; ~vi):
(3:6)
The term on the left hand side is taken to be 0, because the electrons cannot leave a region in a large group without a large charge imbalance. Thus, the forces must be in close balance. For the densities and temperatures of the plasma typical of the plumes of ion and Hall thrusters, the collision term is negligible. The ratio of collision frequency to ion plasma frequency (!epii lnNNDD ), where ND is the number of charged particles in a Debye cube, is used to gauge collisionality. This ratio is presented in Table 3.1 based on the average ion density at each thruster's exit. Since these values are much less than one, the plasma electrons are considered to be collisionless on the time scale of interest.
Table 3.1. Ratio of collision frequency to plasma frequency for each thruster.
Thruster ne(m;3) Te(eV ) Ratio UK-10 6.54e15 1.0 1.7e-3 UK-10 6.54e15 5.0 1.9e-4 SPT-100 2.5e17 3.0 1.9e-3 SPT-100 2.0e17 3.0 1.8e-3 SPT-100 2.5e17 10.0 3.9e-4
The previous analysis only considers ions exchanging momentum with the electrons. Tables 3.2 and 3.3 are presented to show the signi cance of electron collisions with neutrals for the UK-10 plume and the SPT-100 plume. The mean free paths for these collisions are shown along with those for charged particle collisions using
33 the Coulomb cross section ( CC). The electron-neutral cross sections are estimated from a table presented in Ref. 47]. The values are 2 10;20m2 for Te = 1 eV and 22 10;20m2 for Te = 3 eV . The densities and velocities are from Tables 2.1 and 2.3. The mean free paths are larger for electron-neutral collisions than for the other types presented. Thus, the electron-neutral collisions are neglected. The charged particle collisions do not a ect the momentum of the ions and electrons unless they exchange momentum with one another. Although the electron-ion mean free paths are comparable to the charge-exchange cross sections given in Tables 2.2 and 2.4, the ratio of the frequencies in Table 3.1 indicate that collisions are on a longer time scale than the plasma behavior. Thus, it is reasonable to treat the electrons as collisionless in this analysis. Table 3.2. Mean free paths for various collisions in the UK-10 ion thruster plume.
Type mfp(m)
e;-Xef 60.
e;-Xeb 680.
e;-e; 0.35
e;-Xe+
2.
Xe+-Xe+ 6.4 10;4
The magnetic eld term is generally neglected. 7,25,26,27] This appears to be
a reasonable estimate for the far eld, but its signi cance near the thruster may
be important. 28] This reduces the equation to a balance of the electric eld and
the pressure gradient:
neeE~ = ;rp:
(3:7)
34 Table 3.3. Mean free paths for various collisions in the SPT-100 Hall thruster plume.
Type e;-Xef e;-Xeb e;-e; e;-Xe+ e;-Xe++ Xe+-Xe+ Xe++-Xe++
mfp(m) 4. 3. 0.07 0.66 0.50 0.21 0.04
By treating the electrons as a perfect gas and assuming isothermal conditions, one
obtains the Boltzmann relation:
ne
=
ne1
exp
"
e kTe
#
(3:8)
where ne1 is de ned where the potential is zero. Experimental measurements
of the potential in the plume of each thruster are used to set this value for the
simulations. 11,18] This equation is easily inverted to obtain the potential from
the electron density. To prevent the potential from becoming unbounded, when no
charged particles are in a PIC cell, a small perturbation is introduced by adding a
constant inside the logarithm.
The near eld in the plumes of these thrusters has variations in the electron
temperature. The e ects of including this variation are considered for the Hall
thruster. For simplicity, the potential is not calculated in the simulations. Instead,
the electric elds in the axial and radial directions are calculated directly. Assuming
35
adiabatic conditions, this becomes:
Ex
r
=
;keTe
@ log(ne) @(x r)
;
k e
@(@xTer) :
(3:9)
For the SPT-100, an electron temperature is imposed in the domain based on
measurements taken by Kim et al. 19] The electrons are assumed to act as an expanding uid under isentropic conditions. Their density decays as r;2 and the temperature scales with ne;1. Constants are chosen to closely match these measurements. A comparison in Fig. 3.1 indicates the level of agreement obtained with
this approach. The experiments indicate that the electron temperature does not
continue to decay in the far eld. Thus, it is assumed in the simulations that the
electrons are \frozen" at 1:0 eV or 2:0 eV in various Hall thruster simulations.
Te (eV)
10
Model
9
10 mm
50 mm
8
100 mm
Experiment 7 10 mm
6
50 mm 100 mm
5
4
3
2
1
0
0
0.05
0.1
0.15
0.2
Radial Distance (m)
Figure 3.1. Comparison of analytical electron temperature to measurements by Kim.
The assumption of a quasi-neutral plasma (ni ne) allows the ions to be used to determine the electron density. Based on measurements of oating potential
36
and electron energy in the plume of the UK-10 ion thruster, de Boer estimated the degree of non-neutralization ((ni ; ne)=ne). 11] In the beam region, it is estimated ; to be on the order of 10;5 10;4.
Another complication to the near eld is the contribution of the magnetic eld
to the motion of the electrons. The magnetic eld poses problems in trying to
interpret probe measurements of plasma properties 48] as well as adding complexity
to the computer modeling. For completeness a model which includes these e ects
on the electron behavior is considered.
The electron momentum equation yields three scalar equations.
(x^
:) ; enEx + env
Br
;
@p @x
;n
eimvx
=
0
(3:10)
(r^
:)
;
enEr
;
env
Bx
;
@p @r
;
n
eimvr
=
0
(3:11)
(^ :) ; enE
;
@p [email protected]
;n
eimv
; en(vxBr ; vrBx) = 0
(3:12)
Again assuming symmetry in the direction, the collision term is kept to balance
the Lorentz force. The azimuthal velocity is then:
v
=
;
e(vxBr ; vrBx) m ei
(3:13)
For the axial and radial directions, the collision term is ignored. Substitution of
the azimuthal velocity into the axial and radial equations gives:
; Ex r =
eBr x m ei
(
vxBr
vrBx)
;
1 ne
@(nkT ) @(x r)
(3:14)
The ratio of the cyclotron frequency to the collision frequency (meB ) appears in
the electric eld equation. The collision frequency is a function of the Coulomb
logarithm. This factor represents the maximum impact parameter, as a ratio of the
Debye length to the impact parameter averaged over a Maxwellian distribution. 49]
37 If the Coulomb logarithm is assumed constant over the density and temperature range, this ratio can be manipulated to obtain: C BTn1:5 , where C is a constant. For plasma parameters typical of these thrusters, this logarithm only varies slightly. Unfortunately, the ratio that appears in the azimuthal velocity equation gives unreasonably high velocities assuming a Hall thruster radial magnetic eld near the exit of 200 G (see Section 1.1.3), because the ion-electron collision frequency is very low. For the values given in Table 2.4, an electron density of 2:9 1017m;3 and an axial velocity of 17km=s, and assuming an electron temperature of 3 eV, the azimuthal velocity is about 3:4 107m=s. Because this is a relativistic velocity, this model is not used in the present study. Work by Keidar and Boyd 28] in simulating Hall thruster plumes avoids this problem by starting the simulations slightly after the exit plane. To examine the signi cance of this term, or rather what magnitude a ects the plume, the ratio of cyclotron frequency to collision frequency could be used as a parameter and a magnetic eld pro le imposed in the plume. However, this is beyond the scope of the present study. 3.3 Boundary Conditions To completely de ne the electric eld in the computational domain, boundary conditions must be speci ed. To preserve symmetry, the radial electric eld is set to zero on the plume centerline. For the simulations using Neumann conditions, the far eld boundaries also have the electric eld set to zero in the direction normal to the grid boundary. The thruster exit and the walls of the spacecraft are given a speci ed potential for the simulations which depends on the thruster considered. As previously mentioned, the simulated particles enter the domain at the thruster
38 exit. Particles which reach the outer boundaries are removed from the simulations as are those which cross the in ow boundary, the thruster exit plane. Ions which strike the walls of the thruster are neutralized. Thus, to preserve neutrality, it is assumed that an electron strikes the wall as well. These neutrals as well as the propellent neutrals which strike the walls are accommodated to the wall temperature. This temperature is assumed to be 500 K. The exact value has little e ect on the plume pro le, because it only a ects a fraction of the neutrals. 3.3.1 Ion Thruster Neumann conditions are used exclusively for the simulations of the UK-10 ion thruster plume. The potential at the thruster exit plane is set to -340.0 V corresponding to the accelerator grid potential. The thruster walls are assumed to be biased to spacecraft potential, estimated to be kTe=e according to assumptions by Roy. 7] The plumes are negligibly a ected by small variations in this value, because it only a ects the particles that are near the walls. A singularity occurs at the lip of the thruster. Physically, the thruster wall has some thickness, but the thruster wall meets the thruster exit plane at a point in the simulations. To model this singularity a double boundary condition is used. In the axial direction, the wall potential is assigned to this node. In the radial direction, the node is part of the thruster exit plane and is assigned that potential. This modi cation is important, because the wall potential and thruster grid potential are very di erent.
39 3.3.2 Hall thruster
Most of the simulations of the SPT-100 plume use Neumann conditions for the far eld boundaries. Dirichlet conditions are used for simulations that attempt to model the geometry of the experimental facilities at the University of Michigan. 18] The potential for the outer boundaries is set to zero in these simulations. The di erence in the two domains is discussed in the next chapter. At the thruster exit plane, the potential is assumed to oat with the plasma potential. The spacecraft walls are set to 0 V. Again, the exact value is unlikely to be important. Unlike the ion thruster plume simulations, these simulations model the wall thickness above the thruster opening. The geometry of the thruster splits the exit plane into four distinct regions. These regions include the lower wall from the symmetry line to the inner radius, the accelerating chamber's exit, the wall above this exit, and the region above the thruster. The inner radius is 0.028 m, the outer radius is 0.05 m, and the thruster wall extends to 0.075 m. At the thruster walls, the radial electric eld is zero, and the axial electric eld is obtained from forward di erencing. For the simulations that employ the varying electron temperature model, a potential is required at the neighboring node, if the wall is set to zero potential. For these nodes, the potential is de ned as:
w+1 = Tw+1 log nw+1=ne1 + ]
(3:15)
where the subscript indicates the node next to the wall, and is the perturbation necessary for cells that have no ions. The same technique is used in the radial direction for the outer thruster wall. Nodes at the thruster exit use forward di erencing for the axial electric eld regardless of the model used. This eld varies with plasma potential. The radial
40 electric eld is obtained from central di erencing, as is done in the center of the computational domain. 3.4 Cathode ow The main three dimensional e ects on the domain would be due to the neutralizing cathode which is located o -axis, as seen in Figs. 1.2 and 1.3. However, the domain is assumed to be axisymmetric. Therefore, assumptions need to be made about the cathode ow. For the UK-10, three independently controllable propellent feeds are used. One goes to the main ow distributor, one to the hollow cathode, and one to the neutralizer cathode. 5] The neutralizer is not modeled it receives less than 5% of the xenon ow. 11] The two ows inside the chamber are combined for the total ow rate used in the simulations. For the SPT, the total ow rate is split between the main ow distributor (anode) and the cathode. For the operating conditions of these thrusters, the cathode ow is a small fraction (about 7% 18] to 10% 29]) of the total ow. The ow through the cathode is necessary to obtain the desired emission current. In order to match the total ow rate, this cathode ow is assumed to exit from the thruster as neutral xenon atoms. 3.5 Collision Phenomena The collisions that are modeled can be classi ed as neutral-neutral, neutral-ion, or ion-ion collisions. The collision selection procedure developed by Bagano and McDonald 39] is used to obtain the macroscopic collision frequency of each type of
41 collision. This allows for one pairing per time step per particle. The probability of collision is then a function of the collision cross section. This cross section is a function of relative velocity and depends on the type of collision being considered.
3.5.1 Neutral-Neutral
These types of collisions are either between neutrals from the thruster and back-
ground neutrals or involve only propellent neutrals. The cross section for momen-
tum transfer is given by: 29]
=
2:117 10;18 vr0e:2l4
m2:
(3:16)
The exponent 0.24 is 2! where 0:5 + ! is the exponent for the temperature de-
pendence of the viscosity coe cient of xenon in the variable soft sphere (VSS)
model. 50] This expression is taken from Ref. 50] in which cross section parameters
are obtained in the temperature range 300 to 1500 K.
3.5.2 Neutral-Ion
These collisions can lead to either momentum transfer or charge-exchange reac-
tions. Again the neutrals can be from either the thruster or the background.
Momentum transfer collisions slow down the ions while accelerating the neutrals.
The cross section for these collisions is given by: 51]
=
6:416 10;16 vrel
m2
where the constant includes the e ect of the reduced mass
(3:17) and the polarizability
of the atom. This derivation assumes a long-range attraction potential of the
form:
'=;
e2 l4
(3:18)
42
where l is the separation of the ion and atom. For low temperatures this type
of interaction is important. 52] Thus, for low relative velocities the cross section
becomes large.
Charge-exchange reactions can occur when a highly energetic ion collides with
a neutral atom at a thermal speed. A slow ion and a fast neutral are the result
of this reaction, but the exact postcollisional properties are not clearly de ned.
To capture this behavior, the simulated particles in the collision merely exchange
an electron and maintain their precollision velocities. Although physically there
would be some momentum change, the extent to which the paths are altered is
unknown. Comparisons are made between simulations using this simple charge-
exchange method and one in which elastic hard sphere collisions are assumed and
momentum is transferred between particles. In this case, the scattering angle is
chosen randomly from an isotropic distribution as is done with the non-charge-
exchange collisions.
Regardless of the postcollision dynamics, two di erent cross sections are used
for neutral-ion charge-exchange collisions. Most of the simulations employ the
cross section given by Rapp and Francis for singly-charged xenon ions, 53] which
is based on Dalgarno's derivation. 54] This cross section is:
; RF = ( 0:8821 log vrel] + 15:1262)2 10;20m2:
(3:19)
Rapp and Francis conclude that the cross section for these collisions has three
ranges. 53] 103p m=s
According to their to vrel < 106m=s,
ndings, where
This is in
expression is atomic mass
valid units
in the range (amu). For
vrel > xenon
the lower bound is about 120 m/s. The relative energy must be large enough
that the discreteness of the impact parameter can be neglected in the derivation
and small enough that the transfer of electron momentum can be neglected in the
43
collisions. 54] For doubly-charged ions, an approximate t to Hasted and Hussain's data 55] by Oh and Hastings 29] is used for the reference case. With a correction this is:
; HH = ( 2:7038 log vrel] + 35:0061)2 10;20m2:
(3:20)
The experiments were performed for ion energies in the range of 400 to 3600 eV, therefore extrapolation to the relative energies typical of the plumes considered in the present study may not be accurate. Sakabe and Izawa note a stronger dependence on the ionization potential for resonant charge exchange for various species. 56] The expression:
; SI = ( 21:2 log10 vrel] + 140:0)(I=I0);1:5 10;20m2
(3:21)
where 'I' is the ionization potential and 'I0' is the ionization potential of hydrogen, is found to match experimental data well in the range vrel 5 103 ; 1 108m=s. For gases with other ionization potentials, good agreement is also obtained with data. Rapp and Francis also compared their expression ( RF ) with experimental data. 53] However, Sakabe and Izawa nd better agreement using their expression ( SI) than with RF . 56] The expression SI can then be extended to doubly-charged xenon as well, using its ionization potential. This expression is used for the charge-exchange cross section in Hall thruster simulations designed to analyze the collisional behavior of the plume (presented in Chapter 7). The various charge-exchange cross sections are presented for relative velocities in the range 1 to 150 km/s in Fig. 3.2. The Sakabe and Izawa expression is comparable to the Rapp and Francis expression for singly-charged, but it is signi cantly di erent from the Hasted and Hussein expression for doubly-charged ions.
44
Cross section (Angstroms^2)
300 250 200 150 100 50 0 1000
Sigma_RF Sigma_SI(Xe+) Sigma_HH Sigma_SI(Xe++)
10000 Relative Velocity (m/s)
100000
Figure 3.2. Comparison of charge-exchange collision cross sections. 3.5.3 Ion-Ion The ion thruster simulations and most of the Hall thruster simulations do not include these types of collisions. However, simulations are performed which examine their e ects on the Hall thruster plume. It is found that the ion energy distribution is negligibly a ected by Coulomb collisions between ions. Ion-ion charge-exchange collisions, discussed later in this section, only have a noticeable e ect on the ion distribution if an unreasonably large cross section is used. Experimental measurements obtained using an electrostatic analyzer at the University of Michigan provide information about the ion energy distribution function in the plume of Hall thrusters. 18] This technique provides the distribution as a collected current versus the ion energy per charge. An advantage of particle simulations is that they provide the distribution function directly. They also allow the e ects of various collisions on this distribution to be analyzed. In most of the simulations, the charged particle interactions are determined
45
macroscopically by the electric eld. Beyond the Debye length, collective behavior
dominates, and particle-particle interactions are negligible. However, for smaller
length scales, momentum transfer collisions due to the Coulomb force between
charged particles may be important. These Coulomb collisions are modeled in
some Hall thruster simulations to determine if their e ect is signi cant to the ion
energy distribution function in the plume. The cross section is given by integrating
the Rutherford formula for the di erential scattering cross section for the Coulomb force. 57] The integration is over solid angles with a weighting function (1 ; cos ),
the fractional change of momentum on scattering. This yields:
CC = 4 1 20
z1z2e2 vr2el
!2
log(
)m2:
(3:22)
The logarithmic factor is the familiar Coulomb logarithm which is a function of
plasma properties, including density and temperature. Over a wide range of den-
sities and temperatures typical of Hall thruster plumes, this value is near 12. The
scattering angle is then a function of the energy and the impact parameter of closest
approach. Mathematically, this is:
= 2 tan;1
z1z2e2 ! 4 0b
(3:23)
where 'b' is the impact parameter of closest approach. This value is likely larger
than the neutral xenon reference diameter (dref), the \size" of the atom in the vari-
able hard sphere (VHS) collision model 30] de ned at 273 K. In the VHS model
the diameter of the atom is assumed to be a function of relative velocity related
to the attractive potential. The impact parameter should also be smaller than the
Debye length. Beyond this length, binary collisions of this type are not as signi -
cant as collective e ects in the plasma. Coulomb collisions could allow momentum
transfer from the propellent ions to the charge-exchange ions, thus smoothing out
46
180
d_ref
b
160
Debye
140
120
Chi (degrees)
100
80
60
40
20
0
1
10
100
1000
10000
100000
v_rel (m/s)
Figure 3.3. Scattering angle versus relative velocity for various impact parameters.
the multi-peaked energy distribution function. However, these collisions are predominantly small angle collisions which do not transfer much momentum. A plot of scattering angle versus relative velocity in Fig. 3.3 for various impact parameters illustrates that except for low relative velocities, the scattering angle is small. A value of 5:7 10;9m, or 10dref, is employed for the impact parameter. Based on the range of values and Fig. 3.3, this would tend to exaggerate the e ects of these collisions on the distribution. Collisions between doubly-charged and singly-charged ions may also lead to charge exchange. Because the doubly-charged ions are accelerated more signi cantly by an electric eld than the singly-charged ions, these collisions can a ect the distribution function. Experimental data obtained by King 18] shows high energy tail in the distribution which would not be due to doubly-charged ions, because the ion energy is measured per charge (a voltage). King suggests such signatures are due to singly-charged ions exchanging an electron with doubly-charged
47
ions. Collisions of this type can be described by the following reaction:
! Xe+f + Xe+F+ Xe+f + + Xe+F
(3:24)
where the subscript 'F' indicates a higher energy than the subscript 'f'. These types of collisions would result in a singly-charged ion with a higher energy and a doubly-charged ion with a lower energy. These collisions are modeled to see their e ects on the distribution in some Hall thruster simulations. No references could be found with data or theoretical derivations for the cross sections for these reactions. Therefore, the cross section given by Sakabe and Izawa is used to examine the e ects on the distribution. The same expression is increased by two orders of magnitude in some simulations in order to see a signi cant e ect in the plume.
Chapter 4 Computational Modeling This chapter describes the computational modeling techniques used to aid in simulating the plumes of these thrusters. Since both the PIC and DSMC methods rely on particles which represent real atoms or ions, it is necessary to have enough particles to adequately represent the velocity distribution. The following techniques permit larger simulations than would otherwise be feasible. This allows the physical models described in the previous chapter to be examined. 4.1 Grids The length scale used for PIC cells di ers from that used for DSMC cells. Although assuming a quasi-neutral plasma does not require resolving on the scale of the Debye length, the mean free path is larger than the desired PIC cell size. For the plumes of these thrusters, the Debye length is on the order of 10;5 ; 10;4 m. The mean free paths listed in Table 2.1 are about 5 orders of magnitude higher than this. Therefore, separate grids are maintained for the two numerical methods. While DSMC uses cell-centered values, PIC cells assign values to the nodes. This facilitates the approach, because the grid interaction (via the particles) is not disrupted by using two distinct grids rather than one for both techniques. 48
49 4.1.1 DSMC Grid The underlying code (MONACO 58]) is a parallel DSMC code. The particle movement is performed by this part of the code. This DSMC code can use unstructured grids to represent the computational domain. This approach allows the exibility to scale the size of the cells as the density changes in the ow eld. The grids used for these simulations are triangulated with dimensions signi cantly smaller than the mean free path. They have been created using the advancing front method. 59] 4.1.2 PIC Grid The modularity of the DSMC code allows the PIC routines to be easily imported. These routines include one which maps the particles to the nodes for charge density, another that solves for the electric eld, and a third which accelerates the particles. Each of these routines utilizes the PIC grid. The PIC cells are rectangular, because solving for the electric eld requires nite di erencing of the governing equations. A grid comprised of rectangular cells clearly de nes the nodes relative to one another. For nite di erencing, the spatial derivatives at a node are approximated by di erences in the known variables at the neighboring nodes. 60] In the present implementation, the electric eld is obtained from the ion density and the electron temperature at the nodes. The PIC cells are substantially smaller than the DSMC cells, thus it may require a large number of PIC cells to represent the same domain. De ning the grid implicitly can lead to substantial savings in memory.
50
4.1.3 Implicit Scheme
Although rectangular grids limit exibility of cell size, the redundant information
on node position can be exploited to reduce the memory requirements. In ax-
isymmetric simulations, the nodes of the grid require two position variables, the
node-centered volume, charge density, potential, and two variables for the electric
eld. Using the Boltzmann relation for the potential permits using one variable for
both charge density and potential. Using the variable electron temperature model
does not require the potential, but another variable is needed for the electron tem-
perature. The electric eld vector is a necessary node variable, but the implicit
grid scheme avoids storing the node position.
For rectangular cells of equal size, it is apparent that the node location in the
grid can be de ned by dividing the particle position by the cell size. This can be
described by:
i
=
xion x
(4:1)
j
=
yion y
(4:2)
where \xion" and \yion" are the x and y positions of the particle, and x and y
are the cell dimensions. Thus, i and j de ne the node's axial and radial position
in the grid. This representation requires only two variables, namely x and y,
to de ne the node positions.
The present implementation permits rectangular grids with cells of di erent
sizes. This task is accomplished by using two one-dimensional arrays to store the
node position, one for each direction. Reference lengths ( x)ref and ( y)ref are
de ned. By using multiples of these lengths for the node positions, the locations (i
and j) in the grid can be obtained using the arrays. This method is illustrated in
51
Particle (xion,yion) xref i = xion/xref i' = xloc[i] xloc 0 1 1 2 2 3 3 3 ... (2 1d Arrays)
xi' / dxi' 1 / dxi'
Figure 4.1. Implicit grid scheme. Fig.4.1. The arrays, xloc and yloc, store the node positions using multiples of the reference lengths. A cell two reference lengths wide has two entries in the array, a cell three reference lengths wide would have three entries, and so on. Thus, the node position can be found by:
i0 = xloc i]
(4:3)
j0 = yloc j]
(4:4)
where now i' and j' de ne the node position with i and j given by Equations 4.1 and 4.2. This approach requires only two one-dimensional arrays rather than a two-dimensional array to describe the node positions. The node's spatial location is also needed to determine the weighting factors
52
described in Chapter 2. In the axial direction this is:
; Si+1
=
xion (
xi x)i
(4:5)
where the subscript \i" indicates the node position relative to the grid. Thus, \xi"
is the node position relative to the particle. From this equation, it is apparent
that a node needs a position and a cell length associated with it. Two arrays of
structures with two variables are used to de ne the positions of the grid nodes.
The array for the axial direction is shown in the example in Fig. 4.1. To avoid
performing a division, the cell length is stored as its inverse in these arrays. The
node position is stored normalized by the cell length which eliminates a further
multiplication. This is convenient, because the cell lengths are multiples of the
reference length. The grid le is very compact. It has the number of nodes and the
reference lengths in each direction. These numbers are followed by a list of integers
indicating the axial node spatial position as multiples of the reference length. A
similar list is then given for the radial position. The only other variables necessary
are those used to de ne the thruster walls. The location in the grid of the exit
plane and thruster radius are necessary for the ion thruster. Two radial values are
needed for the Hall thruster simulations along with the exit plane position.
4.2 Particle Weighting Scheme The simulated particles in both numerical methods represent a certain number of actual particles, called the particle weight. This number is set when a particle is generated. The value is generally xed, but may vary with cell, species, or even with particle. For a computational grid with cells of di erent volumes, cell weights can allow a more uniform particle count. 30] This is particularly useful in
53 axisymmetric simulations where the cells near the axis are considerably smaller than those further away. In axisymmetric simulations, radial weighting factors are often used which assign a larger weight to particles further from the axis. However, using two di erent grids to represent the computational domain makes this an unattractive option. This would be di cult to code and would require detailed bookkeeping. The following approach begins as a weight by species method, but uses a weight for each particle. A given particle's weight may vary as a result of collisions in the domain. At the thruster exit, the neutral density is signi cant for ion thrusters and Hall thrusters. In the plume, it is the behavior of the ions that is of primary interest. A su cient number of ion particles are needed in each cell to represent the ion density for the PIC method without introducing arti cial gradients. Therefore, to increase the ion particle count without creating an overabundance of neutrals, a particle weighting scheme is used. The method outlined by Bird for species weighting only conserves momentum and energy on the average, not in each collision. 30] This approach is not conducive to ows with few collisions. Karipides uses an overlay technique to account for species densities of di erent magnitudes. 61] However, this scheme is designed for simulating minor species, those with densities orders of magnitude below the total density. For the plumes of the thrusters modeled in the present study, the neutral density is the dominant species only in the near plume. As the neutrals expand, their density becomes less than the ion density, because the ion beam does not diverge as much. A species weighting method developed by Boyd 62] splits the particle with the larger weight into two parts if it collides with a particle with a smaller weight. One of these has its momentum changed according to the collision dynamics and the two parts are recombined afterwards.
54
However, for collision pairs with vastly di erent velocities, this technique fails to
represent the velocity distribution of the species with the larger weight. Charge-
exchange reactions commonly involve collision pairs of this type. For this reason a new particle weighting scheme is used for the charge-exchange reactions, while Boyd's scheme is used for momentum transfer.
The present scheme is generalized for ion-neutral collisions regardless of which
particle has the larger weight, but in these cases, the neutral either has an equal
or larger weight. The reaction can be described by the following equation:
WnXes + WiXe+f ! (Wn ; Wi)Xes + WiXef + WiXe+s :
(4:6)
Here, the superscript \+" indicates either singly- or doubly-charged ions. The ion
species could have di erent weights as well. This scheme has the drawback that a third particle is created in each charge-exchange reaction. However, the mean free path is large enough ( 1 m at the thruster exit) that this is not a concern. Physically, this shows that only Wi of the Wn atoms undergo a collision and the others are unchanged. This leaves a distribution of Wn ; Wi atoms at the same velocity and the other Wi with a much higher velocity. If these two particles were instead recombined using Boyd's scheme, Wn atoms would have a velocity between the extremes which none of these atoms physically would have. A reference weight is set for all particles in the simulations. The Hall thruster
simulations discussed in Chapter 6 use an initial weight for each species. In mul-
tiples of the reference value, a weight of ten is used for the neutrals from the
thruster and a weight of one for the ions. Unfortunately, the ion thruster simula-
tions presented in Chapter 5 did not have this modi cation to the code. The ratio of neutral atoms to ions at the thruster exit is higher for this thruster than for the
Hall thruster modeled.
55 4.3 Chamber Back Pressure
As mentioned in the preceding chapters, the vacuum chambers in which exper-
iments are performed have a nite back pressure. At room temperature, these pressures give densities of comparable magnitude, 1016 ; 1018m;3, to the propel-
lent densities. In regions beyond the near eld, these neutrals are more prevalent
than the propellent neutrals. Therefore, it is important to model them. Rather
than use valuable computational memory and time storing these particles and mov-
ing them through the domain, a method is developed for simulating their e ect on
the other species in the plume.
As mentioned in Chapter 2, temporary particles are created in the cells for
possible collisions, but they are not tracked by the code. The background particles
are assumed to be entirely xenon neutrals. Their properties are sampled from a
Maxwellian distribution at the speci ed temperature and density.
Not all actual particles in a simulation are paired with these temporary back-
ground particles in a given time step. Otherwise, it would be possible for a particle
to collide twice in a time step. In order to reproduce the collision rate, it is sim-
plest to determine the number of collision pairs as if the background particles were
actual particles. The number of foreground particles paired with other foreground
particles is then given by:
Nfore
=
Nprop
Nprop
Nprop + 0:5
Nback :
(4:7)
In this expression \prop" indicates the actual particles, \back" indicates the back- ground particles, and Nprop ;Nfore is the number of particles available for possible collisions with the background particles. The number of background particles in this calculation (Nback) is given by the density, cell volume, and particle weight in
56 the cell. The number of possible collisions between background and actual particles ; (Nprop Nfore) is less than the number of background particles. This approach assumes collisions between background particles occur, but they are not modeled. Thus, the collision rate is reproduced as if the temporary particles were actually simulated directly, because the correct number of collision pairs of background particles are permitted. The analysis is slightly more complicated for the weight by particle scheme. The temperature of the background gas is used to generate an array that can be used for each velocity component. Temporary particles then sample velocity components from this array. Then, the regular collision routine is used and the sample particle is discarded. If a charge-exchange collision takes place, the temporary particle is put in the simulation as an ion while the original ion is discarded. The e ect of the high-energy neutrals produced in this interaction on the background distribution is neglected. Based on a ux balance of charge-exchange ions and fast neutrals, the density of these neutrals should be about three orders of magnitude below that of the bulk background gas. Therefore, it is unnecessary to perform collisions that would involve only background particles. As a test of this routine, a zero-dimensional simulation was performed. Ions were placed in a single computational cell and permitted to collide with the xed background neutrals. Equilibrium conditions were imposed at a temperatures of 500 K. Allowing only charge-exchange collisions using the cross section of Rapp and Francis, 53] it is simple to calculate the expected increase of charge-exchange ions in the cell. The fraction of charge-exchange ions is plotted versus time step for this simulation and from theory in Fig. 4.2. Good agreement is obtained. This scheme models the background gas without many of the problems inherent
57
Fraction of CEX Ions
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 00
Simulation Theory
500
1000
Time step
1500
Figure 4.2. Test for background routine. in other methods. Introducing particles at the boundaries delays a simulation reaching steady state. Using ambient conditions to initialize the domain with particles may create an overabundance of particles, a memory constraint, and requires additional computational time for movement of these particles.
4.4 Computational Domain The computational domain must include the locations where the experimental data sets used for validation were obtained. The boundary conditions for the domain are described in Chapter 3. Each of the simulations assumes azimuthal symmetry about the thruster centerline. The domain must be large enough that the far eld boundary conditions do not a ect the ow eld signi cantly. Three di erent computational domains are used in this study. Two of these have the thruster
58 located near the upstream edge of the domain with a small back ow region. The other attempts to represent the experimental test facility in which the bulk of the Hall thruster data was obtained. 17,18,22] 4.4.1 Ion thruster The ion thruster computational domain is shown in Fig. 4.3. The back ow region is emphasized along with what is referred to as the region above the thruster in the Chapter 5. The back ow region includes the part of the domain that is in the opposite direction of the main ow from the thruster.
0.7 m
Back flow Region
Above the Thruster
Wall
. m
0.05 m Thruster
Symmetry 1.7 m
Figure 4.3. Computational domain for ion thruster simulations.
The computational domain is comprised of 1,400 DSMC cells and 3,900 PIC cells. It extends 1.3 m axially and 0.70 m radially from the thruster exit. To include some of the back ow region, 20 cm precedes the thruster exit plane. The thruster radius is 5 cm. Presented in Fig. 4.4 are the grids used for the computations in Chapter 5. The upper half is the DSMC grid, and the PIC grid is in the lower half. The DSMC grid has triangulated cells above r = 5 cm. However inside the
59
Radial Distance (m)
DSMC
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0.0
0.5
1.0
1.5
Axial Distance (m)
PIC
Figure 4.4. Computational grids for UK-10 simulations. ion beam, below r = 5 cm, rectangular cells are used for the DSMC grid. This allows the nodes of the DSMC cells to coincide with nodes of the PIC grid. Twice as many PIC nodes in the radial direction are used in this region. The ratio of the radial grid spacing to the Debye length, based on the ion number density, is plotted in Fig. 4.5 for two simulations. One simulation (Grid A) had the same number of nodes as the DSMC grid below r = 5 cm. The other simulation (Grid B) employed the grid with two radial nodes per DSMC node. The ratio is plotted versus radial position from values 5 cm axially from the thruster. The sudden increases in this ratio with increasing radial position occur where the radial cell size is increased from its previous value. Both simulations yield a ratio of less than 10 above the thruster radius. However at lower radial positions, the grid chosen for the simulations (Grid B) is a substantial improvement over the other. Any further reduction in cell size would lead to resolution problems, because the particle count would need a corresponding increase.
60
Ratio y : d
80
70 Grid A
60
Grid B
50
40
30
20
10
00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Radial Distance (m)
Figure 4.5. Ratio of radial grid spacing to local Debye length. 4.4.2 Hall thruster (small) The computational domain for most of the Hall thruster simulations also has the thruster at the upstream boundary. It is assumed to be 14 cm long allowing a back ow region. The geometry of this thruster is di erent from that of the ion thruster. The acceleration chamber is an annulus. The inner radius is 28 cm and the outer radius is 50 cm. Figure 4.6 illustrates the computational domain used. This domain extends 1.3 m axially and 1 m radially. It is represented by 1,600 DSMC cells and 9,500 PIC cells. Although the background pressure of the University of Michigan facility is included, the outer boundaries act as a vacuum where the particles are removed from the simulation. Also shown in this gure is the experimental arc along which King's data were obtained. 17,18] Measurements were taken at various angles along an arc 50 cm and 1 m from the thruster.
61 Vacuum & Er = 0
Vacuum & Ex = 0
Vacuum & Ex = 0 Wall
Wall 0.022 m Thruster m 0.028 m
Data Arc R = 0.5 & 1.0 m Symmetry 1.3 m
1.0 m
Figure 4.6. Computational domain. A new computational boundary is introduced to sample various properties as particles cross these boundaries. The heat ux calculations presented in Chapter 6 and the sampling of the ion energy distribution function presented in Chapter 7 are obtained using this procedure. Neither property is obtainable from the variables sampled in a given cell. Treating this boundary as a wall which allows particles to pass through provides a more exact measure of the ux. To work within the existing data structures of the code, 58] the boundary coincides with cell edges. Particle properties are sampled at these edges depending on their direction of motion. Thus, an edge has two distributions, one for each side. This procedure emulates an experimental probe which would only obtain data from particles going towards it. Particle velocities can be sampled in the cells as is done for the ion thruster simulations, but the boundary scheme gives the distribution of particles in the direction aligned with the probe. This technique is employed for both domains used for the simulations of Hall thruster plumes.
62 4.4.3 Hall thruster (full chamber) Larger Hall thruster simulations include the full geometry of the experimental facility used at the University of Michigan. 18] This geometry is used to better represent the experimental conditions in the simulations and to make comparisons in the back ow region. This domain has a length of 9 m and a radius of 3 m. The thruster is centered radially and located between 2.9 m and 3.1 m axially to match the positioning presented in Ref. 18]. The chamber walls are assumed to be at ground potential. Particles which reach the walls are removed, because the background routine again maintains the facility background pressure. 4.5 Parallel Implementation The UK-10 ion thruster simulations presented in Chapter 5 are performed on one processor. The execution time is about 48 hours for a 500,000 particle simulation on an R10000 SGI workstation using the grid described in the previous section. Although the neutral density at the thruster exit is higher for these thrusters than for the Hall thrusters, the ion density is lower for the ion thruster. The increased ion density from the Hall thruster leads to more charge-exchange ions through collisions with background neutrals. An increased number of particles is needed in order to simulate both ion varieties. The weight by particle technique allows the ion count to be increased. Performing simulations on multiple processors allows the total particle count to increase. The underlying DSMC code 58] is parallelized to run on the IBM SP/SP2 architecture. The computational cells are divided among the processors such that each processor maintains its own collection of neighboring cells. Parallel PIC codes
63 also partition by blocks of cells. 23] However, in the present implementation, separate grids are maintained by the two algorithms. The data structures are de ned such that the particles are associated with the DSMC grid. Using the implicit grid scheme, the PIC algorithm can be run in parallel by a simple tree-sum of the charge density at the nodes. The overhead in this method is that each processor has a copy of the entire PIC grid and must calculate the electric eld for the entire domain. However, these physical models only require a central di erencing scheme for calculating the electric eld. By making this simple change to the PIC part of the code, advantage is taken of the underlying parallelization of the DSMC code. The signi cance of the increased particle count this permits is presented in Chapter 6. Simulations have been performed for the SPT-100 Hall thruster using di erent numbers of processors. Most of the SPT-100 simulations presented in Chapter 6 employ 700,000 to 800,000 particles on an IBM SP2 using the smaller domain described in the previous section. These simulations require about 10 hours in serial followed by 50 to 60 CPU hours distributed over four processors, because the charge-exchange ion time scale is large compared to the propellent ions. These slow ions prolong the time taken to reach steady state. A parallel e ciency is de ned for each processor as a ratio of time spent executing the DSMC-PIC algorithm to the total time on a given processor (communication time + DSMC execution time). In these simulations, this value ranges from 80-90% and varies for each processor. The uctuations are due to the particle weighting scheme. As a result of synchronization during a given time step, the total time required is the same for each processor in the time step. 44] Therefore, an overall e ciency can be de ned from these values as the sum of the individual e ciencies divided by the number of processors. Periodically,
64 a load balancing routine 58,63] is used which distributes the particles more evenly among the processors to improve the e ciency. 4.6 Thruster Boundary Conditions The particles enter the ow domain at the thruster exit plane at a de ned density, temperature, and average velocity that vary with species. The potential at the exit plane is set to -340 V, the accelerator grid potential, for the UK-10 ion thruster simulations. The exit plane potential is allowed to vary with particle density for the SPT-100 Hall thruster simulations. Inverting the Boltzmann relation to obtain the potential requires two values. An electron temperature is estimated for each thruster based on experimental data. A value of 1 eV is used for the ion thruster based on de Boer's plume measurements. 11] The Hall thruster simulations assume 3 eV based on Myers measurements, 21] except in simulations where the electron temperature varies in the plume. The other necessary value is ne1 which is de ned where the potential is zero. This value gives the magnitude of the plasma potential, but it has very little e ect on the di erences in potential except at the boundaries. This reference density is obtained by xing the plasma potential at the thruster exit and using the ion density there. The near eld potential measurements of de Boer 11] are interpolated back to the exit on axis. Estimating a voltage of 7.2 V gives a value of 7:75 1012m;3. Using King's plasma potential measurements 18] to interpolate back to the axis is slightly more di cult. Di erent values are used, as discussed in Chapter 6. These values are on the order of 1015m;3. In each simulation, a small perturbation is introduced by adding 1 10;5 to the term inside the logarithm to
65 Table 4.1. Typical plasma parameters for the simulations. Te ne1 exit d Thruster eV] m;3] V] 10;5] UK-10 1 7:8 1012 7.2 9.2 SPT-100 3 6:5 1015 10.3 2.4 prevent the potential from becoming unbounded when a cell is devoid of ions. The value of ne1 cancels in the ow eld in calculating the electric eld, without the small perturbation. However, at the boundaries it has an e ect. In most simulations the thruster walls are set to a xed potential. Thus, the magnitude of the potential at the neighboring nodes a ects the ions. Ions striking the wall are assumed to be neutralized. These neutrals are then accommodated to the wall temperature. This boundary condition is unlikely to have a noticeable e ect on the plume properties. The charge-exchange ions and slow neutrals are the particles which strike the walls. They have low energy, and therefore contribute little to ux values as they move through the domain. The plasma parameters are included along with the Debye length in Table 4.1 for the simulations described in Chapters 5 and 6.
Chapter 5 Ion Thruster Simulations The main purpose of the ion thruster study is to verify that the computer code has the capability to model the plumes of plasma thrusters. Computational results are compared with available experimental data provided by de Boer 11] and by Pollard 13] for the UK-10 ion thruster. 5] The data from de Boer's plume investigations 11] were obtained for the two grid system, whereas Pollard's data 13] were obtained for the three grid system. The present study also serves to establish the aspects of the modeling that are signi cant in de ning the plume pro le. This is accomplished by testing the sensitivity of the ow eld to various assumptions about physical uncertainties. 5.1 Inlet Pro le The speci ed operating conditions include: ow rate (0.73 mg=s), thrust (18 mN), beam current (0.33 A), and beam voltage (1100 V ). These simulations include only singly-charged xenon ions, because doubly-charged ions only constitute a small fraction (about 10%) of the ions. 13] The magnitude of the ion velocity is obtained by assuming acceleration by the full potential di erence between the grids. Neglecting the velocity of the ions before they reach the screen grid, this velocity 66
67
is:
jvij
=
s
2e mi
which is 40.2 km/s for xenon at a voltage of 1100 V . Using this value along with
the current at the exit de nes the average ion number density as:
ni = ejviIj rt2 which is 6:54 1015m;3 for these conditions. The Gaussian pro le mentioned in
Chapter 3 is applied to the ion density in most simulations. The maximum number density (nmax) in this expression is then obtained by integrating the current density at the exit plane. For a uniform ion velocity this is:
I = 2 j j Z0rte vi nmax exp ;(r=rt)2]rdr:
(5:1)
Ion Flux Ч10-20 (m-2 s-1)
5
4.5
4
3.5
3
2.5
2
Gaussian
1.5
Data 5cm
1
0.5
0-6
-4
-2
0
2
4
6
Radial Distance (cm)
Figure 5.1. Gaussian ion ux pro le with near eld measurement.
68
The ion ux given by this ion density and the uniform ion velocity is presented
with near eld measurements (x = 5cm) of ion ux by de Boer 11] in Fig. 5.1.
The agreement is encouraging. The only signi cant di erence is that the Gaussian
density expression does not decay as sharply at the thruster radius (r = 5cm).
The neutral ux is then obtained from the ow rate and ion ux. Sonic conditions
are assumed with a stagnation temperature of 500 K for the neutral atoms which
sets the neutral velocity and density at the thruster exit. The solenoids used to
generate the magnetic eld inside the thruster chamber heat the thruster walls. 1]
For a stagnation temperature of 500 K, the velocity is 198.9 m/s, and the density is 8:22 1017m;3 for xenon.
Once the maximum ion density and the ion velocity are obtained, the functions
for density and velocity are discretized to be used as input conditions for a xed
number of cells. Ten cells of equal radial spacing are used to represent the thruster
at the exit plane. This technique is used for the following inlet pro les as well.
The UK-10 ion thruster has both the screen grid and the accelerator grid dished
inward to provide mechanical stability and to diminish beam divergence. 5] The
extent to which this curvature a ects the exit pro le of the ions is unknown. Using
an estimated dishing depth of 5 mm, simulations are performed which assume these
grids are spherical or parabolic in shape. The dishing depth is actually 6 mm. 64]
The assumptions are illustrated in Fig. 5.2. The distance \R" is for a spherical
dish, and \p" is for a parabolic dish. Geometrically, these are:
R=
2 + rt2 2
p
=
rt2 4
:
Assuming that the ion bulk velocity vector at the thruster exit is perpendicular
69
rt =5cm = 5mm
x
R or p
Figure 5.2. Focusing of concave acceleration grid.
to the accelerator grid, the ion velocity is given by: vi(r) = jvij cos (r)
where
; cos
(r)
=
r(R(R;
)=rt )2+r2
rt2
for the sphere. This focuses the beam at a focal point 25.25 cm away from the
thruster.
Again the magnitude of the ion velocity is given by the potential di erence,
and the ion number density is obtained by integrating the current density. This is:
I
=2
Z rt 0
evi(r)ni(r)rdr:
Using the de nition of ion velocity, this can be manipulated to give:
ejviIj
rt2
=
Z1 0
p;2;+
s
ni(s)ds
(5:2)
where s (r=rt)2 and ; (R ; )=rt (or (p ; )=rt). Regardless of whether n(r)
is constant at the exit plane or not, the value of nmax can be obtained by relating
this integral to the average ion density, as de ned previously. The same analysis
70 holds for a parabolic grid, where the variable \R" is replaced by \p" for parabolic focusing. The focal point for a parabolic dish is 12.5 cm from the thruster. The background pressure in the vacuum tank has been measured to be 2 10;6 torr. 11] As discussed in the previous chapter, this gas is assumed to be comprised entirely of neutral xenon atoms at a temperature of 295 K. As indicated in Table 5.1, a simulation has also been performed assuming expansion into a vacuum. The base case assumes a Gaussian ion density pro le and a uniform ion velocity at the exit. Charge-exchange collisions do not exchange momentum, and the cross section used is that of Rapp and Francis. 53] The electron temperature is set to 1.0 eV, and the facility back pressure is included as well. Variations in the base case include: assuming an electron temperature of 5.0 eV as opposed to 1.0 eV, allowing momentum transfer for charge-exchange collisions, and using a larger charge-exchange cross section. The electron temperature is based on measurements by de Boer for this thruster which show variations from 0.5 eV to 3.0 eV. 11] This temperature varies with position. A value of 1.0 eV is an approximate average temperature from the data, but a value of 5.0 eV is also considered in order to examine the e ect on the plume. Allowing momentum transfer in charge-exchange collisions has been discussed in Chapter 3. Experimental work by Pollard 14] suggests a charge-exchange cross section for xenon which is as much as 3-8 times higher than that predicted by the Rapp and Francis expression. 53] This would lead to a higher production rate of charge-exchange ions. Therefore, a simulation assuming the theoretical value is compared with one using three times this value for the cross section. The various cases are presented in Table 5.1.
71
Table 5.1. Table of various UK-10 simulations.
Case v(r)
Variation
Base uniform
-
SV spherical
v(r)
PV parabolic
v(r)
Te = 5eV uniform
Te = 5eV
Momentum uniform CEX w/ mom. tr.
CEX uniform
cex = 3 RF
No pb uniform No background pressure
5.2 General Pro le The base case is used as a starting point to test the behavior of the code. Presented in Fig. 5.3 is a plot of the ion ux at various axial positions for this base case. The shape of the pro le is maintained throughout the region, but the magnitude at the centerline decreases along the axis, as the beam spreads. Figure 5.4 shows this contour given as the percentage of ion current enclosed below the radial position. A comparison between this case and the one with an electron temperature of 5.0 eV is presented. Comparison of these two simulations indicates that the ion temperature is not the only cause of spreading. The potential gradients are magni ed by a larger electron temperature for similar density gradients according to the Boltzmann relation. These results are not surprising if one considers the magnitude of the electrostatic forces compared to the hydrodynamic forces. As the ion beam passes
72
Ion FluxЧ10-20 (m-2 s-1)
4
5 cm
15 cm
25 cm
35 cm
48 cm
3
60 cm
2
1
0-0.1
-0.05
0
0.05
0.1
Radial Distance (m)
Figure 5.3. Radial pro les of ion ux at various axial locations for the base case.
Radial Distance (cm)
15 Te = 1.0 eV 10
95
5
90
70
50
30
0
10
10
30
50
-5
70
90 95 -10 Te = 5.0 eV
-15
10
20
30
40
50
Axial Distance (cm)
Figure 5.4. Integral ux contours for both electron temperatures.
73
through a cell in the computational domain, it experiences an electrostatic force in the radial direction, since the gradients in density are steepest in that direction. The ratio of this force to the axially directed force gives a measure of the amount of spreading of the beam across a cell. The electrostatic force across a cell is:
FE = NqE = Nq y where is the change in potential and y is the radial grid node spacing. Here N is the number of ions in the volume represented by the cell, and the other variables have their usual de nitions. The hydrodynamic force is:
Fmv = m_ ux = uAu = nmu2xA
where m_ is the ow rate and ux is the axial component of velocity. The area A is
of the cell face and n is the cell's average ion number density. The ratio of the two
forces is given by:
FE Fmv
=
Nq nmu2xA
y
=
q mu2x
xy :
Using known values of q and m and assuming ux 35 km/s due to deceleration
close to the grid, a value of 6 10;4 y= x is obtained. This indicates that the
bending of the beam is only signi cant for large changes in potential across the
cell. Using the Boltzmann relationship, the change in potential varies as Teln(nn12 ), where n1 is the ion density at the lower node and n2 is the ion density at the upper
node. Inside the beam this density ratio is between 1 and 5 which does not make
a very signi cant potential di erence. With a higher electron temperature, the
spreading is more pronounced.
Simulations are performed with and without back pressure for the uniform ion
velocity pro le. Figure 5.5 shows contours of the density of the charge-exchange
74 ions for the two cases. In both cases, it is clear that their production is signi cant. In the near eld region, charge-exchange ions are accelerated back into the grid leading to grid erosion. In the plume, the ions spread due to the electric elds. These forces are strongest on the edge of the beam where density gradients are largest. The beam ions are moving so quickly in the axial direction that they are a ected much less than the charge-exchange ions. The electrostatic forces pull these charge-exchange ions away from the beam. This is most apparent near the thruster exit where most of the charge-exchange reactions take place. Ions not pulled back into the grid are accelerated away from the beam before they spread much thermally. This explains the lobe a few centimeters from the exit, close to the beam. Also, the uniform contours above the beam are not what would be expected from thermally spreading ions. The e ects of back pressure on the magnitude and location of the charge-exchange ions are signi cant. Therefore, any comparison with experimental data should include the e ects of the background pressure of the facility in which the measurements were obtained. The density is lower along the axis without the back pressure, because very few charge-exchange reactions take place in the far plume. Also, fewer charge-exchange ions reach the region behind the thruster exit when the background pressure is not included. The charge-exchange collisions also a ect the neutral distribution. The in uence of the creation of fast neutrals on the distribution of neutral particles from the thruster is presented in Fig. 5.6. This shows a plot of the average axial velocity of the thruster neutrals. A simulation has been performed with only the neutral ow conditions at the thruster exit. This gas expands into the domain with the xed background pressure. The results from this simulation are presented in the lower half of Fig. 5.6, whereas the upper half of the gure is from the base case. The
75
Radial Distance (m)
0.6
pb
2
0.4
1
0.2 0.0 -0.2
76
3
4 6
5
4
3
2
-0.4
1
No pb
0.0 0.2 0.4 0.6 Axial Distance (m)
Lv. Xe+ 7 1E15 6 5E14 5 3E14 4 1E14 3 5E13 2 1E13 1 5E12 0.8
Figure 5.5. Contours of charge-exchange ions for uniform velocity case with and without back pressure.
increase in average velocity for the base case indicates that a substantial fraction of the neutrals have undergone a charge-exchange reaction. A comparison of neutral density shows a corresponding decrease for the case with ions, but the e ect is less pronounced.
5.3 Ion Flux Comparisons Experiments by de Boer provide ion ux pro les at various axial positions. 11] These measurements were obtained using an electric probe consisting of two parallel cylindrical electrodes in the direction perpendicular to the ow. Interpreting the current-voltage characteristics obtained by varying the potential between the electrodes provides the ion ux and electron temperature. This ion ux data is presented in Fig. 5.7. The shape of this pro le resembles the simulated results, but with increasing axial distance discrepancies between simulation (shown in Fig. 5.3) and experiment become more apparent. Experimentally, the peak of the ion ux
76
1
0.6
2
Radial Distance (m)
0.4 0.2 0.0 -0.2
2 2 1
2
44
3 34
3
3
Lv.
Vx (m/s)
6 5000
5 3000
3 5
4
6
4 3
1000 500
2 300
1 100
1 -0.4
-0.6
1
2 2 2
0.0
0.5
1.0
Axial Distance (m)
Figure 5.6. Contours of average axial velocity of thruster neutrals with and without ions.
increases with axial distance until about 15 cm and then decreases. Contours of this ux show a waist near 15 cm from the thruster exit. This e ect is caused by the focusing of the curved acceleration grids of the UK-10. 11] Presented in Fig. 5.8 is the ion ux along the axis for various simulations compared with the experimental data. In this gure, the peak in the data is clear. The simulations with a uniform ion velocity at the exit plane clearly do not reproduce this e ect. Therefore, simulations which describe the curvature geometrically, according to the previous analysis, are performed. Comparisons between the simulation which assumes the spherical velocity prole and ion ux measurements at two axial locations near the waist of the beam are shown in Fig. 5.9. By 15 cm the simulated peak in ux is about an order of magnitude higher than the experiment. This location is before the imposed focal point. Also, the sharp peaks near the axis indicate that the focusing of the thruster grid is much less pronounced than a fully spherical grid predicts. Similar results are obtained for the parabolic pro le simulation. For both pro les, the bulk of the
77
8
5cm
7
15cm
25cm
35cm
6
48cm
60cm
5
Ion Flux Ч10-20 (m-2 s-1)
4
3
2
1
0-10
-5
0
5
10
Radial Distance (cm)
Figure 5.7. Radial pro les of ion ux at various axial locations from measurements by de Boer. 10]
Beam Ion Flux Ч10-20 (m-2 s-1)
10
9
Base Te = 5eV
CEX
8
Momentum
7
No pb Data
6
5
4
3
2
1
0
0
0.25
0.5
0.75
1
Axial Distance (m)
Figure 5.8. Comparisons of ion ux along axis for various cases with data. 10]
78
6
SV 10cm
SV 15cm
5
Data 10cm
Data 15cm
4
Ion Flux Ч10-21 (m-2 s-1)
3
2
1
0
-0.05
0
0.05
0.1
Radial Distance (m)
Figure 5.9. Comparisons of spherical velocity test case with de Boer's ion ux data 10] at 10 cm and 15 cm from the exit.
ions reach the axis, and the waist is therefore very narrow. Hence, it is concluded that the paths of the ions are not determined exactly by the dish geometry. Information about the actual pro les at the thruster exit plane would require an accurate representation of the ow inside the thruster. Various ion optics simulation studies have recently been performed which consider the plasma ow in the inter-grid region of ion thrusters. 65,66] However, these studies have primarily focused on the grid erosion. Lacking information on this exit plane pro le for the ions, attempts are made to better represent the pro le of the ions in the plume. The true exit pro le must lie somewhere between a fully focusing grid and a uniform velocity pro le. To help quantify this departure from full focusing, the spherical and uniform velocity pro les are superimposed. One simulation combines a 5% spherical velocity pro le with 95% uniform velocity where the percentage is by ion density. Another simulation uses 10% spherical and 90% uniform velocity. This
79 new expression is integrated as well to obtain appropriate values for the maximum ion density (nmax). Presented in Fig. 5.10 are results from the 5% case along with the base case and de Boer's experimental data at 15 and 25 cm from the thruster exit. The 5% case agrees well in magnitude with the experimental data at these locations, but it has a very pronounced peak as seen in the fully spherical velocity simulation. The uniform velocity case agrees well in shape with the experiment. Comparisons of the various simulations with the near eld experimental data at 5 cm axially from the thruster are presented in Fig. 5.11. Good agreement with the experiment is obtained for each case, except for the fully spherical case. Therefore, the inlet ion pro les are not unreasonable. These simulations demonstrate that the beam essentially follows its initial velocity pro le, as predicted by the ratio of electrostatic to hydrodynamic forces.
10
Base 15cm
Base 25cm
5% SV 15cm
8
5% SV 25cm
Data 15cm
Data 25cm
6
Ion FluxЧ10-20 (m-2 s-1)
4
2
0
-0.05
0
0.05
0.1
Radial Distance (m)
Figure 5.10. Comparisons of ion ux with de Boer's data 10] at 15 and 25 cm from the exit.
80
7
Base
5% SV
6
10% SV
SV
Data
5
Ion FluxЧ10-20 (m-2 s-1)
4
3
2
1
0
-0.05
0
0.05
0.1
Radial Distance (m)
Figure 5.11. Comparisons of ion ux from various ion velocity pro les to experiment 5 cm from the exit. 10]
Figure 5.12 shows a typical result for the beam ion ux for the di erent uniform velocity cases. The momentum transfer case is nearly identical to the base case, because the depletion of the beam ions by charge-exchange reactions is unaltered. The larger cross section case (CEX) is only slightly lower due to a higher chargeexchange collision rate. Closer to the thruster exit the di erences are even less signi cant. This comparison along with Fig. 5.8 indicates that the ion beam is not changed signi cantly by these parameters.
81
6
Base
CEX
5
Momentum
Data
4
Ion FluxЧ10-20 (m-2 s-1)
3
2
1
0
-0.05
0
0.05
0.1
Radial Distance (m)
Figure 5.12. Comparisons of ion ux 25 cm from the thruster exit for various uniform ion velocity cases with data. 10]
5.4 Ion Density
Thus far, the comparisons have been for ion ux, which includes ion density and velocity. Di erences may be found for these constituent variables. Experimental work by Pollard provides plume measurements of ion density at various angles for xed distances from the thruster exit. 13] This data was obtained for a slightly higher thrust (20 mN) using a Langmuir probe. Comparisons of ion density allow a closer examination of the ion pro le. Results from various uniform velocity simulations are presented in Fig. 5.13 with Pollard's data at 30 cm from the thruster exit. Qualitatively, the simulation results agree quite well with the measured data. Figures 5.14 and 5.15 show comparisons at 61 cm and 122 cm from the exit. In each case, the uniform velocity pro les overpredict the ion density near the axis. Consistent with these peaks, the experiment shows more spreading at large angles
82
1016 1015
Base Te = 5eV CEX Momentum Data
Ion Density (m-3)
1014
1013
1012
-40
-20
0
20
40
Angle (deg)
Figure 5.13. Comparisons of ion density to experimental measurements by Pollard 12] taken at various angles 30 cm from the center of the exit.
away from the axis. The e ect on the charge-exchange ions is apparent in these plots. They dominate the ion density away from the beam in the wings of the pro les. As expected the ion density is highest for the case with a larger chargeexchange cross section (CEX). This case and the base case agree better with the experiment away from the axis than the others. The momentum transfer case predicts a lower density, because the postcollisional velocity of the charge-exchange ions is substantially higher for this type of collision. As with the ion ux comparisons, these ion density comparisons indicate that the beam ions are only slightly a ected by changing the physical parameters. The larger electron temperature predicts more beam spreading (as seen in Fig. 5.4), but the agreement with ion density is not as good at the larger angles. At angles close to the thruster centerline, the higher electron temperature case shows better agreement with data than the other
83
1016 1015
Base Te = 5eV CEX Momentum Data
Ion Density (m-3)
1014
1013
1012
-50
0
50
Angle (deg)
Figure 5.14. Comparisons of ion density to Pollard's data 12] at various angles 61 cm from the center of the exit.
cases. The variable electron temperature model described in Section 3.2 may yield better agreement at all angles, because the actual electron temperature varies in the plume. 11] This model is used for a set of Hall thruster simulations described in the next chapter. An advantage of this approach, which will be discussed further in the next chapter, is that a higher temperature can be used in the beam, while a lower temperature is used further away. This better represents the domain, as indicated by electron temperature measurements by de Boer. 11]
5.5 Charge-Exchange Ion Pro le Although the e ect on the beam ions is not very signi cant, changing the physical parameters does a ect the charge-exchange ion pro le. The di erence in chargeexchange ion ux near the thruster exit is illustrated in Fig. 5.16. The charge-
84
11016 11015 11014
Base Te = 5eV CEX Momentum Data
Ion Density (m-3)
11013
11012
-40
-20
0
20
40
Angle (deg)
Figure 5.15. Comparisons of ion density to Pollard's data 12] at various angles 122 cm from the center of the exit. exchange ion ux is two orders of magnitude below that of the beam ions. This suggests that the measured ux is almost entirely beam ions. The behavior of the charge-exchange ions near the thruster is of particular interest. It is these ions that may be pulled back by the grid's potential and impinge on its surface. This may cause sputtering of grid material which not only degrades grid performance, but this grid material may also be harmful to the spacecraft surfaces. The base case and the cross section case both show a ux of ions back into the grid. Thermal velocities are not strong enough to overcome the electric eld. The magnitude of the di erence in these two cases corresponds to the factor of three di erence in charge-exchange reaction rate. The momentum transfer case has less ions being pulled back, because a substantial number of reactions lead to charge exchange ions with large axially directed velocities. Away from the thruster exit near the axis, the pro les are nearly identical for these cases as well as for the Te = 5 eV
85
CEX Ion FluxЧ10-18 (m-2 s-1)
1 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 0
Base CEX Momentum
0.025
0.05
0.075
0.1
Radial Distance (m)
Figure 5.16. Comparisons of charge-exchange ion ux at the thruster exit for various uniform velocity cases.
case. Figure 5.17 shows a radial pro le of charge-exchange ion ux 20 cm from the thruster exit plane. The ux in the radial direction gives an indication of the depletion of charge-exchange ions from the ion beam. The simulation which uses the larger charge-exchange cross section, labeled CEX, gives a ux higher than the base case, because the ion density is about three times higher. This is apparent from Fig. 5.18 which shows the radial velocity. The higher temperature case has a higher radial velocity than the base case, but the ux is nearly identical. The momentum transfer case shows how the postcollisional velocity of these chargeexchange ions is distributed isotropically. Its ux is comparable to the base case, but the average velocity is considerably higher. The velocity distribution functions of charge-exchange ions from the simula-
86
CEX Ion FluxЧ10-17 (m-2 s-1)
10
9
Base
CEX
8
Momentum
Te = 5eV 7
6
5
4
3
2
1
0
0
0.2
0.4
0.6
Radial Distance (m)
Figure 5.17. Comparisons of charge-exchange ion ux in the radial direction 20 cm from the thruster exit.
15
Base
Te = 5eV CEX
12
Momentum
CEX Radial Velocity (km/s)
9
6
3
00
0.2
0.4
0.6
Radial Distance (m)
Figure 5.18. Comparisons of charge-exchange ion radial velocity 20 cm from the thruster exit.
87 tions are presented in Figs. 5.19-5.21 for the base and momentum transfer cases. Figure 5.19 shows the distribution of axial velocity at the beam edge 10 cm axially from the thruster exit. The base case merely forms charge-exchange ions with the distributions of the neutrals (both from the exit and the background). The thermal spread is very low as indicated by the sharp spike around zero velocity. The momentum transfer case shows a more uniform distribution between 0 and 35000 m/s (the relative velocity of the collisions). Figure 5.20 shows the distribution of radial velocity at the same location. Again the momentum transfer case is more uniform. Also the base case indicates acceleration of these ions by the electric eld. The shift of the peak shows an average velocity around 1000 m/s whereas the thermal velocity is only about 250 m/s. This acceleration is more apparent in Fig. 5.21. This shows the distribution in a computational cell whose center is about 45 cm axially from the exit and 20 cm from the axis of symmetry. It shows a representative distribution away from the ion beam. The peak in the base case is again shifted by acceleration to an average of 2500 m/s and has a wider spread than on the beam edge. The momentum transfer case also shows this acceleration as particles are biased towards a higher velocity. The axial velocity distribution at this location is similar to that shown in Fig. 5.19. Of particular concern is the back ow of charged particles which may be harmful to spacecraft components, such as the solar arrays. In this computational domain (shown in Fig.4.3) only a small part of the back ow region is represented. A plot of the magnitude of the ion current density 10 cm axially from the thruster exit in the back ow region is presented in Fig. 5.22. This plot shows the variation for the di erent simulations. The cross section case is the only one substantially di erent and again shows a higher density of charge-exchange ions than the other
88
Normalized Distribution Function
100 10-1
Base Momentum
10--210
0
10
20
30
Axial Velocity (km/s)
Figure 5.19. Normalized axial velocity distribution functions of charge-exchange ions at the beam edge (r = 5cm), 10 cm away from the thruster exit plane.
Normalized Distribution Function
100 10-1
Base Momentum
10--210
-5
0
5
10 15 20
Radial Velocity (km/s)
Figure 5.20. Normalized radial velocity distribution functions of charge-exchange ions at the beam edge (r = 5cm), 10 cm away from the thruster exit plane.
89
Normalized Distribution Function
0.20 0.15
Base Momentum
0.10
0.05
0.00-10
-5
0
5
10 15 20
Radial Velocity (km/s)
Figure 5.21. Normalized radial velocity distribution functions of charge-exchange ions away from the beam (r 20cm), 45 cm away from the thruster exit plane.
simulations. These ions are not the only concern for interaction with the spacecraft. The grid material sputtered by the charge-exchange ions may also be ionized. These heavy metallic ions would then be a ected by the electric elds created by the charge-exchange ions. The next section examines how these self-consistent elds are a ected by altering the physical parameters.
5.6 Electric Field Results The plasma potential along the axis from various simulations is presented in Fig. 5.23 and compared with experimental measurements of the oating potential taken by de Boer. 11] The di erence in the plasma potential and the oating potential measured by the probe is a factor of order one times the electron temperature and depends on ion energy. 11] For an electron temperature of 1 eV and an ion kinetic energy of 1000 eV, the oating potential is about 2 V less than the plasma potential. Since the electron temperature varies in the plume, the plasma
90
0.006 0.005 0.004
Base Te = 5eV CEX Momentum
J (A/m2)
0.003
0.002
0.001
0.0000.0 0.1 0.2 0.3 0.4 0.5 0.6 Radial Distance (m) Figure 5.22. Comparisons of the magnitude of ion current density 10 cm from the thruster exit in the back ow region.
potential can be obtained from the oating potential measurements using the electron temperature value at the corresponding location. However, this conversion is not straightforward, because there is uncertainty in the value of the electron temperature in the plume. Di erent results are obtained depending upon the technique used to obtain the electron energy. 11] Thus, the simulations are compared to oating potential. If a simple correlation between oating and plasma potential is assumed, the largest discrepancy is again due to the lack of beam focusing in the simulations. The simulation results are time averaged over 1,000 time steps. In the actual simulations these values vary each time step. Given this level of agreement with the data, only estimates of the electric elds are made. A radial pro le 30 cm from the exit shown in Fig. 5.24 shows a strong radial electric eld at the edge of the beam where the gradients in charge density are most pronounced. The sharpness of the peaks indicates the width of the beam edge. The edge is well-de ned, because the beam ions are only moderately a ected by the electrostatic forces. Further away from the axis, the eld is nearly negligible. A
91
10
9.5
Base
CEX
9
Momentum
Data
8.5
Potential (V)
8
7.5
7
6.5
6
10
20
30
40
50
60
Axial Distance (cm)
Figure 5.23. Comparisons of plasma potential along axis for various cases with de Boer's data for oating potential 10]
plot of the axial electric eld shows a similar trend and is not included. It indicates a strong eld at the edge of the ion beam which pulls ions towards the back ow region. Figures 5.25 and 5.26 show the sensitivity of the potential and electric eld to model assumptions. At a xed radial position 15 cm from the thruster centerline, the potential and radial electric eld are compared for various cases. The variations in potential seen in Fig. 5.25 re ect the di erences in ion density in these cases, since log(n). The slope of the lines behind the thruster exit indicate the electric elds that accelerate the ions in the back ow region and towards the spacecraft. The magnitude of this axial electric eld ranges from 10 - 30 V/m for these cases. Further away from the axis, this eld goes to zero. The values for the potential are considerably lower than the 7-12 V found experimentally by Pollard 13], even though the ion density values are comparable. The negative potential is an artifact
92
Radial Electric Field (V/m)
700 600 500 400 300 200 100 0 -100 -200 0
Base CEX Momentum
0.1
0.2
Radial Distance (m)
Figure 5.24. Comparisons of radial electric eld 30 cm away from the thruster exit.
of the Boltzmann relation. This discrepancy suggests that away from the beam the potential may not be de ned by the Boltzmann relation. The radial electric eld shown in Fig. 5.26 has more variation than the previous eld in Fig. 5.24. The peaks are where the edge of the beam is encountered. These again are more pronounced for the momentum transfer case since the charge-exchange ion density is lowest for this simulation.
5.7 Conclusions The computer code performs adequately to address the issues involved in simulating the plumes of plasma thrusters. The expected ion dynamics were captured qualitatively. The beam ions spread more due to electric elds than they would due to purely thermal e ects. The charge-exchange ions were pulled back towards
93
5
0
Plasma Potential (V)
-5
-10
Base
CEX
Momentum
-15
0
0.5
1
Axial Distance From Thruster Exit (m)
Figure 5.25. Comparisons of potential 10 cm above the thruster.
Radial Electric Field (V/m)
250
200
Base
CEX
Momentum
150
100
50
0
-50
0
0.5
1
Axial Distance From Thruster Exit (m)
Figure 5.26. Comparisons of radial electric eld 10 cm above the thruster.
94 the thruster near the exit and away from the beam further away. Analysis of the velocity distribution functions of the charge-exchange ions indicates that the collisional phenomena are captured properly as well. Results from simulations compared satisfactorily with available experimental data for ion ux, ion density, and oating potential in the plume of the UK-10 ion thruster. These comparisons also indicated the signi cance of the initial ion velocity pro le assumed in the computations. Computations where the acceleration grid was assumed to be a section of a sphere indicated that the ion beam, although altered noticeably, followed its initial pro le in the near plume. The various uniform velocity simulations showed very little variation in this region. A better understanding of the trajectory of the ions between the screen grids is needed to determine their pro le at the thruster exit more accurately. Because the focusing of the grids is not as pronounced as the simulations predict, various ion inlet pro les with imposed divergence angles could be tested as is done for the Hall thruster simulations presented in the next chapter. An improved pro le is necessary to obtain better agreement with the experimental data. The pro le of the charge-exchange ions was found to be sensitive to various physical assumptions. Among these were electron temperature, background pressure, charge-exchange collision cross section, and the collision dynamics for chargeexchange reactions. Comparisons of ion density with experimental data indicated the signi cance of the charge-exchange cross section, especially away from the ion beam. These comparisons also demonstrated that the change in momentum after charge-exchange collisions was less pronounced than that predicted by hard sphere collisions. Computations that employed either hard sphere collision dynamics for charge-exchange collisions or an increased electron temperature underpredicted
95 the ion density in comparison with experimental measurements. Moreover, the charge-exchange ion pro les in the near eld and back ow regions depend on these assumptions. The sensitivity of the charge-exchange ions to these assumptions indicates a need for better understanding of these phenomena. Along the axis, the various computations gave potential values in reasonable agreement with experimental data. However, the potential in the plume away from the beam was underpredicted by these computations. This suggests that the Boltzmann relation does not accurately de ne the potential outside the beam.
Chapter 6 Hall Thruster Simulations Simulations are performed in order to model the plumes of Hall thrusters. Using available experimental data for validation, the ability to simulate these plumes accurately is assessed. The plume of the stationary plasma thruster SPT-100 is modeled, because there is a wealth of experimental data available. 15]{ 22] Measurements of current density and ion velocity were reported in Refs. 15], 16] and 19]. Electron density measurements are reported in Ref. 21]. More recently, measurements of ion velocity, ion density, heat ux, neutral ux and plasma potential were obtained. 17,18,22] The present study relies on these data to perform the most comprehensive comparisons between experimental and computational work thus far. Two separate studies are presented | one which addresses the sensitivity of the plume pro le to various inlet pro les, and another which compares results from the two models discussed in Chapter 3. These models include the simple Boltzmann relation and the variable electron temperature model for determining the electric eld. Along with the sensitivity study, simulations are performed to illustrate the advantage of increasing particle count by using a parallel computer. 96
97 6.1 Inlet Pro le
The nominal operating conditions given in the experiments are a discharge current
of 4.5 A, a discharge voltage of 300 V and a total ow rate of 5.2-5.4 mg/s with
7% cathode split. 18] The thrust produced is 85 mN. These parameters de ne the
thruster exit plane conditions, which are used for input values to the simulations.
However, it is necessary to separate these into densities and velocities for each
species.
Using the ow rate and the thrust, the densities and velocities can be obtained
by setting either the current or the thruster e ciency. The anode ow rate is:
;(1
%c)m_ = m2
Z rout rin
spX eciesnj j
(r)vj
(r)rdr
(6:1)
where %c is the cathode split. The thrust is:
Ft
= m2
Z rout rin
spX eciesnj j
(r)vj2(r)rdr:
(6:2)
In these equations, each of the densities and velocities are functions of radial posi-
tion, and the velocities are in the axial direction. These functions can be written
as:
nj(r) = njmaxgj(r) vj(r) = jvjjhj(r)
The fraction of doubly-charged ions in the plume is de ned as:
=
ni
nii + nii
:
These ions would be accelerated more by the electric eld than singly-charged ions.
Assuming a xed constant, this relates the two velocities by:
vii = avi:
98
As with the ion thruster, sonic conditions are assumed at the thruster exit
for the neutrals. They are also assumed to be uniform across the exit, because
the electric eld does not a ect their pro le. The temperature is assumed to be
1000 K, because the thruster walls become heated when the thruster is ring. The
wall temperature is near 800 K. 15] Using 1000 K gives a velocity of 325 m/s for
xenon (va) which agrees well with experiments using a similar Hall thruster. 67]
The ow through the cathode is 7% of the total ow rate. This ow is also assumed
to be entirely neutral xenon at the thruster exit in order to match the total ow
rate. These xenon atoms account for a substantial fraction of the neutrals entering
the domain. As with the ion thruster simulations, the background pressure of
the experimental facilities must be included in the simulations. The background is
assumed to be composed entirely of neutral xenon atoms at 300 K. The background
pressure of the facility in which the bulk of the experimental data were obtained
is 6 mP a. 18]
Equation 6.1 can be written as:
; j j ; (1
%c)m_ = m2
Z rout rin
" nava
+
nimax
vi
f1(r) + 1a
!# f1(r) rdr
(6:3)
where f1(r) = gi(r)hi(r). Equation 6.2 can similarly be expressed as:
j j ; Ft = m2
Z rout rin
" nava2
+
nimax
vi
2
f2(r)
+
a2 1
!# f2(r) rdr
(6:4)
where f2(r) = gi(r)h2i (r). De ning A na, B nimaxjvij, and C jvij, it is clear
that another equation relating the three variables would yield these unknowns. As
with the ion thruster, once these values are known the functions for the densities
and velocities are discretized for a xed number of cells at the thruster exit. These
simulations use eleven cells of equal radial spacing.
99
6.1.1 Variations in Ion Pro le
This study uses e ciency as the third constraint. The ion current at the thruster
exit is then variable. The propulsion e ciency is approximately 0.5-0.55 for the
SPT-100 Hall thruster at these operating conditions. 20] This e ciency is de ned
as a function of ow rate and thrust, which are already used as constraints. In
order to close this system of equations, an e ciency is de ned as:
=
0:5m_ u~2 IV
where I and V are the discharge current and voltage, and u~2 is de ned such that:
; 2 IV = m2
Z rout rin
" Ava3
+
BC
2
f3(r)
+
a3 1
!# f3(r) rdr
(6:5)
where A, B, and C are as de ned previously, and f3(r) = gi(r)h3i (r).
The analysis is simpli ed for these simulations. First, the ion densities are
assumed to be radially uniform. Thus, the current density varies with vi using a
divergence angle. Measurements by Manzella 15] indicate a radial component to
the ion velocity in the near eld. Oh and Hastings inverted an analytic function
to choose a divergence angle for each particle. 25] The present study simpli es this
process by assuming the divergence angle varies linearly with radial distance from
the center of the annulus ring. The ion velocity varies as:
vi(r) = jvij cos k(r ; r)]
where k = rout;rin , and r is the radial position halfway between the inner and outer radius. A half-angle of 10 degrees is used in one case and 20 degrees in another. A uniform ion velocity (0 m/s radial velocity) case is also considered for comparison. The fraction of doubly-charged ions is estimated as 25% by Manzella 16] for an SPT-100. This value is used for this study. Assigning the same velocity to the
100
Table 6.1. Values at the thruster exit for Case 1.
Species Xe Xe+ Xe++
Density Velocity Temperature
1017m;3] km/s]
eV]
1.58
0.325
0.086
1.74
17.6
4
0.58
17.6
4
doubly-charged ions as is used for the singly-charged ions leads to an ion velocity (jvij) that closely matches average ion velocity measurements in the plume obtained by King. 17] Plots of neutral density versus e ciency at a thrust of 85 mN and a ow rate of 5.2 mg/s indicate that the e ciency must be about 0.54 or greater in order for the neutral ux to be positive. A value of 0.54 is used. Table 6.1 shows the densities and velocities for neutrals, singly-charged ions, and doubly-charged ions. These values give a current of 4.4 A and a utilization e ciency of 0.988. In this study, the electron temperature is assumed to be 3 eV based on Ref. 21] and is used to obtain the potential from the Boltzmann relation. Based on plume measurements of plasma potential, 17] the exit plane potential is taken as 10 V. This sets the value of ne1 to 1:0 1016m;3. An ion temperature of 4 eV is assumed at the exit based on Ref. 15] for most of the simulations, and is compared with a 0.4 eV case (labeled Case 3 in the gures). The various simulations performed for this study are presented in Table 6.2.
101
Table 6.2. Ion conditions at thruster exit for sensitivity study.
Case Speed Velocity Temperature
km/s] Variation eV]
1 17.6 g(r) : 10
4
2 17.8 g(r) : 20
4
3 17.6 g(r) : 10
0.4
4 17.5 uniform
4
6.1.2 Ion Pro le for Electron Model Comparisons
This study attempts to improve on the thruster exit plane pro le described in
the previous section. It uses the current to obtain the densities and velocities. The maximum ux (B = nimaxjvij) can be obtained directly from the current by integrating the normalized ux f1(r). Based on near eld current density mea- surements, 19] a current of 4 A is assumed at the exit plane. This data also shows
a more pronounced peak than the uniform ion density model gives. Therefore, a
Gaussian pro le is imposed of the following form:
ni(r)
=
nimax
exp
" ;
(r
; &
r)2
#
and two values of & are used which a ect the width of the current density. The
fraction of on plume
doubly-charged ions ( ) is assumed to be only 10% in measurements by King. 18] Their velocity is assumed
tthoisbsetupdy2
based times
that of the singly-charged ions based on the ratio of charges. The ion velocity and
neutral density can then be obtained from the thrust and ow rate. The thrust
is set to 84.9 mN, and the ow rate used is 5.4 mg/s. Based on results from the
102 Table 6.3. Values at the thruster exit for reference (Ref) case.
Species Xe Xe+ Xe++
Density Velocity Temperature
1017m;3] km/s]
eV]
5.76
0.325
0.086
2.54
16.8
4
0.28
23.8
4
sensitivity study, the half-angle of divergence is assumed to be 10 degrees. Densities and velocities for & = 2 10;4 are presented in Table 6.3. These values give a utilization e ciency of 0.956. The normalized current density at the thruster exit for these assumptions is presented in Fig. 6.1. An ion temperature of 4 eV is again assumed. Table 6.4 shows the various simulations performed to examine the e ects of a varying electron temperature. The value for reference electron number density determines the magnitude of the potential, but has little e ect on the di erences in the potential. The listed value for the electron temperature in the variable electron temperature cases represents the far eld value. A value of 1.0 eV is based on Kim's measurements, 19] whereas 2.0 eV is based on measurements by Myers and Manzella. 21] The imposed pro le discussed in Chapter 3 is used for the electron temperature in the plume. The cases marked with asterisks include a full chamber geometry simulation as well as a smaller domain one. The full chamber simulations allow comparisons with current density measurements at angles up to 180 degrees from the thruster centerline. These also allow examination of any possible chamber e ects.
103
1
Normalized Current Density
0.9
0.8
0.7
0.6
0.5 0.03
0.035
0.04
0.045
0.05
Radius (m)
Figure 6.1. Normalized current density at the thruster exit used as input for simulations.
Table 6.4. Input conditions for modeling study. Case Te(eV ) Te pro le ne1(m;3) Ref* 3.0 constant 6.5e15 VT 1 1.0 varying 6.5e15 VT 2 2.0 varying 6.5e15 VT 3* 2.0 varying 1.3e15
104 6.2 Results from Sensitivity Study In this study, the ion pro le at the thruster exit plane is varied to examine the sensitivity of the simulations to this pro le. Measurements of radial velocity by Manzella 15] indicate an almost linear variation with radial position 11 mm from the thruster exit. A comparison of the simulations with this data is shown in Fig. 6.2. The simulations also give an almost linear variation that is reasonably close to the data at this axial location. There are noticeable discrepancies between simulation and data. Neither case reproduces the asymmetry found experimentally, because it was not included in the thruster exit plane pro le for this investigation. The di erence in the simulations is the imposed divergence angle of the ions. Case 1 assumes a maximum divergence angle of 10 degrees, whereas Case 2 assumes 20 degrees.
8
Ion Radial Velocity (km/s)
6 Case 1
Case 2
4
Data
2
0
-2
-4
-6
-0.005
0
0.005
0.01
Radial Position (m)
Figure 6.2. Comparisons of ion radial velocity with Manzella's measurements 14] at an axial distance of 11 mm.
105 6.2.1 Current Density Comparisons for current density from various simulations with experiments are presented at two di erent radial locations in Figs. 6.3{6.5. The base case (Case 1) shown in Fig. 6.3 agrees well with King's experimental data. 18] The comparisons in Figs. 6.4 and 6.5 indicate the insensitivity of the computed current density to the variations in inlet ion pro les. The variations are more pronounced close to the centerline however, each case is within a factor of about 2 of the data. Comparing Cases 1, 2, and 4 indicates the signi cance of the divergence angle. Clearly, the higher angle leads to a lower peak value at these radial distances. Experiments show a more de ned peak than the simulations. The low ion temperature case (Case 3) gives a more de ned peak, but an underprediction of current density occurs at higher angles.
101
Case1 50cm Case1 1m Data 50cm Data 1m
100
Current Density (mA/cm2)
10-1
10-2 -50
0
50
Angle (deg)
Figure 6.3. Comparisons of current density with King's data 17] at radial distances of 50 cm and 1 m.
106
101 100
Case1 Case 2 Case 3 Case 4 Data
Cuurent Denisty (mA/cm2)
10-1
10-2 -50
0
50
Angle (deg)
Figure 6.4. Comparisons of current density with King's data 17] at a radial distance of 50 cm.
101 100
Case1 Case 2 Case 3 Case 4 Data
Current Density (mA/cm2)
10-1
10-2 -50
0
50
Angle (deg)
Figure 6.5. Comparisons of current density with King's data 17] at a radial distance of 1 m.
107 6.2.2 Ion Velocity and Density In Figs. 6.6{6.8 the current density is separated into density and velocity components. The velocity comparisons in Figs. 6.6 and 6.7 indicate agreement within experimental uncertainty up to about 50 degrees. The simulation results only present the average ion velocity of the propellent ions, neglecting the charge-exchange ions. The experimental data were obtained using a retarding potential analyzer (RPA). 18] Numerical di erentiation of the current versus voltage data approximates the derivative of this curve which is proportional to the energy-per-charge distribution function. For singly-charged ions of uniform mass, the energy distribution function also represents the velocity distribution function. By normalizing and numerically integrating this velocity distribution function, the average velocity is obtained. This process is likely to minimize the e ects of the charge-exchange ions, because they only contribute to the current measurements at low voltages. Thus, the actual contribution of the charge-exchange ions to the average velocity may be obscured. Therefore, comparing the average ion velocity of only the propellent ions from the simulation with the experimental measurements is reasonable. The agreement is encouraging. The ion velocity and current density are used to obtain the experimental number density shown in Fig. 6.8. The simulations allow the ion density to be given directly. Both total density and the propellent ion density of the simulations are shown in Fig. 6.8. The propellent ion density agrees well at low angles, but it is too low at higher angles where the charge-exchange ions represent a higher fraction. The pro le of the total density for Case 1 compares favorably with the experiment, but the magnitude is too high. Following the same reasoning as described earlier, the measured ion density should be between the total and propellent densities from the simulation. In contrast to the ion velocity
108
25
20
Ion Velocity (km/s)
15
10
5
Case 1
Case 2
Data
00
10
20
30
40
50
60
Angle (deg)
Figure 6.6. Comparisons of average ion velocity with King's data 16] at a radial distance of 50 cm.
comparisons in Figs. 6.6 and 6.7, Case 2 does not represent the ion density as well as Case 1. The total ion density from the other two cases, not presented, is also higher than the measured values.
6.2.3 Heat Flux Heat ux measurements were obtained by a heat ux probe consisting of two Schmidt-Boelter thermophile transducers. 18] The probe was biased to ground potential, which was between 5 and 8 volts below the plasma potential, to eliminate electron current heating. Thus, the measured heat ux includes both neutrals and ions. The simulations give the contributions of the ions and the propellent neutrals. The background neutrals contribution is negligible. It is estimated to be on the order of 1 W=m2. In Figs. 6.9{6.11 comparisons of these simulations with heat ux measurements are shown. There is no data at angles between -2 degrees to 2 degrees, because the measured heat ux exceeded the calibrated range of the
109
25
20
Ion Velocity (km/s)
15
10
5
Case 1
Case 2
Data
00
10
20
30
40
50
60
Angle (deg)
Figure 6.7. Comparisons of average ion velocity with King's data 16] at a radial distance of 1 m.
1017
Ion Density (m-3)
1016
1015
1
0
14 0
Case 1 (Total) Case 1 (Prop.) Case 2 (Total) Case 2 (Prop.) Data
10
20
30
40
50
60
Angle (deg)
Figure 6.8. Comparisons of ion density at a radial distance of 50 cm with King's data. 16]
110 sensor. 18] In the simulations, the heat ux is obtained from particles crossing the boundary. The values obtained are therefore assigned to the center of the edge. Thus, no value is available at 0 degrees, because no edge is centered around this location. The simulations match the data well near the thruster centerline at 50 cm. At 1 m the simulated heat ux is within a factor of three of the measurements. At large angles from the centerline, where charge-exchange ions dominate, there is an underprediction by two orders of magnitude at both locations. The chargeexchange ions have much lower velocities than the propellent ions, and the heat ux scales with velocity cubed. The spreading of the propellent ions does not improve the agreement at the higher angles signi cantly, and it leads to poorer agreement near the centerline (see Figs. 6.9 and 6.10). As with the ion density comparisons, a divergence angle ranging from -20 degrees to 20 degrees (Case 2) does not lead to as good agreement as the other cases. The ion temperature has a surprisingly small e ect on the heat ux (see Fig. 6.9). The ion temperature a ects the spreading of the beam without altering the magnitude of the mean velocity. The heat ux at higher angles should increase with ion temperature, because a higher temperature allows more of the propellent ions to diverge from the centerline. To estimate the e ect on the ions of biasing the probe, an energy of 8 eV is added per ion to the heat ux calculations and is labeled 'qtot' in Fig. 6.11. Clearly, it is signi cant only at the larger angles, but it does not make a qualitative di erence. Upon integration of the heat ux from the experiments, this study nds that this heat transfer rate is greater than the power put into the system. Integrated values for the heat ux are shown in Table 6.5. The input electrical power is only 1350 W. The simulations give values indicative of the measured e ciencies ( 50%) at a radial distance of 50 cm. The energy of the ion ow at the thruster
111
105 104
Case 1 Case 2 Case 3 Case 4 Data
103
Heat Flux (W/m2)
102
101
100 -50
0
50
Angle (deg)
Figure 6.9. Comparisons of heat ux with King's data 16] at a radial distance of 50 cm.
104 Case 1 Case 2 Case 4 Data 103
Heat Flux (W/m2)
102
101
100 -50
0
50
Angle (deg)
Figure 6.10. Comparisons of heat ux with King's data 16] at a radial distance of 1 m.
112
105 104
qKE qtot Data
Heat Flux (W/m2)
103
102
101
100 -50
0
50
Angle (deg)
Figure 6.11. Comparisons of heat ux with King's data 16] at a radial distance of 50 cm.
exit is a function of the e ciency and the input power. By 1 m, collisions with the background neutrals have decreased the energy of the ion ow. The simulations which include a varying electron temperature model are also included in the table, but plots of the heat ux show only moderate di erences from the reference case. A comparison of various simulations from the modeling study (see Table 6.4) is presented in Fig. 6.12. The integrated values are slightly higher for the variable electron temperature cases, because the larger electric eld in the near eld region gives an added acceleration to the ions. Given the input power, the experimental data appears to be corrupted. As suggested by the measurements of the ion energy distribution function, 18] presented along with computations in the next chapter, the experimental facility is likely a ecting the plume data.
113
Table 6.5. Integrated heat ux at both radial locations.
Case 50 cm (W) 1 m (W)
King 1692
2339
Case 1 630.0
449.4
Case 2 637.8
460.0
Case 3 609.6
410.1
Case 4 618.3
424.8
Ref
652.5
454.9
VT 1 678.0
468.6
VT 2 678.5
472.6
FG VT 680.3
478.3
Ref
104
VT 1 VT 2
Data
103
Heat Flux (W/m2)
102
101
100 -50
0
50
Angle (deg)
Figure 6.12. Comparisons of heat ux with King's data 16] at a radial distance of 50 cm.
114 6.3 Results from Modeling Study In this study the two models for determining the electric eld are compared. These include the Boltzmann relation and the variable electron temperature model with the imposed electron temperature described in Section 3.2. The temperature in the computational domain is from an analytic expression and is based on near eld measurements. 19] 6.3.1 Electron Density Calculation of electron density from both studies are compared with measurements by Myers and Manzella 21] at a radial distance of 31 cm. Values from the sensitivity study are shown in Fig. 6.13. Cases 1 and 2 show the best agreement. The di erence in the two simulations is that a 20 degree divergence angle is assumed for Case 2. The other simulations have inlet conditions which lead to less spreading of the ions, namely a lower ion temperature (Case 3) and no imposed divergence angle (Case 4). They overpredict the density near the axis. The electron densities from the simulations from the modeling study are shown in Fig. 6.14. Recall the reference case (Ref) uses the Boltzmann relation while the VT simulations use the variable electron temperature model. Including a varying electron temperature a ects the spreading of the ions. The radial electric eld is larger than it would be without a temperature gradient. The case labeled 'Ref' is nearly identical to Case 1. It is clear from Fig. 6.14 that the case labeled 'VT 3', which assumes a far eld electron temperature of 2 eV, agrees best with the measurements. The other varying electron temperature simulation (VT 1) shown in Fig. 6.14 matches the data well at large angles, but gives a higher number density on axis. The di erences
1017 1016
115
Case 1 Case 2 Case 3 Case 4 Data
ne (m-3)
101-550
-25
0
25
50
Angle (deg)
Figure 6.13. Comparisons of electron density from sensitivity study with measurements by Myers 20] at a radial distance of 31 cm.
in the simulations is caused by the far eld electron temperature. The magnitude of the electric eld is proportional to the electron temperature using these models. Therefore, the lower electron temperature (VT 1) leads to less acceleration and a higher plasma density. The comparisons with the data indicate the importance of the far eld value. An advantage of the varying electron temperature is that a smaller temperature can be used in the far eld than in the near eld. The comparisons also suggest that 2.0 eV is a better far eld value than 1.0 eV or 3.0 eV.
6.3.2 Near Field Current Density Comparisons with measurements of current density in the plume near eld by Kim 19] are shown in Figs. 6.15 and 6.16 at three axial locations. At an axial
116 1017
ne (m-3)
1016
Data
Ref
VT 1
101-550
VT 3
-25
0
25
50
Angle (deg)
Figure 6.14. Comparisons of electron density from modeling study with measurements by Myers 20] at a radial distance of 31 cm.
distance of 10 mm, the agreement is reasonable regardless of the model used. However, by 100 mm the e ects of an electron temperature gradient are more apparent. The ion beam reaches the axis more quickly than is predicted by the reference case, because the electric eld attributed to the temperature gradient is in the same direction to that due to the density gradient in these models. The variable electron temperature case more closely matches the shape of the measurements at 100 mm than the reference case. Its peak is about 50% higher than that of the reference case as well. Even by 50 mm the current density is higher at the axis for a varying electron temperature. For these near eld comparisons, there is negligible di erence between the di erent variable electron temperature cases, therefore only the rst is shown.
117
Current Density (mA/cm2)
80
DSMC/PIC
10 mm
70
50 mm 100 mm
60
Data
10 mm
50
50 mm
100 mm
40
30
20
10
0
0
0.05
0.1
0.15
0.2
Radial Distance (m)
Figure 6.15. Comparisons of current density from reference case with measurements by Kim. 18]
Current Density (mA/cm2)
80
DSMC/PIC
10 mm
70
50 mm 100 mm
60
Data
10 mm
50
50 mm
100 mm
40
30
20
10
0
0
0.05
0.1
0.15
0.2
Radial Distance (m)
Figure 6.16. Comparisons of current density from variable electron temperature case with measurements by Kim. 18]
118 Table 6.6. Integrated current (in Amps) at various axial locations. Case 10 mm 25 mm 50 mm 100 mm Kim 3.97 4.92 4.95 4.51 Ref 4.04 { 4.09 4.19 VT 1 4.14 { 4.24 4.33 VT 3 4.13 { 4.22 4.32 Integrated values are compared in Table 6.6. At 10 mm and 100 mm, varying the electron temperature gives only a 4% di erence in integrated current. A simulation which assumed & = 1 10;4 gives a pro le that is much more pronounced than the reference case. It more poorly represented the experimental data. The uniform ion density pro le used in the sensitivity study gives a wider peak at 10 mm, as shown in Fig. 6.17. These comparisons indicate the importance of including the variation in electron temperature in simulations. The spreading of the ion beam is magni ed by this variation. The current density measurements show more focusing than the simulations, but the narrow peak found experimentally nearest the exit indicates that the half-angle of divergence at the thruster exit does not exceed 20 degrees. It appears that some other e ect is responsible for the ion beam focusing at the axis. One possibility is the e ect of the magnetic eld on the electrons in the near plume.
119
Current Density (mA/cm2)
80
70
Base 10mm
Data 10mm 60
50
40
30
20
10
0
0
0.05
0.1
0.15
0.2
Radial Distance (m)
Figure 6.17. Comparisons of current density from Case 1 with measurements by Kim. 18] 6.3.3 CEX Ion Pro le
For contamination concerns, the back ow and far eld regions of the plume are more signi cant than the near eld. The next few gures compare the reference case with the variable electron temperature simulations. Figure 6.18 shows the densities of the various ions 4 cm axially from the thruster in the back ow region. Clearly, the charge-exchange ions are dominant. This is also true well above the thruster in the radial direction, because the propellent ions are predominantly axially directed. Comparisons of the charge-exchange ions both in the back ow region and above the thruster for the various cases indicate the signi cance of the magnitude of the far eld electron temperature. More importantly, the current density shows di erences at these locations as seen in Figs. 6.19 and 6.20. At the locations in the back ow region (Fig. 6.19), the reference case has a magnitude
120
Ion Density (m-3)
1015 1014 1013 1012
Xe+cex Xe++cex Ref VT 1 VT 2 Xe+p
1011
Xe++p
1010
0.2
0.4
0.6
0.8
1
Radial Distance (m)
Figure 6.18. Comparisons of ion density for the various simulations at an axial distance 4 cm from the thruster exit in the back ow region.
that is a factor of two higher than the Te = 1eV case near the thruster and 50% higher near the top of the domain. The negative values indicate the ion ow is in the opposite direction of the main ion ow from the thruster. The results are similar at 70 cm above the thruster (Fig.6.20) as well in both the axial and radial directions.
6.3.4 Full Chamber Geometry Simulations which include the full vacuum chamber geometry allow comparison with experimental data at large angles in the back ow region. The back ow region includes the half of the domain in the direction opposite to which the thruster is ring. This extends from 90 degrees o -axis to 270 degrees. The current density results from the full geometry cases are compared with measurements by King 18]
121
0
Axial Current Density (A/m2)
-0.1
-0.2
Ref VT 1 VT 3
-0.3
-0.4
-0.5
0.2
0.4
0.6
0.8
1
Radial Distance (m)
Figure 6.19. Comparisons of axial current density for the various simulations at an axial distance 4 cm from the thruster exit in the back ow region.
Current Density (A/m2)
0.4 0.3 0.2 0.1 0 -0.1
jr Ref VT 1 VT 2 jx
0
0.5
1
Axial Distance (m)
Figure 6.20. Comparisons of current density for the various simulations at a radial distance 0.7 m above the thruster.
122
Current Density (mA/cm2)
Ref 50cm
Ref 1m
100
Data 50cm
Data 1m
VT 50cm VT 1m
10-1
10-2
10-3 -150 -100
-50
0
50
Angle (deg)
100 150
Figure 6.21. Comparisons of current density with measurements by King. 17] at a radius of 50 cm and 1 m in Fig. 6.21. Even at large angles into this back ow region, good agreement is obtained for both models. These simulations match the data as well as, if not better than, the simulations considered in the sensitivity study. The current density is found to be relatively insensitive to the varying electron temperature. Di erences in electron number density are found for the various simulations. These ndings suggests that the far eld electron temperature is important. The variations of the electron temperature in the near eld do not a ect the far eld current density. However, the densities are dependent on the value of the far eld electron temperature even for simulations with the same near eld electron temperature pro le. The value of the electron temperature a ects the velocity and thus the density of the charged particles. Comparison of the full geometry simulations with their smaller domain counterparts shows only small di erences. For the locations considered in the previous
123
0
-0.1
Axial Current Density (A/m2)
-0.2
-0.3
-0.4
-0.5
0.2
Ref VT 3 Ref FG VT FG
0.4
0.6
0.8
1
Radial Distance (m)
Figure 6.22. Comparisons of axial current density 4 cm from the thruster exit (in the back ow region) for the full geometry cases and the smaller domain cases.
gures, the plasma potential, electron density, and the current density are nearly identical. However, in the back ow region, di erences are noticeable but small (see Fig. 6.22). Again the negative values indicate the ow direction. Only small di erences between the full geometry cases are found in this region for current density. More signi cant di erences are seen about 2.5 m above the thruster in Fig. 6.23. A simulation has also been performed in which the facility background pressure is neglected. Current density results from this simulation are presented with the reference case values in Fig. 6.24. Clearly, the facility back pressure has a signi cant e ect on the measurements in the plume. These background particles collide with the propellent ions. Nearly all of the current density in the back ow region is from the charge-exchange ions created by these collisions.
124
Current Density (A/m2)
0.02 0.015 0.01 0.005 0 -0.005
jy jx
Ref FG VT FG
0
2
4
6
8
Axial Distance (m)
Figure 6.23. Comparisons of current density 2.425 m above the thruster for the full geometry cases.
101
Ref 50cm
Ref 1m
100
No pb 50cm
No pb 1m
10-1
Current Density (mA/cm2)
10-2
10-3
10-4 -150 -100
-50
0
50
Angle (deg)
100 150
Figure 6.24. Comparisons of current density with and without back pressure.
125 6.3.5 Potential The plasma potential is compared with measurements taken by Marrese 22] at an axial distance of 48 cm in Fig. 6.25. In the simulations the potential drops below zero due to the choice of the reference electron density used in inverting the Boltzmann relation. The experimental measurements here and by King 18] (in Fig. 6.26) do not become negative. The third variable electron temperature simulation (VT 3) is shown which uses a lower value for the reference number density based on Marrese's measurements. This value also agrees with the far eld electron density of the simulations. This simulation matches Marrese's data well. The rst two variable electron temperature cases have the same shape as Marrese's data, but have negative values like the reference case. The agreement with King's data shown in Fig. 6.26 is less encouraging. The shape of the potential pro le measured experimentally is explained in Ref. 28] through a detailed analysis of the e ect on plasma potential of the magnetic con guration of various Hall thrusters. 6.4 Parallel Study The purpose of this study is to illustrate the importance of increasing the particle count by using more processors. This study uses input pro les described by the sensitivity study. This base case is run on one, four, and eight processors. The total number of particles varies from 350,000 to 3,000,000. No back pressure is used in these cases so that the charge-exchange collisions are more rare, thus allowing the statistical behavior to be better examined. Figures 6.27{6.29 show results from these simulations. The aggregate ow eld values at the positions considered in the previous gures are not noticeably changed
126
Potential (V)
10
8
6
4
2
0 Data
Ref
-2
VT 1
VT 2 VT 3 -4
0
20
40
R (cm)
Figure 6.25. Comparisons of potential with measurements by Marrese 21] at an axial distance of 48 cm.
15 VT 3 Data 12
Potential (V)
9
6
3
00
10 20 30 40 50 60
Angle (deg)
Figure 6.26. Comparisons of plasma potential with measurements by King 17] at a radial distance of 50 cm.
127 by the number of processors used. The PIC values are more a ected, because the cells are smaller than the DSMC cells and the ions are the only particles contributing to the charge density at the nodes of the PIC grid. Therefore, the e ects on the electric eld are examined. In Fig. 6.27, instantaneous values of potential for the one processor and eight processor cases are compared with a time-averaged value. Both are nearly identical to the average near the thruster centerline, but the one processor simulation exhibits more noise. These uctuations are even more pronounced at the large angles, where charge-exchange ions dominate. The instantaneous potential for the eight processor case has only minor uctuations even at 1 m where the ion density has dropped signi cantly. Figure 6.28 shows uctuations of the axial electric eld for all three cases. However, the eight processor case maximum value never exceeds 50 V/m. This would only change the velocity of a singly-charged ion by about 1 m/s in a time step. Even for the four processor case the noise is noticeable. Similar trends are seen in the radial electric eld. To emphasize the advantage of a large particle count, Fig. 6.29 shows a comparison of the current density of the charge-exchange doubly-charged ions 10 cm before the thruster exit plane. In the back ow region resolution is generally di cult. For this reason, the most minor species is considered. None of the cases are particularly noisy, but the sharp gradients in the one and four processor cases are indicative of noise due to an insu cient number of particles in the cells. The wall is assumed to neutralize the ions, therefore the eight processor pro le is as expected. This study illustrates the importance of including su cient particles to simulate the plumes. This is not only necessary to represent the properties of species in the back ow region but is also important for the charge density in the PIC algorithm. An increase in particle count is made possible by parallel simulations which use
128
5
8 Procs
Average
0
1 Proc
Plasma Potential (V)
-5
-10
-15
-20
-25 0
10
20
30
40
50
60
Angle (deg)
Figure 6.27. Comparisons of instantaneous potential values at a radial distance of 50 cm.
Axial Electric Field (V/m)
200
150
8 Procs 4 Procs
1 Proc
100
50
0
-50
-100 0
10
20
30
40
50
60
Angle (deg)
Figure 6.28. Comparisons of instantaneous axial electric eld values at a radial distance of 1 m.
129
Current Density (A/m2)
2.0E-2 1.5E-2
1 Proc 4 Procs 8 Procs
1.0E-2
5.0E-3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Radial Distance (m) Figure 6.29. Comparisons of charge-exchange Xe++ current density in the back ow region. the background routine and the weight by species technique. This increase better represents the plume and allows various physical models to be examined.
6.5 Conclusions Simulations were performed for the SPT-100 at nominal operating conditions and compared with available experimental data. These data include current density in the near eld, far eld, and back ow regions, ion velocity, ion density, electron density, heat ux, and plasma potential. Based on comparisons with these data, the simulations represented the plume of this thruster reasonably well. Calculations of far eld ion properties are a ected by variations in the ion inlet pro le. However, the comparisons with experimental data show only moderate di erences. Each of the various simulations agreed well with experimental measurements of far eld ion current density. Previous computational work compared results only to current density measurements. 29] The good agreement found in
130 these studies is not surprising. The present study was able to examine the plume in more detail because of the availability of additional properties. 18] Including a varying electron temperature was found to be important in the modeling of Hall thruster plumes. The electron temperature gradient signi cantly a ects the near eld of the plume. Simulations incorporating the varying electron temperature model gave better agreement with the near eld current density measurements than those using the isothermal Boltzmann relation. Results suggest that outside of this region the plasma behavior is not dictated by this near eld behavior, but rather by the value of the electron temperature. Thus, the far eld electron temperature is also important. This temperature is likely to be lower than the near eld value. The simulations also agreed well with the measurements of electron number density. Agreement with the plasma potential was dependent upon the value for the reference density used in the physical models. However, it is the di erences in the potential in the plume that are important for the ion pro le. In order to use measurements in the back ow region (from 90 degrees o -axis to the opposite side of the thruster exit) for code validation, it is important to simulate the full geometry of the experimental facility. These simulations can also aid in determining any possible e ects from the facility. The simulations were found to represent the current density very well. Even at large angles in the back ow region good agreement with experiment was obtained.
Chapter 7 Analysis of Collision Phenomena The results of the previous chapter suggest that the code has been validated through comprehensive comparisons with available experimental data. The data include current density in the near and far eld, ion density and velocity, electron density, and plasma potential. 15,17,18,19,21,22] Although the level of agreement suggests that the model does a reasonable job of representing the plume of this thruster, the simulations described in this chapter examine this conclusion in more detail. Available data of the ion energy distribution allows a more in depth assessment. Experimental work at the University of Michigan by King provides not only the current density in the plume but also the microscopic distribution of ion current as a function of energy. 18] This data can provide information about the collisional behavior in the plume. The study presented in this chapter examines in detail the collisional behavior in the plume of the SPT-100. The ion energy distribution is sampled in the numerical simulations to examine the e ects of ion-neutral and ion-ion collision types. The various collisions have been discussed in detail in Chapter 3. 131
132 7.1 Inlet Pro le Flow conditions for these simulations are again based on nominal operating conditions of the SPT-100. Using a ow rate of 5.2 mg/s (with 7% cathode split) and a thrust of 84.9 mN, the densities and velocities for each species can be obtained from Equations 6.1 and 6.2 once a third equation is assumed. Unlike the previous studies, the ion kinetic energy is used as the third constraint. The peak of the ion energy distribution functions occurs near 250 V for most of the experimental data taken by King. 18] The physics in the acceleration chamber of the Hall thruster allows the ions to be accelerated to nearly the applied potential of 300 V as described in Chapter 1. Using the potential where the peak occurs, the average ion velocity jvij is obtained directly. Other values are found by assuming the same Gaussian ion density pro le used in the modeling study (Section 6.1.2) with & = 2 10;4 and setting the fraction of doubly-charged ions . The half-angle of divergence is assumed to be 10 degrees in most of the simulations, because it led to good agreement in the previous chapter. Other simulations use a 20 degree angle. In this study, the fraction is assumed to be 11% based on King's molecular beam mass spectroscopy (MBMS) measurements. 18] As with the previous studies, the ion temperature is assumed to be 4 eV at the thruster exit, and the background pressure of the test facility, pb = 6 mPa, is included. The densities and velocities based on these assumptions are presented in Table 7.1. These values give a current of 3.5 A and a utilization e ciency of 86.5% (based on the main ow rate) at the thruster exit.
133
Table 7.1. Values at the thruster exit for reference (Ref) case.
Species Xe Xe+ Xe++
Density Velocity Temperature
1017m;3] km/s]
eV]
17.1 0.325
0.086
1.91
19.1
4
0.236 27.1
4
7.2 Simulations The previous simulations included only neutral-neutral collisions and ion-neutral charge-exchange and momentum transfer collisions. As discussed in Chapter 3, each collision type can be described by its cross section. Two expressions for ionneutral charge exchange are discussed { the Rapp and Francis model 53] used in the previous studies and the Sakabe and Izawa model 56] for comparison. The Sakabe and Izawa expression can be extended to doubly-charged ion collisions with neutrals, whereas the Rapp and Francis model is used together with that of Hasted and Hussein 55] for doubly-charged ions. Ion-ion collisions are also considered in this study. Simulations are performed which model momentum transfer collisions using the Rutherford model. 57] Other simulations model ion-ion chargeexchange collisions using the Sakabe and Izawa expression for ion-neutral charge exchange. 56] Most of these simulations use this expression scaled by two orders of magnitude in order to see any signi cant e ect on the distribution. Both types of ion-ion collisions are included in some simulations in order to see if the momentum transfer collisions smooth out the multi-peaked distribution caused by the ion-ion
134 Table 7.2. Input conditions for various simulations.
Case n;i i;i MISC.
Ref* RF 0.0
-
VT* RF 0.0
rTe
CC RF 0.0
CC
CX1* SI 100 SI
-
CX2 SI 100 SI 20 div.
ALL* SI 100 SI
CC
r VT2* SI 100 SI 20 div., Te
charge-exchange reactions. The cross sections for these collisions can be found in Chapter 3. The distribution of the high energy ions in the plume depends on their divergence from the annular acceleration chamber. Therefore, the various collision models are included in simulations with and without a varying electron temperature. Two half-angles of divergence are also considered. A complete list of simulations performed is found in Table 7.2. The asterisks indicate that the computational domain included the full chamber geometry.
7.3 Results Using the Sakabe and Izawa expression instead of that by Rapp and Francis for the ion-neutral charge-exchange cross section does not signi cantly alter the ion distribution function. Therefore, results are presented for the reference case, labeled
135 \Ref" in Table 7.2 and the gures. Presented in Fig. 7.1 are the distribution functions from the reference case at a radial distance of 50 cm from the thruster on the centerline and at 20 degrees o axis along with the corresponding measurements. In this gure as well as the subsequent comparisons, the results are presented normalized, because the peaks are dependent upon the area over which the current is measured and the discretization of the energy. The location of the maximum current near 250 V for the simulation and the experimental data along with the agreement in the width of this peak justi es the assumptions made in the computations for the thruster exit plane conditions. A peak is not found at 500 V as may be expected for the double ions, because the sampling of the distribution using an electrostatic energy analyzer (ESA) gives the energy per charge of the ion. 18] The signi cant tail extending from 400 to 600 V found experimentally is likely caused by collisional behavior inside the plume. The structure of the experimental data is not seen in the simulation. King suggests that the doubly-charged ions undergo charge-exchange interactions with the singly-charged ions. 18] The gradual increase in current with energy before the peak is presumed to be due to slower doubly-charged ions resulting from these collisions. Therefore, simulations are performed in which these ion collisions are permitted. A simulation in which the ion-neutral charge-exchange cross section is used for these collisions is not signi cantly di erent from the reference case. The frequency of the ion-neutral charge exchange is much larger than that for ion-ion charge exchange, when the cross sections are comparable. Both frequencies are proportional to the ion density, but only the ion-neutral frequency scales with neutral density. The background neutrals become the dominant species in the plume as the ion density decays. Simulations are thus performed using the Sakabe
136
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
Ref 0° Ref 20° Data 0° Data 20° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.1. Comparison of reference case with experimental data at r = 50 cm and = 0 and 20 .
and Izawa cross section increased by two orders of magnitude to see the e ects of ion-ion charge exchange. Results from a simulation using this cross section for these reactions are presented in Fig. 7.2 with the reference case on the centerline. These results have been normalized by the reference case peak. The low peak near 125 V is from the singly-charged ions that become doubly-charged ions after charge exchange. The tail centered around 500 V consists of the singly-charged ions which had been doubly-charged ions. A simulation using a 20 degree half-angle of divergence at the thruster exit is presented in Fig. 7.3 with experimental data. The structure is quite similar to case \CX1" which uses a 10 degree divergence angle. The peak at 20 degrees o axis is lower than that found experimentally, but it is higher than that for the reference case. The higher divergence angle leads to more ion spreading. The experimental data do not show a spike near 125 V, but instead show a more gradual
137
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
Ref 0° CX1 0° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.2. Comparison of reference case with one allowing ion-ion CEX collisions at r = 50 cm and = 0 .
increase in the current. This same smooth structure is seen near 400 V where the high energy tail begins. An ion temperature larger than the value assumed in the simulations could lead to this structure. However, the width of the main peaks are in good agreement, suggesting that this is not the case. Another possibility is momentum transfer between various ions. As mentioned previously, the e ects of these types of collisions are examined. Comparisons of the reference case with one which includes Coulomb collisions, in Fig. 7.4, indicates an insensitivity of the ion distribution to these collisions. Similar comparisons are made between the case labeled \CX1" and one which includes the Coulomb collisions with the ion-ion charge-exchange collisions in Fig. 7.5. Again these momentum transfer collisions are not signi cant enough to alter the pro le. Both comparisons are at 20 degrees o axis at the axis, the plots would be nearly identical. Because the Coulomb force is a long-range force, the collision rate is quite high. However, these collisions
138
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
CX2 0° CX2 20° Data 0° Data 20° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.3. Comparison of case with larger divergence angle with experimental data at r = 50 cm and = 0 and 20 .
predominantly result in low scattering angles. Therefore, the negligible e ect on the distribution function is to be expected. The in uence of the divergence angle on the ion distribution in the plume is apparent in Figs. 7.6 and 7.7. Results from the larger divergence angle simulation are presented along with a simulation using the 10 degree half-angle. The larger divergence angle causes more of the beam ions to spread from the axis, thus more high energy ions are sampled 30 degrees from the axis. The other simulation leads to a substantial number of charge-exchange ions being sampled, and the peak is not due to propellent ions. Comparisons at 1 m are similar. By 40 degrees away from the axis, the peak is due to the low energy charge-exchange ions for both simulations. However, results from the larger divergence angle case do show a substantial number of ions around 250 V. Also included in Fig. 7.7 is a simulation using a 20 degree divergence angle which includes the variable electron temperature. This
139
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
Ref 20° CC 20° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.4. Comparison of case with and without Coulomb collisions at r = 50 cm and = 20 .
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
CX1 20° All 20° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.5. Comparison of case with and without Coulomb collisions at r = 50 cm and = 20 .
140
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 00
CX1 30° CX2 30° Data 30° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.6. Comparison of case with larger divergence angle (20 ) to 10 case at r = 50 cm and = 30 .
simulation produces a peak near that of the data. The tail found experimentally is not produced in any of the simulations. By 60 degrees the ion distribution in each of the simulations is dominated by charge-exchange ions, and beam ions are negligible. The experimental data do not show a signi cant number of charge-exchange ions at low energy. As presented in Chapter 6, current density measurements obtained in the near eld by Kim et al 19] show a narrow peak at 10 mm. If too large of a divergence angle is assumed, the resulting peak would be too wide. The agreement at an axial distance of 10 mm in Fig. 7.8 indicates that 20 degrees is in fact not too large. Also, at 100 mm the agreement with the data is very good. Comparisons with current density measurements by King 17] are similar to those found in Chapter 6. Comparisons at a radial distance of 50 cm are presented in Fig. 7.9. Assuming a larger divergence angle does lead to better agreement. However, by 1 m, the
141
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
CX1 40° CX2 40° VT2 40° Data 40° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.7. Comparison of cases with larger divergence angle (20 ) to 10 case at r = 50 cm and = 40 .
computed peak is lower than the data for this simulation. Close examination of these comparisons indicates that the peaks near 250 V in the distribution functions at angles greater than 20 degrees o axis for the simulations should be below those of the experiment. This result can be seen in Figs. 7.1 and 7.3. The agreement for current density is reasonable, but at angles between 20 and 50 degrees the discrepancies are apparent in the distribution functions peaks. In order to illustrate the divergence of the ion beam, a comparison of various simulations at 30 degrees from the axis is presented in Fig. 7.10. The variable electron temperature leads to more spreading than the reference case, but an imposed divergence, as in cases \CX2" and \VT2", has more of an e ect. It has been found that the experimental data do not show as much in uence of the low energy ions created by charge-exchange reactions as the simulations do. Therefore, the decay of the magnitude of the peak current with angle from the axis is examined
142
Current Density (mA/cm2)
80
CX1 10 mm
70
CX1 100 mm
CX2 10 mm
60
CX2 100 mm
Data 10 mm
50
Data 100 mm
40
30
20
10
0
0
0.05
0.1
0.15
0.2
Radial Distance (m)
Figure 7.8. Comparisons of current density from both divergence angle simulations with measurements by Kim.
101 CX1 CX2 Data 100
Current Density (mA/cm2)
10-1
10-2
-50
0
50
Angle (deg)
Figure 7.9. Comparisons of current density from cases allowing ion-ion CEX collisions with measurements by King.
143
f (V)
4 3.5 3 2.5 2 1.5 1 0.5 00
Ref 30° VT 30° CX2 30° VT2 30° 100 200 300 400 500 600 Ion Energy/q (Volts)
Figure 7.10. Comparison of various cases at r = 50 cm and = 30 . in Fig. 7.11. This shows the peak of the ion distribution function in the range 200 to 300 V normalized by the value on the axis at a radial distance of 50 cm for the case labeled \CX2". The experimental data has been normalized by the value measured at this location. At angles above 40 degrees from the centerline, the simulations predict less high energy ions. The case with a 10 degree divergence angle gives results even further from the data. The data show high energy ions up to 100 degrees from the thruster centerline. It is unclear what mechanism causes these propellent ions to divert from an axiallydirected path. Comparison of current density between simulation and experiment in the back ow region shows good agreement. The distribution functions at large angles behind the thruster are presented in Fig. 7.12. At both radial distances 150 degrees away from the axis, the reference case is shown with experimental measurements. There is very little di erence in the two radial positions in the
144
100
CX1
CX2
VT2
Data 10-1
fpeak (V)
10-2
10-3
10-4 0
10
20
30
40
50
60
Angle (deg)
Figure 7.11. Peaks in distribution function vs. angle for various cases with experimental data at r = 50 cm. simulation. The simulation agrees well with the data at 50 cm but not at 1 m. The data shows a peak with higher energy at 1 m than the peak at 50 cm. This suggests an acceleration of the ions. A simulation which includes ion-ion chargeexchange collisions is not noticeably di erent from the reference case. This is not surprising, since it is the charge-exchange ions formed from ion-neutral collisions that reach the region the back ow region.
7.4 Conclusions The ion energy distribution function has been compared with experimental data at various angles from the thruster centerline. The agreement between the simulations and experiments is reasonable at low angles from the centerline for most cases. However, at the larger angles, the agreement is poor. A high energy tail is found
145
f (V)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
Ref 50cm Ref 1m Data 50cm Data 1m
25
50
75
100
Ion Energy/q (Volts)
Figure 7.12. Comparison of reference case with experimental data at r = 50 cm and 1 m and = 150 .
experimentally, but not computationally, at most angles in front of the thruster. This suggests that some physics is not being captured in the simulations. The computed ion distribution function is insensitive to the frequent Coulomb collisions, because they result predominantly in small scattering angles. Also, a cross section for ion-ion charge exchange of comparable magnitude to that of ionneutral charge exchange was insu cient to show the structure in the distribution function found in the measured data. A much larger cross section for ion-ion charge exchange resulted in structure similar to these data. However, this cross section seems unreasonably large. The Coulomb force is likely to prevent these collisions, because the ions would need to be very close to exchange an electron. In a recent study, Hofer et al. list possible e ects on the distribution function by the proximity of the experimental set-up to the chamber walls. 68] Among these are enhanced ow rates and arti cially created charge-exchange collisions due to a local increase
146 in pressure near the set-up. These collisions may lead to signatures in the ion distribution function not indicative of the plume of the thruster. It is also possible that ions are scattered from the chamber walls. The properties of these ions would also alter the structure of the distribution. Spreading of the propellent ions is determined by the ion temperature, the ion divergence angle at the thruster exit, and the electrostatic eld in the plume. The width of the distribution function at low angles agreed well with that found experimentally. Thus, the assumed ion temperature at the thruster exit is reasonable. Simulations using a larger divergence angle led to better agreement with the data away from the thruster centerline. However, it did not lead to good agreement at the higher angles in front of the thruster. Comparisons with near eld current density measurements indicated that spreading of the ion beam does occur in the plume, but the divergence half-angle is unlikely to exceed 20 degrees. A variable electron temperature model had the e ect of spreading the ions in the far plume while permitting a narrow beam in the near eld, as the measured data suggest. More detailed modeling combined with a variable electron temperature is necessary to reproduce the data in this study. As mentioned in Chapter 6, the magnetic eld in the near plume is likely to a ect the electrons. The ion behavior is in uenced by the electrons, thus the magnetic eld may alter the ion pro le from its shape at the thruster exit. Modeling of the conditions found in the vacuum chamber where the experimental data were obtained is also necessary. High energy ions may scatter from the chamber walls. Some of the ions which are scattered could be collected by the MBMS system. Another possibility for the structure found experimentally is that the pro le of the ion beam is more complicated at the thruster exit than has been assumed in the
147 simulations. The ions in the acceleration chamber have an azimuthal component to their velocity that has been neglected. The ion velocity distribution may not be a Maxwellian at the exit. Since all ions are not created in the same place in the chamber, the properties may vary signi cantly. Although most of the ion beam is enclosed within a xed divergence angle, some ions may leave the thruster with a large radial velocity.
Chapter 8 Conclusions A computer code which combines the direct simulation Monte Carlo and the Particle-in-Cell methods was developed for computation of plumes from ion thrusters and Hall thrusters. This parallelized code uses a weight by particle technique to simulate directly the properties of ions and neutral atoms. Both charge-exchange and momentum transfer collisions between these species are included. An implicit grid scheme permits a large number of PIC grid cells without inordinate memory costs. The use of separate grids for the two numerical methods permits only limited interaction between the algorithms. Because the grids must share the particles, the velocity and position of the particles must be updated in the DSMC grid based on PIC information, and the particle mapping for charge density at the PIC nodes becomes slightly more cumbersome if the simulation uses multiple processors. A routine was developed to include the e ects of the back pressure of the vacuum chamber in which experimental data were obtained. These features allowed physical models to be tested and o er the capability for further investigations. In this study, calculations were compared with experimental data. These data were obtained for the plumes of the UK-10 ion thruster and the SPT-100 Hall thruster. The ion thruster plume study successfully tested the computer code. Qualitatively the computer calculations are encouraging, but di erences between the simulations and the data are apparent. The Hall thruster simulations agreed 148
149 well with the plume data. This study also illustrated some important features in the plume. 8.1 Summary of Conclusions 8.1.1 Ion thruster The ion thruster simulations captured the ion dynamics qualitatively. The beam ions spread more due to the electric elds than purely thermal e ects would cause. The charge-exchange ions were pulled back towards the thruster near the exit and away from the beam further away. This study found a dependence on various physical assumptions. The plume is sensitive to the electron temperature value, the charge-exchange cross section, and the collision dynamics of charge-exchange collisions. Comparison with experimental data indicated the importance of the ion pro le at the thruster exit plane. The focusing of the ion beam by the concave acceleration grids leads to a more structured ion beam than the simulations show. Charge-exchange collisions between the grids may be important. Ions may also be re ected specularly by the grids leading to large scattering angles. A better understanding of the e ect of the accelerator grids is required to de ne the ion pro le and thus to improve the agreement with the data. 8.1.2 Hall thruster Calculations of far eld ion properties are a ected by variations in the ion inlet pro le. However, the comparisons with experimental data show only moderate di erences. Each of the various simulations agreed well with experimental measurements of far eld ion current density. These variations included half-angles of
150 divergence of 0, 10, and 20 degree and two di erent ion temperatures. Including a varying electron temperature was found to be signi cant in the modeling of Hall thruster plumes. The electron temperature gradient signi cantly a ects the near eld of the plume. Simulations incorporating the varying electron temperature model gave better agreement with the near eld current density measurements than those using the isothermal Boltzmann relation. These simulations also led to better agreement with the electron number density and the slope of the potential measurements. Discrepancies between simulations and experiment in the near eld also suggest a missing feature in the modeling. The experimental data indicate an ion beam focusing earlier at the axis than the simulations predict. In the simulations, the half-angle of divergence can be set to any value to focus the beam at the axis. However, comparison of current density in the near and far eld indicates this value is unlikely to exceed 20 degrees. The magnitude of the electron temperature a ects the spreading of the ions in the plume. Moreover, results suggest that the plasma behavior in the far eld region is not dictated by the near eld behavior, but rather by the value of the electron temperature. Therefore, the far eld electron temperature is important. This temperature is likely to be lower than the near eld value. The simulations agreed well with the measurements of electron number density when lower values of electron temperature were used for the far eld. Agreement with plasma potential was dependent upon the value for the reference density used in the physical models as well as the electron temperature. Simulating the full geometry of the experimental facility is important in order to use measurements of the back ow region for code validation. The simulations
151 were found to represent the current density very well. Even at large angles behind the thruster good agreement was obtained with experiment. The back pressure of the experimental facilities a ects the plume. It is necessary to include this interaction in order to reproduce the current density measured in the back ow region. The level of agreement indicates that the routine used to simulate the background pressure is adequate. Analysis of the ion energy distribution function reveals information about the collisions in the plume. The ion thruster study indicated the signi cance of the charge-exchange cross section on the ion density. This study also found that simulations assuming no momentum transfer for charge-exchange collisions match ion density measurements better than a simulation in which momentum transfer was assumed. The Hall thruster analysis led to more information about ion collision phenomena. Inclusion of ion-ion momentum transfer collisions in the simulations had little e ect on the computed ion distribution in the plume, because they were predominantly small angle collisions. Comparisons of the distribution function with experimental measurements showed discrepancies. The structure found experimentally was not reproduced in the simulations. This structure includes a large tail near twice the peak voltage. Even at large angles o axis, these high energy ions are found. Therefore, some e ect of the distribution obtained experimentally is not being captured. The structure of these data is unlikely due to ion-ion chargeexchange collisions. Simulations using an unrealistically large cross section gave only minimal similarity, and the structure was not found in the simulations at large angles o axis. Based on the results from the simulations, the following parameters and models are recommended to compute the Hall thruster plume. The variable electron
152 temperature model with a far eld electron temperature of 2 eV should be used. The ion inlet pro le presented in Chapter 7 is recommended along with this model. This pro le included using an ion kinetic energy of 250 eV and the modi ed Gaussian for the density with a maximum half-angle of divergence of 10 degrees to de ne the densities and velocities at the thruster exit. The charge-exchange cross section for doubly-charged ions given by Sakabe and Izawa 56] is more realistic than the Hasted and Hussain 55] expression in this energy range. There is very little di erence in using the Rapp and Francis 53] expression or the Sakabe and Izawa 56] expression for the singly-charged ions. If the simulations are to match conditions in an experimental facility, the background pressure of the facility must be included. 8.2 Future Work Although the simulations did a reasonable job of representing the plumes of the Hall thruster, comparisons with the experimental data also indicate a need for additional modeling e orts. One area of concern is the electron modeling. Improvements can be made to this modeling without simulating the electrons as particles. The imposed electron temperature indicated that the variations signi cantly a ect the plume, especially in the near eld. However, the values employed in these simulations are only representative of the temperature. Moreover, they only apply to this particular Hall thruster at these operating conditions. To properly obtain the electron temperature variation would require calculating the temperature in the plume. With additional computational e ort the electron energy equation can be used to obtain the temperature.
153 The near eld current density calculations di er from the experimental data. It is likely that the assumptions for the electron modeling are less valid in the near eld region. For example, the magnetic eld lines do not end at the thruster exit. The e ect of this magnetic eld in the plume is an issue which should be addressed. As discussed in Chapter 3, the signi cance of this term can be examined by varying the ratio of the cyclotron frequency to the collision frequency. However, this approach is not as physical as a more in depth analysis of the behavior of the magnetized electrons. The computer code used in this study o ers the exibility to include these modi cations to the electron modeling. Routines can be added to solve for the electron temperature, or a more complex equation can replace the simpli ed electron momentum equation employed in the simulations. The e ect of the acceleration chamber of the Hall thruster on the ions is another area for further investigation. As discussed in Chapter 7, the ion pro le at the exit may exhibit some of the structure found experimentally in the distribution function measured in the plume. Plasma oscillations inside the chamber are observed experimentally using optical bers for the SPT-100. 69,70] At frequencies on the order of tens of kHz, uctuations are observed in ion density and ion mean energy. 69] Even the measured ion energy distribution function in the plume exhibits a time dependence. 70] To better understand the e ects on the ion pro le, a detailed analysis of the interior of the Hall thruster is needed. Performing full chamber geometry simulations allows the chamber e ects on the experimental data to be examined. It was found that the background pressure of the vacuum chambers must be included in the simulations. The background gas leads to more charge-exchange ions in the back ow region than would otherwise be there. Although the chamber walls were set to ground potential in the simulations,
154 their e ect on the domain needs to be examined in more detail. Other wall e ects may include sputtering of material by high energy ions, scattering of ions back into the plume, and heat transfer with the gas. These e ects have not been modeled. Determining how the test facilities a ect the data is important. If the extent of the alteration can be quanti ed, the data can be adjusted accordingly in order to estimate the plume properties that would be observed without the facility in uence. Simulations allow possible chamber e ects to be examined. They can also be used to try to represent conditions found in space. If the chamber e ects on the plume are determined, predictions of the plume behavior in space can be made both by simulation and experiment. However, it is easier to specify background conditions in the simulations than it would be to control the chamber conditions.
Bibliography 1] Crofton M.W., \Evaluation of the United Kingdom Ion Thruster," Journal of Spacecraft and Rockets, Vol. 33, No. 5, 1996, pp. 739-747. 2] Sutton, G.P., Rocket Propulsion Elements, 6th Edition, Wiley-Interscience, 1992. 3] Polk, J.E., Kakuda, R.Y., Anderson, J.R., and Brophy, J.R., \Validation of the NSTAR Ion Propulsion System on the Deep Space One Missin: Overview and Initial Results," AIAA Paper 99-2274, June 1999. 4] Stuhlinger, E., Ion Propulsion for Space Flight, McGraw-Hill, 1964. 5] Fearn, D.G., Martin, A.R., and Smith, P., \Ion Propulsion Research and Development in the UK," Journal of British Interplanetary Society, Vol. 43, No. 10, 1990. 6] Jahn, R.G., Physics of Electric Propulsion, McGraw-Hill, 1968. 7] Roy, R.S., \Numerical Simulation of Ion Thruster Plume Back ow for Spacecraft Contamination Assessment," Ph. D. thesis, MIT, 1995. 8] Ashkenazy, J., Fruchtman, A., Raitses, Y., and Fisch, N.J., \Modelling the Behaviour of a Hall Current Plasma Accelerator," Plasma Physics and Controlled Fusion, Vol. 41, 1999, pp. A357-A364. 9] Zhurin, V.V., Kaufman, H.R., and Robinson, R.S., \Physics of Closed Drift Thrusters," IEPC Paper 97-191, August 1997. 10] Shumlak, U., Jarboe, T.R., and Sprenger, R.A., \Physics of the Hall Thruster," AIAA Paper 97-3048, July 1997. 11] de Boer, P.C.T., \Electric Probe Measurements in the Plume of an Ion Thruster," Journal of Propulsion and Power, Vol. 12, No. 1, Jan.-Feb. 1996, pp. 95-104. 155
156 12] de Boer, P.C.T., \Measurements of Electron Density in the Charge Exchange Plasma of an Ion Thruster," Journal of Propulsion and Power, Vol. 13, No. 6, November-December 1997, pp. 783-788. 13] Pollard, J.E., \Plume Angular, Energy, and Mass Spectral Measurements with the T5 Ion Engine," AIAA Paper 95-2920, July 1995. 14] Pollard, J.E., \Pro ling The Beam of The T5 Ion Engine," IEPC Paper 97019, August 1997. 15] Manzella, D.H., \Stationary Plasma Thruster Ion Velocity Distribution," AIAA Paper 94-3141, August 1994. 16] Manzella, D.H. and Sankovic, J.M., \Hall Thruster Ion Beam Characterization," AIAA Paper 95-2927, July 1995. 17] King, L.B., Gallimore, A.D., and Marrese, C.M., \Transport-Property Measurements in the Plume of an SPT-100 Hall Thruster," Journal of Propulsion and Power, Vol. 14, No. 3, 1998, pp. 327-335. 18] King, L.B., \Transport-Property and Mass Spectral Measurements in the Plasma Exhaust Plume of a Hall-E ect Space Propulsion System," Ph. D. thesis, University of Michigan, 1998. 19] Kim, S., Foster, J.E., and Gallimore, A.D., \Very-Near-Field Plume Study of a 1.35 kW SPT-100," AIAA Paper 96-2972, July 1996. 20] Sankovic, J.M., Hamley, J.A., and Haag, T.W., \Performance Evaluation of the Russian SPT-100 Thruster at NASA LeRC," IEPC Paper 93-094, September 1993. 21] Myers, R.M. and Manzella, D.H., \Stationary Plasma Thruster Plume Characteristics," IEPC Paper 93-096, July 1993. 22] Marrese, C. and Gallimore, A., private communications. 23] Roy, R.S., Hastings, D.E., and Taylor, S., \Three-Dimensional Plasma Particle-in-Cell Calculations of Ion Thruster Back ow Contamination," Journal of Computational Physics, Vol. 128, 1996, pp. 6-18. 24] Oh, D. and Hastings, D.E., \Three Dimensional PIC-DSMC Simulations of Hall Thruster Plumes and Analysis for Realistic Spacecraft Con gurations," AIAA Paper 96-3299, July 1996. 25] Oh, D. and Hastings, D.E., \Axisymmetric PIC-DSMC Simulations of SPT Plumes," IEPC Paper 95-160, September 1995. 26] Gatsonis, N.A. and Yin, X., \Axisymmetric DSMC/PIC Simulation of Quasineutral Partially ionized Jets," AIAA Paper 97-2535, July 1997.
157 27] Wang, J., Brophy, J., and Brinza, D., \3-D Simulations of NSTAR Ion Thruster Plasma Environment," AIAA Paper 96-3202, July 1996. 28] Keidar, M. and Boyd, I.D., \E ects of a Magnetic Field on the Plasma Plume from Hall Thrusters," Journal of Applied Physics, Vol. 86, No. 9, 1999, pp. 4786-4791. 29] Oh, D. and Hastings, D.E., \Experimental Veri cation of a PIC-DSMC Model for Hall Thruster Plumes," AIAA Paper 96-3196, July 1996. 30] Bird, G.A., Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, 1994. 31] Birdsall, C.K. and Langdon, A.B., Plasma Physics Via Computer Simulation, Adam Hilger, 1991. 32] Boyd, I.D. and Gokcen, T., \Computation of Axisymmetric and Ionized Hypersonic Flows Using Particle and Continuum Methods," American Institute of Aeronautics and Astronautics Journal, Vol. 32, No. 9, 1994, pp. 1825-1835. 33] Boyd, I.D., Candler, G.V., and Levin, D.A., \Dissociation Modeling in Low Density Hypersonic Flows of Air," Physics of Fluids, Vol. 7, No. 7, 1995, pp. 1757-1763. 34] Boyd, I.D., VanGilder, D.B., and Beiting, E.J., \Computational and Experimental Investigations of Rare ed Flows in Small Nozzles," AIAA Journal, Vol. 34 (11), 1996, pp. 2320-2326. 35] Oh, C.K., Oran, E.S., and Cybyk, Z.C., \Microchannel Flow Computed with the DSMC-MLG," AIAA Paper 95-2090, June 1995. 36] VanGilder, D.B., Font, G.I., and Boyd, I.D., \Hybrid Monte Carlo-Particle-inCell Simulation of an Ion Thruster Plume," Journal of Propulsion and Power, Vol. 15 (4), 1999, pp. 530-538. 37] Boyd, I.D., VanGilder, D.B., and Liu, X., \Monte Carlo Simulation of Neutral Xenon Flows in Electric Propulsion Devices," Journal of Propulsion and Power, Vol. 14 (6), 1998, pp. 1009-1015. 38] Chen, G., Boyd, I.D., Roadman, S.E., and Engstrom, J.R., \Monte Carlo Analysis of a Hyperthermal Silicon Deposition Process," Journal of Vacuum Science and Technology A-Vacuum Surfaces and Films, Vol. 16, No. 2, 1998, pp. 689-699. 39] Bagano , D. and McDonald, J.D., \A Collision Selection Rule for a Particle Simulation Method Suited to Vector Computers," Physics of Fluids, Vol. 2, No. 7, 1990, pp. 1248-1259.
158 40] Boyd, I.D., \Rotational-Translational Energy Transfer in Rare ed Nonequilibrium Flows," Physics of Fluids A, Vol. 2, No. 3, 1990, pp. 447-452. 41] Boyd, I.D., \Analysis of Vibrational-Translational Energy Transfer Using the Direct Simulation Monte Carlo Method," Physics of Fluids A, Vol. 3, No. 7, 1991, pp. 1785-1791. 42] Abe, T., \Inelastic Collision Model for Vibrational-Translational and Vibrational-Vibrational Energy Transfer in the Direct Simulation Monte Carlo Method," Physics of Fluids, Vol. 6, No. 9, 1994, pp. 3175-3179. 43] Chen, G. and Boyd, I.D., \Statistical Error Analysis for the Direct Simulation Monte Carlo Technique," Journal of Computational Physics, Vol. 126, No. 2, 1996, pp. 434-448. 44] Kannenberg, K., \Computational Methods for the Direct Simulation Monte Carlo Technique with Application to Plume Impingement," Ph. D. thesis, Cornell University, May 1998. 45] Ruyten, W.M., \Density-Conserving Shape Factors for Particle Simulations in Cylindrical and Spherical Coordinates," Journal of Computational Physics, Vol. 105, 1993, pp. 224-232. 46] Myers, R.M., Pencil, E.J., Rawlin, V.K., Kussmaul, M., and Oden, K., \NSTAR Ion Thruster Plume Impacts Assessments," AIAA Paper 95-2825, July 1995. 47] Mitchner, M. and Kruger, C.H. Jr., Partially Ionized Gases, John Wiley and Sons, 1973. 48] Hutchinson, I.H., Principles of Plasma Diagnostics, Cambridge University Press, 1987. 49] Chen, F.F., Introduction to Plasma Physics and Controlled Fusion, Plenum Press, 1984. 50] Koura, K. and Matsumoto, H., \Variable Soft Sphere Molecular Model for Air Species," Physics of Fluids A, Vol. 4, No. 5, May 1992, pp. 1083-1085. 51] Dalgarno, A., McDowell, R.C., and Williams, A., \The Mobilities of Ions in Unlike Gases," Phil. Trans. of the Royal Society A, Vol. 250, 1958, pp. 411-425. 52] Banks, P., \Collision Frequencies and Energy Transfer Ions," Planetary and Space Science, Vol. 14, 1966, pp. 1105-1122. 53] Rapp, D. and Francis, W.E., \Charge Exchange Between Gaseous Ions and Atoms," Journal of Chemical Physics, Vol. 37, No. 11, 1962, pp. 2631-2645.
159 54] Dalgarno, A., \The Mobilities of Ions in Their Parent Gases," Phil. Trans. of the Royal Society A, Vol. 250, 1958, pp. 426-439. 55] Hasted, J.B. and Hussain, M., \Electron Capture by Multiply Charged Ions," Proceedings of the Physical Society, Vol. 83, 1964, pp. 911-924. 56] Sakabe, S. and Izawa, Y., \Simple Formula for the Cross Sections of Resonant Charge Transfer between Atoms and Their Positive Ions at Low Impact Velocity," Physical Review A, Vol. 45, No. 3, February 1992, pp. 2086-2089. 57] Ichimaru, S., Statistical Plasma Physics, Addison-Wesley, 1992. 58] Dietrich, S. and Boyd, I.D., \Scalar and Parallel Optimized Implementation of the Direct Simulation Monte Carlo Method," Journal of Computational Physics, Vol. 126, 1996, pp. 228{243. 59] Peraire, J., Vahdati, M., Morgan, K., and Zienkiewicz, O.C., \Adaptive Remeshing for Compressible flow computations," Journal of Computational Physics, Vol. 72, 1987. 60] Jaluria, Y. and Torrance, K.E., Computational Heat Transfer, Hemisphere Publishing Corporation, 1986. 61] Karipides, D.P., \Detailed Investigation of Spacecraft Glow," Ph. D. thesis, Cornell, May 1999. 62] Boyd, I.D. \Conservative Species Weighting Scheme for the Direct Simulation Monte Carlo Method," Journal of Thermophysics and Heat Transfer, Vol. 10, No. 4, October-December 1996, pp. 579{585. 63] Dietrich, S. and Boyd, I.D., \Loadbalancing for the DSMC Method on Parallel Computers," Proceedings of the 19th International Symposium on Rare ed Gas Dynamics, Oxford University Press, 1995. 64] Crofton, M.W., \Measurement of Neutral Xenon Density in an Ion Thruster Plume," AIAA Paper 96-2290, June 1996. 65] Wang, J., Anderson, J., and Polk, J., \Three-Dimensional Particle Simulations of Ion Optics Plasma Flow," AIAA Paper 98-3799, July 1998. 66] Peng, X., Keefer, D., and Ruyten, W., \Plasma Particle Simulation of Electrostatic Ion Thrusters," Journal of Propulsion and Power, Vol. 8, No. 2, 1992, pp. 361-366. 67] Cedolin, R.J., Hargus Jr., W.A., Hanson, R.K., and Capelli, M.A., \LaserInduced Fluorescence Diagnostics for Xenon Hall Thrusters," AIAA Paper 96-2986, July 1996.
160 68] Hofer, R.R., Haas, J.M., and Gallimore, A.D., \Development of a 45-Degree Parallel-Plate Electrostatic Energy Analyzer for Hall Thruster Plume Studies: Preliminary Data," IEPC Paper 99-113, October 1999. 69] Darnon, F., Kadlec-Philippe, C., Bouchoule, A., and Lyszyk, M, \Dynamic Plasma and Plume Behavior of SPT Thrusters," AIAA Paper 98-3644, July 1998. 70] Pagnon, D., Darnon, F., Roche, S., Bechu, S., Magne, L., Bouchoule, A., Touzeau, M., and Lasgorceix, \Time Resolved Characterization of the Plasma and the Plume of a SPT Thruster," AIAA Paper 99-2428, June 1999.

DB VanGilder

File: numerical-simulations-of-the-plumes-of-electric-propulsion-thrusters.pdf
Title: thesis.dvi
Author: DB VanGilder
Published: Thu Oct 31 07:45:32 2002
Pages: 180
File size: 1.18 Mb


Essential Elements, 3 pages, 0.17 Mb

TRIZ-what is TRIZ, 3 pages, 0.12 Mb

, pages, 0 Mb

Antioxidant Superfoods, 28 pages, 0.8 Mb
Copyright © 2018 doc.uments.com