Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UFE0050292/00001
## Material Information- Title:
- The Rarefied Gas Electro Jet (RGEJ) Micro-Thruster for Propulsion of Small Satellites
- Creator:
- Blanco, Ariel
- Place of Publication:
- [Gainesville, Fla.]
Florida - Publisher:
- University of Florida
- Publication Date:
- 2017
- Language:
- english
- Physical Description:
- 1 online resource (176 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Aerospace Engineering
Mechanical and Aerospace Engineering - Committee Chair:
- ROY,SUBRATA
- Committee Co-Chair:
- FITZ-COY,NORMAN G
- Committee Members:
- BALACHANDAR,SIVARAMAKRISHNAN
MOORE,ROBERT C
## Subjects- Subjects / Keywords:
- drift-diffusion-approximation -- electrothermal-thrusters -- micro-propulsion -- plasma -- rarefied-gas
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF - Genre:
- bibliography ( marcgt )
theses ( marcgt ) government publication (state, provincial, terriorial, dependent) ( marcgt ) born-digital ( sobekcm ) Electronic Thesis or Dissertation Aerospace Engineering thesis, Ph.D.
## Notes- Abstract:
- Small satellites have multiple applications such as commercial, military, and space science missions. Further development of micro-propulsion concepts is needed. In this dissertation, a micro-thrusters for small satellites called the rarefied gas electro jet (RGEJ) is investigated numerically. The RGEJ design aims to increase the efficiency of the energy exchange from the electrical source to the propellant, avoiding the high thermal dissipation problem associated with other technologies. In order to simulate the RGEJ, an existing and previously described in literature, in-house modular Multi-scale Ionized Gas (MIG) flow solver platform is used. The rarefied gas is modeled using density-based compressible flow equations with rarefied boundary conditions and the ionized gas is modeled using local mean energy approximation (LMEA). Benchmarking and validation of the rarefied gas and ionized gas modules are performed, respectively. Several thruster cases are studied with different discharge voltages to investigate the discharge characteristics and the performance parameters of the thruster such as mass flow rate, thrust, power consumption, and specific impulse. Finally, a thermal analysis is performed to investigate the thermal losses. For the highest voltage tested (750 V), operating with argon propellant at plenum pressures of 1 Torr, the RGEJ thruster requires a total electrical power of 406 mW to heat the flow. Under these operation conditions, the RGEJ produces a thrust of 1.843 mN and a specific impulse (Isp) of 60.7 s with a thrust effectiveness of 1240 (micro-N/W). This Isp is a 37.6 percent improvement over the Isp of a cold gas thruster operating with the same geometry, propellant, and plenum pressure. The heat loss due to conduction and radiation is expected to be less than 7.5 percent of the total electrical power. The Isp of the RGEJ thruster was found to be 35 percent higher than the competing technology called the free molecule micro-resistojet (FMMR). Base on the thermal analysis, the RGEJ is predicted to have significantly fewer heat losses than FMMR. To the knowledge of the author, a numerical simulation of a thruster design using plasma-aided technology to replace the energy transfer mechanism of the FMMR thruster and decrease its inherent heat losses has never been attempted before. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2017.
- Local:
- Adviser: ROY,SUBRATA.
- Local:
- Co-adviser: FITZ-COY,NORMAN G.
- Electronic Access:
- RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2018-06-30
- Statement of Responsibility:
- by Ariel Blanco.
## Record Information- Source Institution:
- UFRGP
- Rights Management:
- Applicable rights reserved.
- Embargo Date:
- 6/30/2018
- Classification:
- LD1780 2017 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

PAGE 1 THE RAREFIED GAS ELECTRO JET (RGEJ ) MICRO THRUSTER FOR PROPULSION OF SMALL SATELLITES By ARIEL BLANCO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 201 7 PAGE 2 201 7 Ariel Blanco PAGE 3 To my family, friends, and mentors PAGE 4 4 ACKNOWLEDGMENTS I thank my parents and maternal grandmother for their continuous support and encouragement, for teaching me to persevere and value the pursuit of knowledge. I am extremely thankful to my friends and lab colleagues wh o have helped me so much, not only academically but also by giving me moral support and solidarity. Finally, yet importantly, it is a genuine pleasure to express my deep sense of thanks and gratitude to all my mentors, especially my advisor, Dr. Subrata Ro y. Without his guidance, mentoring, and support, this study would not have been completed. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 ABSTRACT ................................ ................................ ................................ ................... 12 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 14 The Rarefied Gas Electro Jet (RGEJ) Design and Motivation ................................ 18 Performance Parameters ................................ ................................ ........................ 25 2 MICROPROPULSION TECHNOLOGIES ................................ ............................... 29 Introduction to Micro propulsion Devices ................................ ................................ 29 Chemical Propulsion Devices for Small Satellites ................................ ................... 29 Cold Gas Thrusters ................................ ................................ .......................... 30 Mono Prope llant Thrusters ................................ ................................ ............... 31 Bipropellant Thrusters ................................ ................................ ...................... 33 Electric Propulsion Devices for Small Satellites ................................ ...................... 34 Electrothermal Thrusters ................................ ................................ .................. 34 Free molecule micro resistojet (FMMR) ................................ ...................... 35 Low power DC arcjets ................................ ................................ ................. 37 Microcavity discharge thrusters ................................ ................................ ... 38 Micro plasma thruster (MPT) ................................ ................................ ...... 39 Electrostatic Thrusters ................................ ................................ ...................... 40 Miniature ion engines ................................ ................................ .................. 41 Miniature hall thrusters ................................ ................................ ................ 43 Coll oidal and field emission electric propulsion (FEEP) thrusters ............... 46 Electromagnetic Thrusters ................................ ................................ ................ 49 Conclusion of Micro propulsion Technologies Survey ................................ ............ 50 3 GAS DISCHARGE PHYSICS ................................ ................................ ................. 52 Elec trical Breakdown of Gases ................................ ................................ ............... 53 Dark Discharge ................................ ................................ ................................ 54 Glow Discharge ................................ ................................ ................................ 56 Arc Discharge ................................ ................................ ................................ ... 57 Qualitative Ch aracteristics of DC Glow Discharge ................................ .................. 57 ................................ ................................ ................................ ......... 62 PAGE 6 6 Sputtering ................................ ................................ ................................ ............... 64 4 GOVERNING EQUATIONS ................................ ................................ .................... 66 Hydrod ynamic Plasma Model ................................ ................................ ................. 66 Boltzmann Equation ................................ ................................ ......................... 66 Distribution Function Moments Equation ................................ .......................... 68 Continuity Equation ................................ ................................ .......................... 70 Momentum Equation ................................ ................................ ........................ 72 Energy Equation ................................ ................................ ............................... 74 Drift Diffusion Approximation ................................ ................................ ............ 76 ................................ ................................ ......................... 81 Transport Properties ................................ ................................ ......................... 84 Boundary Conditions ................................ ................................ ........................ 87 Open surface boundary ................................ ................................ ............... 88 Anode boundary ................................ ................................ .......................... 89 Cathode bounda ry ................................ ................................ ...................... 90 Dielectric wall ................................ ................................ .............................. 90 Plasma Electrostatic Force and Electrothermal Heating Source ...................... 91 Rarefied Gas Equations ................................ ................................ .......................... 92 Heat Transfer Governing Equations ................................ ................................ ....... 97 Thermal Radiation Heat Transfer ................................ ................................ ..... 97 View Factor Formulas ................................ ................................ ...................... 98 Conduction Heat Transfer ................................ ................................ .............. 100 Co nduction Heat Transfer Boundary Conditions ................................ ............ 101 5 NUMERICAL MODEL ................................ ................................ ........................... 102 The Ionized Gas Module ................................ ................................ ....................... 103 The Rarefied Gas Module ................................ ................................ ..................... 107 Gale rkin Weak Statement ................................ ................................ ............... 108 Finite Element Basis Function ................................ ................................ ........ 110 Bilinear Basis ................................ ................................ ................................ .. 112 Newton Raphson Scheme ................................ ................................ .............. 115 Decoup led Heat Transfer Codes ................................ ................................ .......... 117 Conduction Numerical Code ................................ ................................ ........... 117 Thermal Radiation Heat Transfer Numerical Code ................................ ......... 117 6 BENCHMARKING AND VALIDATION ................................ ................................ .. 119 Ionized Gas Module (IGM) ................................ ................................ .................... 119 Verification and Grid Convergence Study of the IGM ................................ ..... 119 Validation of the IGM ................................ ................................ ...................... 122 Rarefied Gas Module (RGM) ................................ ................................ ................ 125 Ve rification of the RGM using plane Poiseuille flow ................................ ....... 125 Benchmarking/Validation of the RGM ................................ ............................ 129 PAGE 7 7 7 CASES STUDIED ................................ ................................ ................................ 132 Geometry and Grid ................................ ................................ ............................... 132 Cold gas t hruster and constant thermal heating source thrusters results and comparison ................................ ................................ ................................ ........ 134 RGEJ thruster performance ................................ ................................ .................. 145 RGEJ thruster discharge characteristics ................................ ............................... 156 Thermal analysis ................................ ................................ ................................ ... 159 8 SUMMARY ................................ ................................ ................................ ........... 162 Summary of the Results ................................ ................................ ........................ 162 The Numerical Model ................................ ................................ ............................ 165 Future D irections to Continue Improving RGEJ ................................ .................... 166 LIST OF REFERENCES ................................ ................................ ............................. 167 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 176 PAGE 8 8 LIST OF TABLES Table page 1 1 Mas s classification of Satellites ................................ ................................ ......... 17 2 1 Cold gas propellant performance. ................................ ................................ ....... 31 2 2 Micro ion engine performance. ................................ ................................ ........... 43 2 3 Small and micro Hall thruster performances. ................................ ...................... 46 2 4 ST 7 colloid thruster cluster performance characteristics. ................................ .. 47 2 5 FEEP Performance Characteristic s. ................................ ................................ ... 48 2 6 Micro PPT Performances. ................................ ................................ .................. 50 3 1 Breakdown voltage constants A and B. ................................ .............................. 63 3 2 Limit of (pd) and minimum breakdown voltage values for argon. ........................ 63 4 1 Relevant Reactions. ................................ ................................ ........................... 86 6 1 Coefficients for the grid refinement study. ................................ ........................ 121 6 2 Compariso n between experiment and numerical model with ( sec =0.01 ) for (300 V < V Pk Pk < 800 V). ................................ ................................ ................... 124 6 3 Microchannel Properties of Fluid for Subsonic Gas Flows. .............................. 130 6 4 Grid dependence test of the centerline u velocity (m/s) at different x locations. ................................ ................................ ................................ .......... 130 6 5 Grid dependence test done using MIG of the centerline u velocity (m/s) at different x locations. ................................ ................................ ......................... 130 6 6 Comparison of centerline pressure ratios at different x locations. .................... 131 7 1 Values of total thermal heating source ( ), mass flow rate ( ), thrust ( F Thrust ), shear force ( ), specific impulse ( I sp ), specific impulse increase ( I sp Inc.) compared to the base case exit plane Mach number ( M Exit ) at the centerline, and exit Knudsen number ( Kn Exit ) at the centerline are displayed .. 136 7 2 Values of: voltage (V), mass flow rate ( ), thrust ( F Thrust ), shear force ( ) specific impulse ( I sp ), specific impulse percent increase ( I sp Inc.), current ( I ), PAGE 9 9 total thermal heating source ( ), total electrical power ( P W ), and fraction of the total electrical power converted into thermal heating source ( /P W ). ........ 145 PAGE 10 10 LIST OF FIGURES Figure page 2 1 Cold gas thruster from Moog, model 58 125A ................................ .................. 30 2 2 JPL Hydrazine milli Newton Thrusters (Hm NT), a monopropelant technology ... 32 2 3 MIT micro biprop ellant engine concept ................................ ............................... 34 2 4 FMMR thruster design ................................ ................................ ........................ 35 2 5 FMMR thruster ................................ ................................ ................................ .... 37 2 6 VAT thruster ................................ ................................ ................................ ....... 38 2 7 Electro thermal thrusters ................................ ................................ .................... 39 2 8 Schematic of the micro plasma thruster (MPT) ................................ .................. 40 2 9 Ion thrusters ................................ ................................ ................................ ........ 42 2 10 Hall thrusters ................................ ................................ ................................ ...... 45 2 11 Colloid thruster ................................ ................................ ................................ ... 47 2 12 FEEP 5 thruster by Centrospazio/Alta ................................ ................................ 48 2 13 PPT thruster ................................ ................................ ................................ ....... 49 3 1 Voltage Current Characteristics of a DC gas discharge. ................................ .... 54 3 2 Qualitative characteristics of DC glow discharge and electric potential distribution at low pressures (0.1 10 Torr) and room temperature ................... 58 3 3 Paschen cu rves for argon at different values of sec ................................ ........... 64 4 1 Gas flow regimes depending on Kn ................................ ................................ .... 94 4 2 View Factors ................................ ................................ ................................ ...... 98 5 1 Schematic of the new MIG code with rarefied gas and ionized gas modules ... 103 5 2 Generic rectangular parallelepiped element ................................ ............... 11 2 6 1 Order of accuracy for Scharfetter Gummel flux discretization scheme for Case 1 and Case 2 ................................ ................................ .......................... 121 6 2 Comparison between numerical model and experiment ................................ ... 124 PAGE 11 11 6 3 Plane Poiseuille flow in a (L=10, H=1) channel with a (Re = 100) using a (2111) nodes mesh ................................ ................................ ......................... 128 7 1 Geometry and design of the RGEJ thruster using several stack slots similarly to the FMMR thruster ................................ ................................ ........................ 132 7 2 Domain region of the RGEJ numerically simulated ................................ .......... 133 7 3 Regions of applied ................................ ................................ ...................... 134 7 4 Graphical representation of different parameters ................................ ............. 137 7 5 Comparisons between case Q2 base case and case Q4 600 ................................ ................................ ................................ .... 140 7 6 Comparison of base case case Q2 600 and case Q4 600 ............................. 141 7 7 Parameters of the RGEJ thrust er for the given voltage regime ........................ 146 7 8 Left plots show a comparison of the rarefied gas results between the base case (on the bottom half) and case 750V (on the top half) ............................... 147 7 9 Shear stress at the wall for all cases in Table 7 2 ................................ ........... 150 7 10 Comparison of plasma discharge results for case 450V and case 750V (on the left), and comparison of plasma discharge centerline results for all cases (on the right) ................................ ................................ ................................ ..... 151 7 11 Comparison of plasma discharge results for case 450V and case 750V (on the left), and comparison of plasma discharge centerline resul ts for all cases (on the right) ................................ ................................ ................................ ..... 152 7 12 Case 750V ................................ ................................ ................................ ........ 155 7 13 Parameters of the RGEJ thrust er for the given voltage regime ........................ 157 7 14 Conduction heat loss analysis for Case 750V assuming 1 cm of width, with the temperature distribution given along the wall as T =f(x). ............................. 159 7 15 Radiation heat loss analysis for Case 750V assuming 1 cm of width, with the temperature distribution of given along the wall as T =f(x). ............................... 160 PAGE 12 12 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for th e Degree of Doctor of Philosophy THE RAREFIED GAS ELECTRO JET (RGEJ) MICRO THRUSTER FOR PROPULSION OF SMALL SATELLITES By Ariel Blanco Dec ember 201 7 Chair: Subrata Roy Major: Aerospace Engineering Small satellites have multiple applications such as commercial, military, and space science missions. F urther development of micro propulsion concepts is needed. In this dissertation, a micro thrusters for small satellites called the rarefied gas electro j et (RG EJ) is investigated numerically. The RGEJ design aims to increase the efficiency of the energy exchange from the elec trical source to the propellant avoiding the high thermal dissipation problem associated with other technologies. In order to simula te the RG EJ, an existing and previously described in literature, in house modular Multi scale Ionized Gas (MIG) flow solver platform is used. The rarefied gas is modeled using density based compressible flow equations with rarefied boundary conditions and the ionized gas is modeled using local mean energy approximation (LMEA). Benchmarking and validation of the rarefied gas and ionized gas modules are performed, respectively. Several thruster cases are studied with different discharge PAGE 13 13 parameters such as mass flow rate, thrust, power consumption, and specifi c impulse. Finally, a thermal analysis is performed to investigate the thermal losses. For the highest voltage tested (750 V) operating with argon propellant at plenum pressures of 1 Torr the RGEJ thruster requires a total electrical power of 406 mW to heat the flow Under the s e operation conditions, the RGEJ produces a thrust of 1.843 mN and a specific impulse (I sp ) of 60.7 s with a thrust effectiveness of 1240 ( N/W) This I sp is a 37.6% improvement over the I sp of a cold gas thruster operating with the same geometry, propellant, and plenum pressure. The heat loss due to conduction and radiation is expected to be less than 7.5% of the total electrical power. The I sp of the RGEJ thruster was found to be 35% higher than the competing technology called the free molecule micro resistojet (FMMR) Base on the thermal analysis, t he RGEJ is predicted to have significantly fewer heat losses than FMMR. a numerical simulation of a thruster design using p lasma aided technology to replace the energy transfer mechanism of the FMMR thruster and decrease its inherent heat losses h as never been attempted before. PAGE 14 14 CHAPTER 1 INTRODUCTION An artificial, or man made, satellite is a machine that orbits the Earth or another body in space. Their purpose varies greatly; some take pictures and map terrestrial data to help scientist accurately predict the weather patterns tracking hurricanes, deforestation, global warming, etc., whil e others look outward and take pictures and collect data from extraterrestrial objects like stars black holes, dark matter, and search for planets in faraway galaxies for a better understanding of our universe. Other satellites focus primarily on communic ations and the Global Positioning System (GPS), which is composed of a group of more than 20 satellites that allow the determination of the exact location of a GPS receiver. Satellites orbit the Earth at different heights, speeds, and paths. Orbits can be c ategorize d in terms of their centric, altitude, inclination, eccentricity, and synchronicity classifications. A geocentric orbit is an orbit around Earth. with an altitude range of 16 0 to 2000 km, 2000 to 35786 k m, and >35786 km, respectively. The eccentricity of an orbit determines whether it is a clo sed (periodic) orbit or an open (escape) orbit. Circular orbits, with an eccentricity of zero, and elliptical orbits with an eccentricity greater than zero and less than one, are closed orbits. P arabolic and hyperbolic orbits are open while radial orbits can be either. A synchronous orbit has a period that is a rational multiple of the average rotational period of the body being orbited and this orbit is in the same direction as the rotation of the body being orbited. An example of a synchronous orbit is a Geosynchronous Orbit (GEO) with 1:1 ratio of the satellite period to the average rotational period of the body PAGE 15 15 being orbited. A circular Geosynchronous Orbit (GEO) is a geocentric orbit with an altitude of 35786 km and its period equals one sidereal day (23 hours, 56 minutes, 4.091 seconds) The orbital velocity necessary to achieve GEO is approximately 30 7 0 m/s. If a GEO has an inclination of zero degrees with respect to the equatorial plane then it is a geostationary orbit. The two most common types of orbits are geostationary and polar. These denominations refer to the inclination of an orbit with respect to the equatorial plane. In a geostationary orbit, a satellite travels from west to east over the equator in the same direction and rate as the Earth orbit, the satellite travels in a north south direction from pole to pole and it can scan the entire globe one strip at a time as the Earth spins. Once artificial satellites are launched and placed i n their desired nominal orbit, they require some form of attitude control and orbit station keeping to change or adjust their orientation and to compensate for the effects of drag from the thin atmosphere. A thruster can adjust the orientation of a satell ite or compensate for the drag effects by generating an acceleration in the desired direction. This acceleration is generated by expelling the propellant mass at a high velocity, imparting momentum to the spacecraft. Several choices of propellant/energy source ( such as chemical, electric, and nuclear ) are available for powering propulsion devices. These categories can be further subdivided. For example, c hemical propulsion includes cold gas, monopropellant, and bipropellant thrusters. Cold gas thrusters produce thrust using an inert gas and their thermal source is the thermal energy stored in the heat capacity of the gas. PAGE 16 16 Monopropellant and bipropellant th rusters depend on a chemical reaction. For monopropellant thruster s, the energy needed to propel the satellite is stored within the chemical bonds of the propellant molecules and a catalytic reaction is used to release the energy. Bipropellant engines use two types of propellants : a fuel and an oxidizer. Electric propulsion can be subdivided into three categories: electrothermal thrusters ( e.g. resistojets) that heat the propellant using electrical energy, electrostatic thrusters (e.g. I on and Hall thruster s) that accelerate the propellant using the Coulomb force, and electromagnetic thrusters ( e.g. Pulse Plasma and Magnetoplasmadynamic thrusters ) that accelerate the propellant using the Lorentz force. Nuclear powered devices use nuclear energy to generate e lectricity and thermal energy using a reactor, or kinetic energy by means of a nuclear pulse. However, nuclear powered devices are avoided for LEO missions and may be suitable for future interstellar missions Most satellites have chemical thrusters (usual ly hydrazine) for station keeping [1] Some have used electric propulsion for north south station keeping and orbit raising, while interplanetary missions typically use chemical thrusters with the exception of a few missions whe re ele ctric propulsion (primarily ion and hall effect thrusters) has been used with great success [2] The purpose of a propulsion system is to change the momentum of a spacecraft; the amount of change in momentum is called the impulse. The velocity change caused by an impuls e to the vehicle is termed the D elta V ( V ). The amount of impulse that can be obt ained from a fixed amount of reaction mass (propellant mass) is determined by the specific impulse (I sp ). The I sp is a key design parameter in space propulsion. It is defined as the change in momentum per unit of propellant weight spent. The higher the I sp the lesser the mass of the propellant PAGE 17 17 required to achieve a certain velocity increment. The I sp is the metric of comparison between different propulsion technologies and depends strongly on the exit velocity of the propellant for each different propulsio n device. The cost of a satellite is proportional to its mass as shown in the following table by Barnhart et al. [3] Table 1 1 Mass classification of Satellites [3] Cost and weight reductions can be achieved by partitioning the function of a single large satellite into a cluster of small ones, each with a mass less than 100 kg; this concept is called formation flying Formation flying can be used for multiple applicat ion such as to make 3 D views of hurricanes by having multiple satellites orbiting on the same path, separated by a specific time interval with different viewing angles [4] Small satellites have several advantages such as sub stantial reduction in the overall life cycle cost by making satellites less costly to construct and reduced launch costs. Furthermore, the utilization of many smaller satellites to construct a constellation allows for a graceful degradation of the system c apability as individual satellites are lost. PAGE 18 18 In recent years, major satellite manufacturers have presented development programs for small multi mission satellites designed to operate on a Low Earth orbit to reduce the mission cost [5] The size and power consumption of optical and radar instruments scale with orbital altitude for a given instrument performance [5] Therefore, a low operational altitude using small and cheap commercial satellites gre atly reduces the mission cost. Unless a suitable station keeping program is adopted, the large atmospheric drag forces present at LEO altitudes can result in a severe perturbation of the orbital geometry and rapid decay of the orbit [5] For example, the 6 kg Dove 1 spacecraft was released in a 250 km orbit and reentered the atmosphere after 6 days [6] A practical micro propulsion system is the main element for the development of a successful smal l spacecraft, which must perform a variety of mission options such as attitude control, station maintenance, altitude raising, plane changes, and deorbiting [4] Since the field of micro propulsion is still in its infancy, fur ther development of micro propulsion concepts is needed. The Rarefied Gas Electro Jet (RGEJ) Design and Motivation As explained by Micci et al. [4] the following issues are of great importance for micro propulsion systems: material compatibility between the propellant and the surface, contamination problems from propellant ablation and vaporization, valve leakage, passage clogging that could result in a single point failure, system reliability and durability issues, manufact uring complexity, and integration complexity. Thrusters for small satellites must produce minimum impulse bits (N s), determined by attitude control requirements [4] An impulse bit is the single firing impulse imparted to a satellite PAGE 19 19 by a thruster and reflects the level of precision of the propulsion system Small impulse bits produce larger time intervals between thruster firings, which reduce propellant consumption and allow for longer duration quiescent spacecraft opera tion, enabling unperturbed scientific measurements [4] In contrast, the thrust requirements for slew maneuvers extend into the (mN) range very large when compared to the impulse bits requirements for attitude control [4] Propulsion systems for small satellites must overcome these issues while being lightweight, compact, low power, efficient, and inexpensive. Cold gas thrusters have some desirable advantages, such as low complexity, small impul se bit (~ 0.1 mN s) and no spacecraft contamination problems when a benig n propellant is used (e.g., N 2 ) [4] However, valve leakage for light propellants and low specific impulse for heavier propellants limit cold gas application as primary propulsion systems unless the ( V ) requ irements are less than ~100 m/s [4] Leakage in cold gas thrusters is the result of contaminants on the thruster valve seat, low propellant viscosity, a nd h igh pressure propellant storage [4] The I sp of these thrusters is typically low unless very light gases are used such as hydrogen or helium, which produce experimentally measurable specific impulses of 272 and 165 s, resp ectively [4] Hydrogen and helium are commonly avoided due to storage problems. Nitrogen is the most frequently used cold gas thruster propellant with an I sp of 73 s due to reasonable storage density, performance, and lack of contamination concerns; while argon produces an I sp of 52 s [4] In order to increase the performance of cold gas thrusters, while keeping their advantageous characteristics, other designs such as the free molecule micro resi stojet PAGE 20 20 (FMMR) have been proposed [4] FMMR consists of a thin film heating element at constant temperature and a long (~ 1 cm), narrow (1 slot instead of a small nozzle to avoid catastrophically plugging the throat with contaminants [4] FMMR operates with low stagnation pressures (50 to 500 Pa) and low exit Knudsen number ( Kn ~1) [4] Due to the low pressures operating conditions of the FMMR direct interaction with the heating element heats the propellant molecules increasing their kinetic energy before exiting the expansion slots. Gas heating occurs primarily by conduction, as in termolecular collisions are negligible [7] FMMR operating on water propellant with a heating element temperature of 600 K, and propellant storage pressure of ( ~ 100 Pa) can produce an I sp of 68 s [4] If argon propellant i s used instead, the I sp is 45 s [4] These I sp values may not look as an improvement, but other micro propulsion technologies require heavy storage tanks to store the propellant. A nitrogen propellant stored at 20 MPa and 300 K in a spherical titanium tank with a safety factor of 2.0 requires a tank t o propellant mass ratio (M t /M p ) of 1.0, decreasing the effective specific impulse (I sp eff ) of a typical cold gas thruster to half the value of the intrinsic specific impulse (I sp ) [4] As described by Micci et al. [4] t he I sp eff is the effective specific impulse after taking into account the extra mass associated with the minimum operating pressure of the thruster the propellant loss due to valve leakage, and the storage tanks. D ue to the unusually low plenum pressures (50 500 Pa) used in the FMMR, some propellant s can be stored as a solid or a liquid (e.g., H 2 O) For example, t he vapor pressure of ice is ~ 50 Pa at 245 K and if a higher plenum pressure is needed to produce higher thrust, waste heat from the spacecraft could be used to increase the vapor pressure in the storage tank [4] Since the FMMR can operate at the PAGE 21 21 vapor pressure of ice, the propellant could be vaporized on demand to generate thrust The low pressure requirement allows the use of very light storage tanks and decreases valve leakage which depends strongly on the pressure difference between the tank and the ambient and the molecular mass of the propellant [4] Therefore, the I sp eff of FMMR is nearly equal to its I sp Unfortunately, FMMR suffers from significant thermal energy losses due to thermal dissipation inherent in the working principle [8] An estimate of the heat loss per heating element by Micci and Ketsdever [4] predicts 2 0 mW due to radiation and 15 0 m W due to conduction through the pedestal at 300 K. FMMR operating with 40 expansion slots cm) requires 6 000 t o 8 000 m W to heat the propellant gas to obtain the expected performance [4] In outer space, the estimated heat loss of this device would be 6 8 00 m W, which is about the same amount needed to heat the propellant. The total powe r required to operate FMMR would be 12800 to 14800 mW. In contrast, the proposed rarefied gas electro jet ( RGEJ [9] ) micro thruster aims to improve cold gas thruster technology by increasing the I sp while retaining the advantageous characteristics of both cold gas thrusters and FMMR without the excessive heat loss associated with the FMMR design. A different technology called the micro plasma thruster (MPT [10] ), developed t o improve the performance of cold gas thrusters, consists of a cylindrical geometry comprising a constant area flow section combined with a diverging exit nozzle with dimensions of the (~ 100 m) [11] The MPT operates in the abnormal glow regime [11] A direct current micro discharge is used to add thermal energy to the gas in order to increase the performance of cold gas thrusters. The MPT operates in the slip flow PAGE 22 22 regime with ( Kn ~ 0.01 ) at the inlet and ( Kn ~ 0.08 ) at the exit plane and requires a minimum operating pressure of 100 Torr, two orders of magnitude higher than the RG EJ [11] Although the Kn numbers in the MPT between the inlet and the exit plane are compa rable to RG EJ, the expansion and heating of the flow in the MPT occurs predominantly in the diverging nozzle This explains why there is only a 2 3% thrust difference between cases with no slip and slip boundary conditions in the MPT, since m ost of the device (the constant area section ) experiences low Kn numbers (0.01 < Kn < 0.03) [11] Since the device has such small volume (V) to surface (S) ratio (V/S = 39.2 ) the thermal losses due to conduction to the walls are expected to be significa nt. For comparison, the RG EJ has a (V/S = 1.5 mm) and operates in the slip flow regime with ( Kn ~ 0.01 ) at the inlet and ( Kn ~ 0.1 ) at the exit plane, with more than 50% of the device experiencing ( Kn > 0.05). Using no slip boundary conditions would produce a drastical ly different solution for the RG EJ. The sp is 74 s for the 750 V case, an improvement over cold gas thrusters with the same geometry and operating conditions by a factor of (~ 1 .5) [11] However, thrust effectiveness of 50 is very low due to thermal losses through the isothermal walls at 300 K [10] [11] The thrust effectiveness is defined in this study as the difference between the thrust with and without plasma aided technology divided by the total electrical power consumed. Other designs improving the I sp of cold gas thrusters using plasma aided technologies are the radio frequency electrothermal thruster (RFET [12] ), and the microwave electrothermal thruster (MET [13] ). The RFET consist of a cylindrical geometry, 18 mm in length and 4.2 mm in diameter, compose d of alumina with copper PAGE 23 23 electrodes operating at a frequency of 13.56 MHz and an electric potential difference of 240 V. For an argon gas plenum pressure of 1.5 (Torr), the power consumption is 10 W, causing an ionization degree of 0.44% and a maximum elec tron number density of 210 18 (m 3 ) [12] The predicted mass flow rate and thrust are 100 SCCM and 2.619 mN [14] The MET has a cylindrical dielectric chamber, 10 mm long and 1.5 mm in diameter, wit h a metal rod antenna on the axis to produce microwave signals at 4 GHz that generate the plasma and heat the gas. The heated argon gas in the chamber at high pressure (10 50 kPa) is expanded through a de Laval nozzle. An experiment performed with 31 kPa p lenum pressure produced a 60 SCCM mass flow rate and a 1.4 (mN) thrust. Both of these thrusters operate with higher power consumption than the target power budget of RGEJ (< 5 W), but they produce higher specific impulse ( I sp > 80 s) than the RGEJ and the MPT. These designs help illustrate the many different examples of plasma aided technology thrusters currently under investigation. The RGEJ design aims to increase the efficiency of the energy exchange from the electrical source to the propellant avoiding the high thermal dissipation problem associated with by heating the gas directly instead of using heater chips while keeping the advantageous character istics of FMMR The RG EJ concept involves localized embedding of electrodes with a DC or an RF applied potential difference along the dielectric surface of its walls to produce a glow discharge plasma [8] The charged particles in the plasma are ac celerated by the electric field, heating the propellant. Since the deposition of energy is localized via the electrodes, the energy budget should reduce, increasing the thruster efficiency. The presence of charged particles implies that both electric and m agnetic fields c an be used PAGE 24 24 to produce thrust. In this study, magnetic fields were not applied and the thruster operates similarly to an electrothermal thruster, using a DC discharge to heat the gas. RG EJ operates in the slip flow regime, which decreases vi scous losses and heat transfer to the wall. The low minimum required operating pressure (1 Torr), a hundredth (100 Torr) is selected by design to match the vapor pressure at nominal storage temperatures of many liquid or so lid propellants (e.g. H 2 O). The goal is to develop eventually a phase change thruster concept to avoid heavy storage tank and valve leakage problems, vaporizing the propella nt on demand to generate thrust similarly to FMMR [4] RG EJ is predicted to have higher specific impulse than cold gas thrusters as well as significant advantages in size, weight, and power (SWaP) over other micro propulsion devices despite its design simplicity, not requiring any neutralizers, heating eleme nts, or a heavy fuel storage tank if a liquid or solid propellant is used [8] In t his dissertation a finite element and a finite difference based numerical analysis of the rarefied gas and the ionized gas in the FMEJ micro th ruster is implemented to calculate the performance parameters of the device. The cases presented, display a non optimized design with a rectangular channel geometry and adiabatic conditions at the wall with low plenum pressure (1 Torr) operating at differe nt DC voltages to study the effect of power deposition from the plasma to the gas The main difference between RG EJ and FMMR consist of the power deposition method. No heating elemen t is required in RG EJ. Instead, the gas is heated by the volumetric therma l heating source increasing the efficiency of the energy transfer. There are many differences between the RG EJ and the MPT. The MPT has higher plenum pressure PAGE 25 25 (100 Torr) which requires a complex argon chemistry model, assumes isothermal walls, has a cylindrical geometry with an exit expansion nozzle, and a smaller volume/surface ratio, which increases particle loss to the walls a s well as thermal losses. The RG EJ has lower plenum pressure (1 Torr) to develop eventually a phase change thruster concept similar to FMMR, which has the added advantage of combining microelectromechanical systems (MEMS) fabrication techniques with a simple, lightweight design d ue to the low plenum pressures, which causes a smaller pressure difference between the thruster inl et and the ambient than the MPT RGEJ is design ed with insulated walls, and has a rectangular slot geometry to prevent catastrophic plugging of the throat, and does not have an expansion nozzle since for many micro propulsion devices the Reynolds number (R e) is so low (~ 100) that the viscous losses counteract any improvement in I sp It is important to emphasize that simulating the thruster with adiabatic conditions requires a more robust iterative procedure between the loosely coupled modules of the code t han if isothermal walls are used. Although argon was u sed as the propellant for the RG EJ for simplicity of the numerical simulation, other propellants (e.g. H 2 O or NH 3 ) may be preferable since they can be stored in solid or liquid form to reduce storage ta nk weight and valve leakage. Performance Parameters In order to evaluate the performance of the micro thrusters, several parameters must be calculated and compared. The mass flow rate ( ) is a measurement of how much mass per time flows through the thruster. This parameter can be used to calculate the total propellant mass necessary for a certain amount of operating time of the thruster. The thrust ( F Thrust ) is the propulsive force of a propulsion system and in space, it dictates the achievable acceleration of a space vehicle given the mass of that vehicle PAGE 26 26 As previously mentioned, the amount of spacecraft momentum change caused by a propulsion system operation is called the impulse. The minimum finite change of momentum that can be delivered by a thruster (the minimum impulse) is defined as the minimum impulse bit ( I bit ). The velocity change caused by an impulse to the vehi cle is termed the Delta V), which is a measurement of the effort require d for making an orbital maneuver These values depend on the valve speed and the mass of the spacecraft in question, respectively. The minimum impulse bit determines the ability of a thruster to do minute momentum adjustments for attitude and orbit control, while the Delta V produced by a propulsion system determines its range of application for a given satellite. The specific impulse ( I sp ) is defined as the change in momentum per unit of propellant weight ( using Earth's frame of reference ) spent and it is a measure of the efficiency of a propulsive system in terms of how much mass is need ed to produce a certain amount of thrust The total efficiency ( ) of an electrically powered thruster is defined as the jet power divided by the total electrical power into the thruster. It indicates how efficiently the power is converted to kinetic energy of the propellant in an electric propulsion system. The thrust effectiveness ( ) is the ratio of thrust vs. total electrical power into the thruster. The following equations are used to calculate thes e parameters at the exit plane of the thruster [11] [2] ( 1 1 ) PAGE 27 27 where the variables u, and p are the density, velocity component in the x direction, and pressure of the gas at the exit plane of the thruster while P indicates the power H and W are the height and width of the thruster, respectively, use to calculate the area ( A ) where W is assum ed to be 1.0 cm. The gravitational constant g 0 is 9.81 m /s 2 The dissertation has the following structure: Chapter 2 provides a review of m icro propulsion d evices Several micro propulsion technologies are described and their pe rformance characteristics are provided to compare their advantages and disadvantages for different applications Chapter 3 provides a background on gas discharge physics. The voltage current characteristics of a DC discharge are examined, the glow discharge mode is explained, and a methodology to estimate sputtering is described ( 1 2 ) ( 1 3 ) ( 1 4 ) ( 1 5 ) ( 1 6 ) ( 1 7 ) PAGE 28 28 Chapter 4 explains the governing equations use d in each numerical module and the thermal analysis. Chapter 5 explains the numerical models in each module and the thermal analysis. Chapter 6 shows the validation/benchma rking of the numerical modules. Chapter 7 shows and explains the results obtained by the numerical simulations and estimates the heat loss through conduction and radiation from the thrusters to justify the assumption of adiabati c walls. Chapter 8 gives a summary of the results and contributions to the field of science; additionally, it proposes future directions to continue improving the RGEJ. PAGE 29 29 CHAPTER 2 MICROPROPULSION TECHNOLOGIES Introduction to Micro propulsion Devices Some of the competing micro propulsion technologies were mentioned in the previous chapter along with their strength and weaknesses. In order to understand the different designs of micro propulsion systems and their limitations, a survey of existing techno logies is presented in this chapter. This survey is not intended to be comprehensive; instead, it aims to provide enough information to appreciate the complications that come with the miniaturization of existing technologies. For larger satellites, the av ailable conventional technologies have been used reliably. For smaller satellite, successful implementation of conventional technologies become increasingly more difficult due to size and weight limitations In the last two decades, several prototypes have been designed and tested with different degrees of success. Extensive reviews and similar surveys of the available thruster designs and their performances are available in litera ture [4] [8] [11] [15] [16] [17] [18] [19] In this survey, the focus is on chemical and electric propulsion categories of thruster designs and their subtypes based on the working principle, design details, and propellant types. Chemical Propulsion Devices for Small Satellites Chemical propulsion o ptions are usually characterized by limited specific impulse [16] Their thrust levels may exceed several Newton for micro propulsion designs. Chemical propulsion options are of particular interest when rapid maneuvering is requ ired, such as rapid plane change or orbit change, or proximity operations near other PAGE 30 30 spacecraft, such as for inspection purposes, were sufficient agility is required for collision avoidance [16] Cold Gas Thrusters The cold gas thruster is one of the most simple propulsion systems. A gas from a high pressure supply tank that expands through a valve and a nozzle produces thrust Cold gas thrusters offer low impulse bits and can operate with contamination free, non toxic propellan t (e.g., N 2 ) [4] Thrust values range from 4.5 mN to 100 N and some s thruster model 58 125A can weight as little as 9 grams and are as small as ( 12 35 mm ) [15] [19] Figure 2 1 Cold gas thruster from Moog, model 58 125A [15] This cold gas thruster uses N 2 for propellant, has a thrust of 4.4 mN, a mass of 9 g, a size of (11.9 34.7 mm), operational pressure from 0 to 50 psia, an I sp of 65 s, an impulse bit of 100 N s, and it requires 10 W to operate its valve. Most of the cold thruster designs in existen ce were not designed for micro spacecraft applications. The power leve ls required for valve actuation (~ 10 W for the Moog 58 125A model) exce [15] Cold gas thrusters have low specific impulse, I sp < 7 5 s for N 2 propellant [4] leading to large propellant mass fractions. The larger propellant mass fraction increases the propellant system mass since a heavy tank is required to contain the high pressure PAGE 31 3 1 gas supply ( ~ 3 500 psia) [4] The large storage pressure also creates a risk of propellant leakage across valve seats. Several approaches have been considered to avoid some of the problems with the cold gas thruster and are elaborated in det ail by Mueller et al. [15] The following table lists the typical cold gas performances, depending on propellant type using data from Micci et al. [4] Table 2 1 Cold gas propellant performance Assuming propellant is at 25 C and the flow expansion is to zero pressure in the case of the theoretical value [4] Propellant Molecular weight (kg/kmol) Density (3500 psia, 0 C) (g/cm 3 ) Isp (s) Theoretical Measured Hydrogen 2.0 0.02 296 272 Helium 4.0 0.04 179 165 Neon 20.4 0.19 82 75 Nitrogen 28.0 0.28 80 73 Argon 39.9 0.44 57 52 Krypton 83.8 1.08 39 37 Xenon 131.3 2.74 31 28 Freon 12 121 46 37 Freon 14 88 0.96 55 45 Methane 16 0.19 114 105 Ammonia 17 Liquid 105 96 Nitrous oxide 44 67 61 Carbon dioxide 44 Liquid 67 61 Neither hydrogen nor helium is commonly used, due to the requirement of a large heavy tank caused by the low gas densities and leakage problems due to the low molecular weight of these gases [4] Mono Propellant Thrusters Monopropellant thruster s use only one working propellant T he propellant undergoes an exothermic chemical reaction that provides the power to generate thrust producing a typical speci fic impulse of ~220 s [4] If hydrazine is used as a propellant it decompose s catalytically in a catalyst bed consisting of iridium coated alumina pellets. PAGE 32 32 The decomposition is highly exothermic and produces a 1000 C gas tha t is a mixture of nitrogen, hydrogen, and ammonia [4] Ammonia decomposes into hydrogen and nitrogen and the degree of decomposition depends on the feed pressure, catalyst type, and geometry, among other factors. The smallest hydrazine thrusters have a thrust in the (0.9 4.45 N) range [4] Hydrazine thrusters are suitable as the primary propulsion system for intermediate to low ( V ) maneuvers ~1000 m/s and some designs, such as the Primex MR 103 C/D model, can produce small impulse bits (50 100 N s) for fine attitude control in small satellites [4] A hydrazine milli Newton thruster (HmNT) at JPL may fit the mass and volume envelope of a CubeSat for small to intermediate ( V ) applications [15] [20] Figure 2 2 JPL Hydrazine milli Newton Thrusters (HmNT), a monopropelant technology [20] It has a thrust of 129 mN, a mass of 40 g, a volume of 8 cc, an I sp of 150 s, and it requires 8 W (open)/ 1 W (hold) to operate its valve. It can produce minimum impulse bits of 50 N s [4], [15] Two disadvantage s of hydrazine engines are its toxicity and flammability, which results in complexity and cost of ground handling and propellant loading procedures [4] Another technology availa ble uses hydrogen peroxide, purified to 90% or higher concentration. Hydrogen peroxide is self decomposing at high temperatures or when a catalyst is presented. Hydrogen peroxide thrusters may fill a similar function as the HmNT hydrazine thrusters. The ma in advantage of h ydrogen peroxide is tha t it is non PAGE 33 33 toxic unlike hydrazine [15] The disadvantage is that hydrogen peroxide slowly decomposes when it is heated or exposed to a catalyst, and almost any organic residue can act as a catalyst [15] Other monopropellant technologies are explained by Soni [8] and Micci et al. [4] The main problems with the monopropellant thrusters for small satel lite applications seem to be the toxicity of the propellant the storage of the propellant, and the heat loss through the thruster structure along with degradation of the catalyst bed depending on the propellant and design used. Bipropellant Thrusters Bipr opellant thrusters require two types of propellants, a fuel and an oxidizer, and are pressure regulated to maintain the mixture ratio. The added complexity and dry mass associated with the design is a tradeoff to obtain higher specific impulse (~ 300 s) th an the one in monopropellant thrusters [4] [15] Due to the additional components in bipropellant thrusters, the monopropellant thrusters outperform them in terms of wet propulsion mass below V s of approximately 500 1000 m/s, depending on spacecraft mass [15] The added complexity increases the cost of bipropellant systems, usually making them more expensive than monopropellant thrusters [4] A good example of a miniaturized bipropellant thruster is the MEMS micro fabricated bipropellant rocket engine achiev ed by London et al. [21] The thruster chip consists of an integrated thrust chamber, nozzle, turbine p umps, and inlet valves [22] Some of the challenges associated with miniaturizing bipropellant engines include managing the losses in combustion efficiency due to reduced mixing and vaporization and dealing with an increase in heat losses into the engine structure resulting in thermal control issues [4] PAGE 34 34 Figure 2 3 MIT micro bipropellant engine concept [22] Testing has been performed w ith gaseous oxygen and ethanol, the thrust obtained is 2.7 N, with an Isp of 300 s. The entire assembly of the thruster chip is about 20 5 mm and consists of an integrated thrust chamber, nozzle, turbine pumps, and inlet valves [15] Electric Propulsion Devices for Small Satellites T he acceleration method producing the thrust is typically used to describe the electric thruster These methods are separ a ted into three categories [2] : E lectrothermal E lectrostatic E lectromagnetic Electrothermal Thrusters Electro thermal thrusters are a class of electric propulsion devices that convert electrical energy into thermal energy of the propellant. The n, the propellant is expanded through a nozzl e or through a long thin channel, to produce thrust. Elevating the propellant temperature increases the specific impulse. PAGE 35 35 Free m olecule m icro r esistojet (FMMR) In this concept, propellant from a plenum at very low stagnation pressures (50 to 500 Pa) flows through very thin micro slots (1 such that the flow is in the free molecular regime [4] C ollisions over the characteristic thruster dimensions only take place with the hot wall surfaces, not among the molecules of th e propellant flow due to the degree of rarefication of the gas Since there is no momentum exchange between the flow particles, all the particles colliding with the heater wall surface will acquire the kinetic energy equivalent to the heater wall temperatu re [15] The increase in propellant kinetic energy increases the I sp of the thruster over similar cold gas thrusters. There are different designs for the FMMR thruster they are all based on the same concept but the geometry is adapted to improve heat transfer to the gas and minimize heat losses. In the design shown in Micci et al. [4] the FMMR geometry is as follow s A B Figure 2 4 FMMR thruster design [4] A) Schematic cross section showing the heating element arrangement with the expansion slot, and B) multi slot (w = 100 m) configuration with a 0.5 0.5 cm cross sectional area. PAGE 36 36 where the design requirement is to arra nge t he thin film heater at a given stagnation temperature to be the last surface contacted by a prop ellant molecule before it is ejected from the expansion slot. Using this design, with ten different 100 m 8 mm slots and expansion angle of 54.74 degree s the FMMR is able to produce 0.25 mN of thrust at an I sp of 45 s for argon propellant and thin film heater temperature of 600 K The same design working under the same parameters using water as a propellant is able to produce an I sp of 68 s. The effect of propellant gas (molecular mass, MM) on the I sp has been investigated and shown in Micci et al. [4] for thin film heater temperature of 600 K, where the Gases with smaller MM produce higher I sp (e.g. h elium, I sp ~ 140 s [4] ), a similar effect is observed in cold gas thrusters H owever propellants that can be stored in liquid or solid states at standard temperature and pressure (STP) such as water or ammonia are favo red to decrease propellant tank mass and valve leakage. The heat loss analysis shown previously in chapter 1 is performed using this A different design has the slots etched into the heater chip, with the resistive e lements patterned around the periphery of the slots as shown in Figure 2 5 [23] Experimental testing with a heater chip temperature of 580 K, using water propellant, achieved thrust values and specific impulse of 129 N and 79.2 s, respectively [24] This design has higher I sp at lower heater temperature than the one in Micci et al. [4] However, the main problem associated with the FMMR concept is the low energy conversion e fficiency due to heat loss which increases the power consumption and may present a significant demand on the power system of a small satellite [8] [15] PAGE 37 37 A B Figure 2 5 FMMR thruster [23] A) Scanning electron microscope side image of cleaved expansion slots, and B) heater chip installed in Teflon flight plenum. Low p ower DC a rcjets Arcjet thrusters use an electr ical arc discharge to impart additional energy to the propellant in order to increase the I sp at the cost of increased power consumption. The design typically has the electrodes placed at the no zzle exit to minimize heat loss. I t can operate with gaseous, liquid, or solid propellant. If the device uses gaseous propellant, it functions as a warm gas thruster with the arc discharge heating the flow. If the design uses solid propellant, the arc ablates a small amount of material from the propellant block whic h sometimes is designed to act as the cathode. Subsequently, the vaporized material is accelerated using thermal expansion [8] T he Vacuum Arc Thruster (VAT) is an example of an ablative pulse plasma DC arcjet thruster [15] A h igh voltage between two metal electrodes causes high electric field near surface roughness spikes. The high power density in the spikes leads to vaporization, electric breakdown, and plasma generation in an arc discha rge. The plasma plume expands into a vacuum and produces a thrust pulse Test on a JPL micro Newton thrust stand demonstrated up to 100 N s impulse bits for 100 Hz pulse trains using Cr as the propellant block material [25] T his technology has demonstrated PAGE 38 38 the ability to achieve very high specific impulse of 1000 3000 s [26] [25] The thruster heads weigh ~ 100 g and the power consumption range between 1 100 W with an eff ectiveness of up to 10 N/W [15] [25] A B Figure 2 6 VAT thruster. A) VAT f iring on thrust stand, and B) power processing unit (PPU) configured for CubeSat with mass of 150 g capable of supplying four V AT thruster heads with power using a 12 24 V bus voltage [15] [25] The advantages of the VAT are its very small impulse bits for attitude contro l, its so lid propellant storage and integrated propellant design which allows for a compact, modular propulsion system The negative characteristics of the design are that large V ) applications may not be possible due to the small impulse bits that limit the t hrust and impulse values [15] Rapid maneuvering may not be possible due to low thrust values. The VAT main applications could be precision pointing and attitude control. Micro cavity discharge thrusters The Micro Cavity Discharge Thru ster (MCD) is a type of electro thermal thruster concept that utilizes a micro discharge in a cavity as small as 10 m in diameter to heat the gas as high as 1500 K [15] The gas is heated by an applied voltage ( 40 0 1200 V AC ) to the electrodes at a frequency of 50 150 kHz which causes the gas to partially ionized PAGE 39 39 with a low degree of ionization (<<1%) [15] The following figure shows the thruster concept and etched micro nozzle [15] [27] A B Figure 2 7 Electro thermal thrusters A) S chematic of MCD electro thermal thruster concept, and B) micro nozzle [15] [27] For a n array of 4 cavities, each with 120 m throat and 210 m exit plane diameter, the mass flow rate and thrust values obtained were between 0.99 5.22 mg/s and 0.6 2.7 mN, respectively, operating with plenum pressure s of 120 240 kPa and 0.25 W per cavity [15] These values produce an I sp of 50 60 s. The MCD thrusters operate at higher gas temperatures than conventional res istojets, promising higher I sp Their design has higher thrust and effectiveness than electrospray arrays, yet lower specific impulse values [15] Their possible uses include attitude control, rapid slew maneuvers, proximity operations, and low ( V ) applications [15] Micro p lasma thruster (MPT) The micro plasma thruster (MPT [10] ), developed to improve the performance of cold gas thrusters, consists of a cylindrical geometry comprising a constant area flow section combined with a divergin g exit nozzle with dimensions of (~ 100 m) [10] [11] PAGE 40 40 Figure 2 8 Schematic of the micro plasma thruster (MPT) [11] The geometry consists of an axisymmetric constant area section a diverging respectively. The radiuses of the upstream constant area section and exit section are 50 he i nlet total pressure is 100 Torr and a small but nonzero outlet pressure domain. The temperature was constant at 300 K for the inlet, the outlet, and the solid walls. T he MPT has a mass flow rate of 0.14 mg/s, a thrust of 67.4 N for the cold gas thruster case, and a thrust of 100 N for the 750 (V) and 650 (mW) plasma aided case which produces an I sp of ~ 74 (s) [10] [11] effectiveness and total efficiency equal to 50 ( N/W) and 0.6%, respectively. Electrostatic Thrusters Electrostatic thruster s accelerate the propellant to high exit velocities using mainly Coulomb force. The design uses electrostatic forces to accelerate heavy charged ions towards the thruster exit. Electrostatic thrusters are classified into three main types: gridded ion engines, Hall effect thrusters and colloid thrusters. PAGE 41 41 Miniature i on e ngines In an ion engine, the propellant is ionized in a discharge chamber and ions are extracted from this chamber using an electrostatic grid system. Usually, heavy ine rt gases are used as a propellant (e.g. x enon ) Two main types of ion engines are available: DC electron bombardment (Kaufmann type thruster) and Radio Frequency (RF) ion engine [15] In a DC electron bombardment ion engine a plasma is generated in a hollow cathode. From this plasma, electrons are extracted and injected into the discharge c hamber to ionize the propellant gas by accelerating the electrons towards an anode surface [15] Typically, the design uses permanent magnets to increase the e ngine efficiency by reduc ing electron loss to the walls and increas i ng the electron path through the discharge. In an RF ion engine, an electromagnetic coil is wrapped around a dielectric discharge chamber to generate an oscillating axial magnetic field that generates an oscillating azimuthal electric field [15] Free e lectrons always present in a gas are accelerated with enough speed to cause ionization of the neutral particles which increases the number of free electrons, causing an electron avalanche that forms the plasma. The effectivenes s, the thrust per power, of the DC discharge engines is typically higher than for RF discharges, but no life limiting internal cathode is needed for the second design, making RF engines very desirable for long term operations. Bo th designs require a neutra lizer to avoid spacecraft charging, consisting of another hollow cathode outside the thruster that neutralizes the positively charged ion beam exiting the thruster. PAGE 42 42 In recent years, different miniature ion engines have been developed for formation flying applications of space telescopes, such as the JPL Miniature Xenon Ion Thruster (MiXI [28] ), the NRIT 2.5 Miniature Ion Thruster [29] and the Miniature Radio Frequency Ion Thruster (MRIT [15] ). A B C Figure 2 9 Ion thrusters A) JPL Miniature Xenon Ion thruster (MiXI) [28] B) the NRIT 2.5 miniature ion thruster of the University of Giessen [29] and C) Miniature Radio Frequency Ion Thruster (MRIT) of the Pennsylvania State University [15] Miniature ion engines have multiple beneficial characteristics as space telescopes propulsion systems such as the use of non contaminating propellant, which is of great importance near sensitive optical surfaces, and the ability to modulate the thrust amplitude without inducing jitter [8] The following table shows th e performance parameters of these three designs [15] PAGE 43 43 Table 2 2 Micro ion engine p erformance [15] Parameter Units MiXI NRIT 2.5 MRIT Propellant Xe Xe Ar Thrust (mN) 0.01 1.5 0.05 0.6 0.001 0.06 Specific Impulse (s) 2500 3200 2861 5480 Power (W) 13 50 13 34 Electrical Efficiency (%) >40 4 47 15 Diameter (cm) 3 2.5 2 Mass (g) 200 210 As observe d in the previous table these engines have a high specific impulse, but high er power consumption and mass requirements in comparison to other technologies. In addition, performance deterioration due to grid erosion caused by heavy ion bombardment is a problem encountered in gri dded ion thrusters [8] In order to adapt micro ion engines to small satellites, their size as well as the size of the power processing unit would have to be drastically reduce d Miniature hall thrusters In a Hall thruster, the external cathode, outside the thruster, emits electrons that are accelerated towards the anode located inside an annular channel. The electrons encounter a radial magnetic field near the exit plane of the thruster T he Lorentz force which is caused by the combined electric and magnetic forces on a point charge due to electromagnetic fields, makes the electrons gyrate around the magnetic field lines and drift azimuthally through the annul ar channel causing ionization. The ions are accelerated out of the thruster by the electric field generated by the electric potential difference between anode and cathode. The ion flow exiting the thruster generates the majority of the thrust. T he Larmor radius, the radius of circular motion of a charged particle in the presence of a uniform magnetic field, is larger for ions than for electrons due to the higher inertia of the ions Since the ions are affected less by the magnetic PAGE 44 44 field, the ions move out of the channel without the magnetic field affecting their path The external cathode neutralize s the ions outside of the channel by providing a source of electrons As a result, Hall thrusters are more compact than ion engines for a given ion beam strength [16] The main difference between Hall thrusters and ion engines is that Hall thrusters do not require grid systems. A Hall thruster design can be categorized as either a stationary plasma thruster (SPT) or a thruster with anode layer ( TAL ) design [2] The main difference is that SPT design ha s a wide acceleration zone and uses dielectric insulating walls such as alumina (Al 2 O 3 ) with low sputtering yield and low secondary electron emission coefficient T he TAL design ha s a narrow acc eleration zone and metallic conducting walls with exit plane metallic guard rings biased at cathode potential to reduce the electron loss along the field lines The SPT type was used in the Meteor spacecraft in 1971. Unfortunately, n either of these designs is suitable for miniaturization. The miniaturization of H all thrusters faces many challenges since the reduction in size requires smaller electron Larmor radius (stronger magnetic fields) to avoid excessive electron losses to the wall s and the associ ated decrease in electrical efficiency [15] A micro Hall thruster requires stronger magnets H owever the smaller devices suffer from increased heating of the magnets by the plasma, which could lead to demagnetization and the loss of electron containment causing the electric efficiency of the device to drop drastically. A miniature 50 W Hall thruster conc ept by MIT suffered from this particular problem, having a poor efficiency of 6% [30] PAGE 45 45 A B C Figure 2 10 Hall thrusters A) MIT 50 W Hall Thrus ter, B) 3 cm CHT device, and C) Schematic of 3 cm CHT device [15] A new type of miniature Hall thruster has been developed by Princeton Plasma Physics Laboratory (PPPL) to improve the efficiency observed in previous versions [15] This new design uses a c ylindrical Hall thruster (CHT [31] [32] ) geometry with a largely axial magnetic field, instead of the annular flow channel with a radial magnetic field [33] In th e CHT design, the axial magnetic field is concentrated at the anode end of the channel causing electrons to be mirrored back, away from the anode region [15] The CHT effic iency, ~30% is higher than for other micro Hall thrusters. However, it is significantly less than for conventional H all thrusters, which may suggest performance deterioration as a result of downscaling [2] Table 2 3 shows the performance parameters of several micro Hall thrusters and their range of efficiency [15] The improvement in specific impulse and electrical efficiency of the new design looks promising. The miniature Hall thrusters have some advantages over ion engines since they require a power processing unit capable of producing 300 V instead of 1 kV and fewer power supplies than thei r counterparts do. PAGE 46 46 Table 2 3 Small and micro Hall thruster performances [15] Parameter Units BHT 200 SPT 30 MIT PPPL CHT 2.6 PPPL CHT 3.0 Propellant Xe Xe Xe Xe Xe Thrust (mN) 4 17 5.6 13 1.8 2.5 12 3 6 Specific Impulse (s) 1200 1600 576 1370 865 1100 1650 Power (W) 100 300 99 258 126 50 300 90 185 Electrical Efficiency (%) 20 45 16 34 6 15 32 20 27 Diameter (cm) 2.1 3 0.4 2.6 3 Mass (g) <1 ~1 However, H all thrusters are heavier due to the ir magnet mass which increases the total mass of the propulsion system [15] Colloid al and f ield e mission e lectric p ropulsion (FEEP) thrusters Colloid al thrusters are a type of electrostatic thruster that do es not require a plasma discharge eliminating the efficiency losses associated with the miniaturization of the plasma chamber found in ion engine and hall thruster designs. Colloidal thrusters emit ions or charged droplets directly from a liquid propellant column when a strong electrostatic field is applied dry powder propellants have also been used [16] In an electrospray thruster, the propellant is injected through an emitter. The emitter can be either a sharp needle, a capi llary, or a narrow slit. An electric field is applied between the emitter and an opposing electrode. The sharp emitter tip causes the field to intensify. Due to the balance of surface tension and electrostatic forces at the tip of the emitter, the conducti ve propellant is distorted into a sharp Taylor cone whose shape intensifies the electric field strength even further [15] Depending on the choice of propellant either ions or charged liquid droplets are emitted. Thrusters that work with liquid metal propellants and emit ions are typically refer red to as Field Emission Electric Propulsion (FEEP) [15] Thruster s that use primarily droplet PAGE 47 47 emissions, typically working with glycerol or ionic liquid prop ellants, are refer red to as colloidal thrusters [15] It is important to note that the difference between these two thruster types is not well defined in the case of the thrusters using ionic liquid, m olten salts at room temperature since they can emit either ions or droplets depending on the operating conditions [15] Figure 2 11 (A, B) sho w s the ST 7/LISA colloid propulsion system consisting of four heads per cluster and nine individually machined emitters per head [15] [34] A B Figure 2 11 Colloid thruster A) I ndividual ST 7 colloid thruster emitter, actual thruster is composed of multiple emitters, and B) cluster head (4 per cluster) [15] [34] The f ollowing table show s the operational parameters of the ST 7 colloid thruster [15] [35] The ST 7 colloid thruster cluster is too large for CubeSat applications, but each thruster head could form a stand alone thruster. The biggest miniaturization challenge of colloid Table 2 4 ST 7 c olloid t hruster c luster p erformance c haracteristics [15] [35] Parameter Units Demonstrated Performance Thrust ( N) 5 35.8 Thrust Noise ( N) < 0.01 Specific Impulse (s) 240 Max Power (W) 24.6 (2 W per head) Nominal Power (W) 16 Beam Divergence ( 1/2 angle) ( d egree) <23 Mass (kg) 15 PAGE 48 48 propulsion systems is not the thruster head itself, but the required miniaturization of the dedicated feed and power processing unit (PPU) subsystems T he FEEP 5 created by the University of Pisa and Centrospazio is shown in the following figure, which is a field emission elect ric propulsion thruster design [15] Figure 2 12 FEEP 5 thruster by Centrospazio/Alta [15] One of the main differences between colloid and FEEP thrusters, besides the latter emitting predominantly ions instead of charged liquid droplets, is that FEEP thrusters are not pressure fed [15] Capillary forces supply the propellant to the emitter, simplifying the design of th e thruster by not requiring valves or propellant tanks. The following table shows the performance parameters of several FEEP thrusters [15] FEEP thrusters allow for extremely high specific impulse values, however, at reduced thrust to power ratios, requiring either long transfer times or high power levels Table 2 5 FEEP Performance Characteristics [15] Parameter Units In FEEP 100 GOCE MTA In FEEP 1000 FEEP 5 FEEP 50 Thrust (mN) 0.001 0.1 0.002 0.65 0.001 1 0.04 1.4 I sp ( 10 3 s) 8 12 8 12 8 12 9 9 Power (W) 0.5 10 6 52 2 80 2.7 93 PPU Mass (kg) Incl. Incl. Incl. 1 1.2 Mass (kg) 0.3 3.5 1.5 0.6 1.2 PAGE 49 49 [15] Cesium, t he most common propellant used, is very corrosive and may clog emitter slits. A more suitable propellant is indium, requiring heater power to liquefy it. Elect ro magnetic Thrusters Electro magnetic thrusters accelerate the flow utilizing electrical current and magnetic fields. They accelerate the ionized propellant using the Lorentz force or they accelerate the ionized propellant by the effect of an electromagnetic field where the electric field is not in the direction of the acceleration. Pulse plasma thru sters (PPTs) : In a PPT, a capacitor connected to two electrodes is charged and discharge d [15] T he arc produced between the two electrodes by the capacitor discharge ablates the solid Teflon used as a fuel rod. The ablated material accelerates due to the Lorentz force produced by the cur rent flowing through it and its interaction with the self induced magnetic field The following figure shows the Micro Liquid Pulsed Plasma Thruster (MILIPULT), containi ng a small liquid reservoir that allows a flow to creep between the electrodes [8] [15] [36] A B Figure 2 13 PPT thruster. A) S chematic of a liquid propellant PPT [8] and B) APL MILIPULT liquid micro PPT [15] [36] PPTs offer very small impulse bits, solid propellant storage in the form of T eflon fuel rods, modularity, and proven operation. The following table shows the performance of different designs [15] PAGE 50 50 The power electronics and storage capacitors account for 40 45% of the system mass [15] Miniaturization of PPT technology should focus on reducing the mass of the thruster head as well as supporting power electronics. The thruster itself faces challenges d uring miniaturization for small satellite applications. For example if the dischar ge energy, because of overall system miniaturization, is too small then the carbon neutrals in the plasma arc can return to the fuel rod surface an d result in charring which can cause shorting of the thruster electrodes. Conclusion of Micro propulsion Tech nologies Survey The above survey was presented i n order to provide some documentation on the existing competing technologies and to depict the strengths and weaknesses of each design. Due to the vast number of micro propulsion technologies under investigat ion, only some of the existing designs are presented here and o nly a few of these designs have actually been commissioned and successfully flow n as part of the propulsion system of small satellites [8] Many of these technologies available cause concerns Table 2 6 Micro PPT Performances [15] Parameter Units Dawgstar AFRL MILIPULT Propellant Teflon Teflon Teflon Water Impulse bit s) 60 2 2 25 0.4 0.6 Thrust 61 264 Specific Impulse (s) 266 Mass (g) 4200 500 30 (5 g prop ) 13.5 Power (W) 15.6 36 1 20 0.5 4 Capacitance 1.3 (per module) 0.42 2 6 0.937 Voltage (kV) 2.8 2.45 5.36 0.2 0.7 Efficiency (%) 1.8 PAGE 51 51 related to their reliability, weight, power requirements, and propellant dispensation, or thrust resolution and repeatability [8] The goal of the Rarefied Gas Electro Jet ( RG EJ) thruster presented in this dissertation is to improve cold gas thruster technologies and provide a viable alternative propulsion system for small satellites. PAGE 52 52 CHAPTER 3 GAS DISCHARGE PHYSICS S ome basic plasma physics concepts are explained in this chapter in order to understand the physics used to model numerically the RG EJ thruster Plasma refers to a system of charged particles, which demonstrate collective behavior [37] Irwin Langmuir chose the similar to positive ions and electrons in a plasma [38] In the real world, p lasmas are compos ed of more than two particles They may also contain negative ions as well as metastable atoms or molecules The plasma state is known as the fourth state of matter due to its unique properties that separate it from solids, liquids, and neutral gases. The charged particles that plasma contain s make it electrically conductive allowing it to interact with external electromagnetic fields and to produce its own electromagnetic fields due to space charge effects and charge carrier currents [39] [40] [41] It is estimated that 99% of the matter in the entire visible universe is plasma making it the most prevalent observable state of matter. In a laboratory, plasma is typically produced by the coupling of electromagnetic energy into a neutral gas through a discharge [42] [43] There are many types of dischar ges at low pressure such as glow discharges, capacitive l y coupled plasma (CCP) discharges inductively coupled plasma (ICP) discharges etc. [11] G low discharges use a Direct Current (DC) or a low frequency Radio Frequency (RF) electric field (<100 kHz) applied to a gap between ele ctrodes Capacitive ly coupled plasma (CCP) discharges use high frequency RF electric fields (~10 MHz). I nductively PAGE 53 53 coupled plasma (ICP) discharges are similar to CCP with the difference that electrodes consist of a coil wrapped around the plasma chamber In this study, DC glow discharge s are of particular interest. Glow discharges are generated due to an electrical breakdown of gases by an externally applied potential [43] The mechanism of electrical breakdown of gases is best described in terms of voltage current characteristics of the discharge [8] which is explained further in the following section. Electrical Breakdown of Gases Plasmas are generated when electrons are ejected from neutral atoms when energetic photons or free electrons with sufficient energy collide with the neutral atoms (or molecules) generating ions and more free electrons [8] This effect c ould be achieve d in the laboratory by applying an external electric potential to a gas [43] The mechanism of electric breakdown can be better understood by considering what happens to the current f lowing through a low pressure gas (1 10 Torr) contained in a sealed glass tube with electrodes at either end as the voltage across the electrodes is gradually increased [41] Figure 3 1 shows the voltage current characteristics of a DC gas discharge in a cold cathode tube [44] The discharge has three main regions: dark discharge glow discharge and arc discharg e The voltage current characteristics plot of a DC gas discharge is nonlinear and it depends on the geometry of the electrodes, the gas used, the pressure, and the electrode material. PAGE 54 54 Figure 3 1 Voltage Current Characteristics of a DC gas discharge [44] Dark Discharge In the dark discharge region, the regime between A and E in Figure 3 1 the discharge remains invisible to the human eye except for corona discharges and the breakdown itself. In this regime, the current is less than 10 5 A [43] The dark discharge is subdivided into three distinct sub regions : the background ionization stage, the saturation regime, and the Townsend regime. In the A B sub region, called the background ionization stage, the background radiation generated by cosmic rays and other sources causes a constant degree of ionizat ion of the neutral particles in the gas. T he applied electric field sweeps out the ions and electrons produced by the ionization process Increasing the voltage between electrodes causes an increasing fraction of these charged particles to be sweep toward the electrodes producing a weak electric current. In the B C sub region, called the saturation region, the current remains constant with increasing voltage. T he voltage between electrodes is high enough at point B to PAGE 55 55 sweep all available electrons and ion s away. The current in this region depends linearly on the background radiation source strength and reaches a saturation point at point B The electric field is too weak to provide enough energy for electrons produced by the background radiation to ionize other neutral particles. In the C E sub region, known as the Townsend discharge, as the voltage is increased the current profile is characterized by an exponential increase. The electrons present in the gas are now able to gain enough energy from the electric field to ionize the neutral particles. As the voltage is further increase d the secondary electrons produce during ionization are more likely to become part of the ionization process along wit h the electrons produced by background radiation leading to an avalanche of electron and ion production. T he D E sub region the c orona discharge is part of the Townsend discharge regime and occurs in areas of high electric field near sharp points, edges, or wires preceding to electric breakdown ( at point E). If the current is high enough, corona discharges can become visible to the human eye emitting light in the UV spectral region and the visible blue spectral region, but for low currents they are dark discharges. At p oint E electrical breakdown occurs in the Townsend regime due to the production of secondary electrons generated by ion or photon impacts on the cathode electrode At point E, when the voltage reaches the sparking potential ( V B ), the current increases by a factor of 10 4 to 10 8 Only the internal resistance of the power supply limits the current. When the discharge tube cannot draw enough current to break down the gas, due to high internal resistance, the discharge remains in the corona regime. Otherwise, the discharge transition s into the normal glow regime due to electron PAGE 56 56 avalanche processes The pressure distance (pd) product and the secondary electron coefficient of the electrode material ( sec ) determine the breakdown voltage ( V B ) f or a [41] Glow Discharge In the glow discharge region, the electron energy and electron number density are high enough for excitation collisions to generate visible light causing the plasma to be luminous In this regime, the current is usually between 10 5 A and 1 A for a typical laboratory cold cathode tube The transition between point E and F is discontinuous. In the F G sub region, called the normal glow regime the voltage is almost independent of the current and the electrode current density is independent of the total current because at low currents the plasma is in contact with only a small area of the cathode. As the current increases, the cathode area in contact with the plasma increases, until point G is reached and the entire cathode area is in contact with the plasma. The transition from normal to abnormal glow discharge happens between current values of 10 3 to 10 2 A. In sub regio n G H, called the abnormal glow regime t he plasma covers all the cathode area, a higher current is only achievable by having a higher current density. The voltage increases significantly with the increasing current. The abnormal glow regime discharge is brighter than the normal glo w regime, and the structures of the discharge near the cathode can blend into one another, forming a uniform glow. A form of hysteresis, a time and past inputs, is observed in the voltage current characteri stics when moving to the PAGE 57 57 left of point G. The discharge reaches a lower current value and current density than at point F, making a discontinuous transition back to the Arc Discharge In the H K region, the arc discharge region after point H the cathode becomes hot enough to emit electrons in a process called thermionic emissions. The discharge undergoes glow to arc transition, in the H I regime when the power supply has low internal resistance at currents of ~ 1 A In t he I J regime, known as the low intensity arc, the discharge voltage decreases as the current increases the arc column is governed by thermal conduction losses [45] and thermionic emission of electrons from the electrodes support t he arc. After point J, which is at approximately 20 to 50 A, the voltage increases slowly as the current increases radiation becomes the dominant loss mechanism yielding the positive gradient characteristic [45] Qualitative C haracteristics of DC Glow Discharge DC glow discharges produce weakly ionized ( low degree of ionization 10 8 10 4 ) non others) are in thermodynamic equilibrium and are characterized by a single temper ature. In a non thermal plasma, the electron temperature is much higher (~ 10,000 K to ~ 100,000 K) than the temperature of the heavy particles which is roughly the same as the room temperature [11] The degree of ionization re fers to the proportion of neutral particles that are ionized into charged particles and it is defined as, PAGE 58 58 where n i and n n represent ion and neutral number densities, respectively [8] The degree of ionization for DC glow discharges are in the range of 10 8 10 4 Due to the low degree of ionization, recombination is of minor importance and most of the charged particles are lost (neutralized) by transport to the solid surfaces. The usual pressure range of operation is between 10 mTorr and 10 Torr and typically, only a few hundred s of volts between cathode and anode is required to maintain the discharge [41] In a laboratory discharge, the cathode is an electrical conductor, usually a metal, with a secondary emission coefficient; it em its electrons when bombarded by incoming positive ions. The behavior of a DC glow discharge can be quite complicated, with many light and dark regions. Figure 3 2 shows the qualitative characteristics of DC glow discharge in a glass tube containing low pressure gas [46] Figure 3 2 Qualitative characteristics of DC glow discharge and electric potential distribution at low pressures (0.1 10 Torr) and room temperature [46] ( 3 1 ) PAGE 59 59 All of the regions are gas, pressure, and voltage dependent in their size and intensity, with some of the smaller features being essentially absent over various parameter ranges [41] The Aston dark s pace is a thin region next to the cathode containing a layer of negative charge and a strong electric fie ld. Electrons are accelerated through this space away from the cathode. The plasma appears dark since the electron density and/or energy is too low to excite the gas In this region, the combined electrons number density of stray initial electrons and secondary electrons from the cathode outnumber the ion number density causing the negative charge In the cathode g low, the electrons are energetic enough to excite the n eutral atoms resulting in the luminous cathode glow. This region has relatively high ion density and sometimes masks the Aston dark space by clinging to the cathode. Its axial length depends on the carrier gas, the pressure, and temperature. Cathode dark s pace (Crooks dark s pace) is a positive space charge layer that has a moderately strong electric field and a relatively high ion density [47] Its axial length depends on the pressure and the applied voltage. In this region, el ectrons are accelerated by the strong electric field from the cathode to the anode, while the positive ions are accelerated in the opposite direction. The se accelerated ions cause the pulverization of the cathode material and the emission of secondary elec trons. The secondary electrons accelerate and cause the creation of new ions through collision with neutral atoms. The majority of the potential drop happens in this region [41] From the end of the PAGE 60 60 cathode dark space to the anode, the plasma experiences only low values of electric fields [47] The negative glow region has the brightest intensity in the whole discharge. Electrons that have been accelerated in the cathode region to high speeds produce ionization, while slower electrons that have had experienced inelastic collisions produce excitations [41] The slower electrons are responsible for the glow although other processes also play a role [47] In this region the electrons predominantly carry the current due to their higher mobility At the end of this region, the electric field must decrease rapidly and the electrons slow down by elastic and inelastic collisions The electric field can reverse in this region to keep the electron current balance [41] The field reversal excludes ions from a region of the column; this region of low ion density prevents the negative glow from joining directly onto the positive column [41] In the Faraday dark space, a transition region between the negative glow and the positive colu mn, the electrons do not have e nough energy to produce excitations, therefore, no glow is observed [41] In this region, the electrons begin to feel the effects of the positive anode potential and are accelerated toward it, alt hough there is only a relatively low field strength [47] The positive column is typically the physically largest component of a normal discharge. If the discharge tube is increased at constant pressure, the cathode structures remain the same size, but the positive column lengthens. In this region, the plasma is quasi neutral and the electric field is weak. Usually, the positive column is a long, uniform glow mode discharge, except when standing or moving striations (traveling o r stationary perturbations in the electron number density) or other PAGE 61 61 instabilities are triggered by disturbance [43] Striations occur due to electron absorption during ionization and the formation of slow secondary electrons that require acceleration before they can cause ionization. The Anode glow is usually brighter than the positive column, and it is not always pre sent. It is the boundar y of the anode sheath The Anode dark space ( anode sheath ), is a single layer of negative space net charge density that separates the anode glow from the anode itself This region has stronger electric fie ld than the positive column Figure 3 2 shows the typical potential and light distribution when the distance between electrodes is large in comparison to the size of the electrodes. If the electrodes are brought closer together, the positive column shrinks and eventually is completely eliminate d. In this case, the anode is in contact with the Faraday dark space and the anode is at equipotential with the bulk plasma. If the distance between electrodes is further decreased, the Faraday dark space also disappears and the anode is in contact with th e negative glow, which is called obstructed glow [48] In this case, the negative glow region has the most positive potential and the anode zone has a negative potential fall that repels electrons and attracts positive ions [48] The o bstructed glow is the most common mode of operation for plasma processing used in industry [48] In the previous discussion, the typical characteristics in the normal glow mode have been considered, which usually occurs between the range of current densities of 10 5 to 10 3 A/cm 2 [41] The abnormal glow discharge appears very similar to the normal glow discharge with the exception o f being more intensely luminous, a s the normal glow transitions to abnormal glow, the cathode voltage drop increases rapidly, and the PAGE 62 62 cathode dark space shrinks. Sometimes, the structures near t he cathode blend into one another, providing a uniform glow. For further information about the electric breakdown of gases and the qualitative characteristics of DC glow discharge s, consult Conrads & Schmidt [49] Fridman [50] Howatson [51] Raizer [43] and Madou [48] Previously, the was describe d as the equation determining the breakdown voltage of a given gas. The introduced by Friedrich Paschen to predict the breakdown voltage (V B ) the voltage necessary to start a discharge, between two electrodes in a gas as a function of the pressure and the gap length (pd) The formula that describes the breakdown voltage is defined as [41] For large values of (pd) V B increases essentially linearly with (pd) There is a limiting value of (pd) below which breakdown cannot occur according to The breakdown voltage is a minimum ( V min ) at some intermediate value of ( pd ) ( 3 2 ) ( 3 3 ) ( 3 4 ) ( 3 5 ) PAGE 63 63 The breakdown voltage curve is called the Paschen curve, and it is a function of the gas and weakly a function of the electrode material [41] The following table shows the values of A and B coefficients for several gases obtained from page 546 in Lieberman and Lichtenberg [41] For a rgon, the following table gives the V min values and the values of ( pd ) at which V min occur s for several secondary electron emission coefficients ( sec ) as well as the value of (pd) Limit below which the Paschen curve can no longer predict the breakdown voltage Table 3 1 Breakdown voltage constants A and B [41] Gas A (cm 1 Torr 1 ) B (Vcm 1 Torr 1 ) Range of E/ p (Vcm 1 Torr 1 ) He 2.8 77 30 250 Ne 4.4 111 100 400 Ar 11.5 176 100 600 Kr 15.6 220 100 1000 Xe 24 330 200 800 H 2 4.8 136 15 600 N 2 11.8 325 100 600 O 2 6.5 190 50 130 CH 4 17 300 150 1000 CF 4 11 213 25 200 Table 3 2 Limit of (pd) and minimum breakdown voltage va lues for argon. ( sec ) (pd) Limit (pd) Vmin V min (Torr cm) (Torr cm) (V) 0.01 0.401 1.091 192 .0 0.0 5 0.265 0.720 126.7 0.07 0.237 0.645 113.4 PAGE 64 64 Figure 3 3 Paschen curves for argon at different values of sec As observed in Figure 3 3 the breakdown voltage depends strongly on the strength of the applied electric field, increasing with increasing electrode s pacing for (pd) values higher than (pd) Vmin For (pd) values lower than (pd) Vmin the breakdown voltage increases with decreasing electrode gap due to a decrease in the probability of collisions between electrons and atoms. For the cases shown in this stud y, the distance between electrodes is chosen in order to satisfy (pd) Vmin at the given inlet pressures and the initial voltage of each numerical simulation is chosen higher than V min Sputtering At energies above a threshold of thr = 20 50 e V [41] heavy particles can sputter atoms from a surface. The sputtering yield ( sput ) represents the number of atoms sputtered per incident ion. Sputtering can cause erosion of electrodes, which eventually can damage a thruster. In a n argon glow discharge, operating at 1.0 Torr and PAGE 65 65 ~950 V voltage difference, the average ion bombardment energy is ~60 eV, which produces a primary sputter yield of ( sput = 0.1) for copper electrodes [52] The sputtering rate, R sput can be calculated from the following equation [41] w here is the ion current density in (A/cm 2 ) and n Target is the number density of the target material in (cm 3 ) For example, the number density of Cu is equal to 8.49 10 22 cm 3 and the primary sputter yield is equal to 0.1 For an ion current density of 5 .0 10 3 ( A/cm 2 ), the peak value of the ion current density for Case 750V in this study the R sput is approximated as 3 68 10 8 cm/s. For copper electrodes, the cathode region experiencing peak ion current density will be eroded 1 mm every 755 hours of continuous thruster operation. Thus the cathode must be replaced when the erosion track becomes comparable to the cathode thickness or t he thruster will be impaired. The erosion rate can be decreased by using more suited materials with lower sputtering yield such as tungsten. ( 3 6 ) PAGE 66 66 CHAPTER 4 GOVERNING EQUATIONS In order to simulate the RGEJ, two different numerical models are used: an ionized gas model, which uses a hydrodynamic plasma model called the Local Mean Energy Approximation (LMEA) to simulate the plasma (charged and excited particles) and a rarefied gas model to simulate the gas itself (the aggregate of all particles as a single mixture) The ionized gas model is used to calculate the volumetric plasma induced electrostatic force components ( F x F y ) and the volumetric electro thermal heating source ( ) needed by the rarefied gas model. The rarefied gas model is used to calculate the gas density ( ), co mponents of velocity ( u, v ), and temperature ( T ) needed by the ionized gas model. In this chapter an overview of the governing equa tions, boundary conditions and equations used to obtain F x F y and Hydrodynamic Plasma Model Boltzmann Equation Since a plasma is an assembly of a great number of moving and interacting particles, plasma problems can be solved by methods of physical kinetics [53] The distribution function is the main statistical characteristic of an ass embly and to find it the whole configuration space and the entire range of particle velocities must be divided into small intervals such that the particle density variation within each interval could be neglected. The distribution function ( F ) can be repre sented as the product of the density of the particles in the configuration space ( n ) and the velocity distribution function ( f ) as described by Golant et al. [53] PAGE 67 67 The equation describing the time evolution of the velocity dis tribution function ( f ) for a particle specie s ( j ) is known as the kinetic (or Boltzmann ) equation [53] [54] ( 4 1 ) The vector x describes the spatial coordinates, v describes the velocity coordinates, t is the time, E and B are the electric and magnetic field s, and n j is the number density [54] Z j represents the ionizity (or valence) of the charged particles and m j is the species mass. The sign of Z j denotes the polarity of the charge. The R.H.S. term represents the rate of change in distribution function due to the particle collisions [53] Since the different collisions occur independently of each other, the rate of change in distribution function because of various types of collisio n can be summed. F or type particles ( 4 2 ) w here the summation is extended to all species of type particles and all types of p collisions. Since f j ( x v j t ) varies in up to six dimensions (3 D and 3 V) and time, solving this equation along with the additional coupled governing equations of electric and magnetic fields is computationally prohibited [54] To reduce the high dimensionality, the velocity distribution function is projected and integrated in velocity space. The necessary normalization condition is [53] ( 4 3 ) PAGE 68 68 The distribution function makes it possible to average out any value over the particle distribution. T he average value of a quantity in velocity space can be defined as follows (using velocity as an example) [53] [54] ( 4 4 ) These two equations previously shown are used to integrate the Boltzmann equation and allow more realistic fluid and plasma problems to be solved. Distribution Function Moments Equation In order to obtain the general form of the moments equation t he Boltzmann equation is multiplied by some combination of the velocity projections denoted by g corresponding to a certain time instant, and then each term in the equation is integrated over the entire velocity space [53] T he distribution function moment correspondi ng to this combination is ( 4 5 ) After ( 4 1 ) is multiplied by g and integrated over the velocity space the first term in the equation is ( 4 6 ) PAGE 69 69 The order of integration and differentiation can be reversed since time t and velocity comp onents v are independent variables. T he coordinates in configuration space x is also independent of v therefore the second term is transform ed as ( 4 7 ) where k is the number of components of v The third term, which is related to the Lorenz force, can be separated in to two parts: the one proportional to the electric field and the other one proportional to the magnetic field. Using integration by parts and taking into account that the distribution function must vanish as each part becomes ( 4 8 ) ( 4 9 ) PAGE 70 70 In ( 4 9 ) the term is equal to zero since depends only on the velocity component perpendicular to v k The final term is obtained by integrating the collision term, ( 4 10 ) The general form of the moments equations is [53] ( 4 11 ) Using ( 4 11 ) the zero moment equation, defined as the equation for the concentration (density), can be derived by setting g equal to 1, the first moment equation is derived by setting g equal to v k in the se cond moment equation g equal s v k v l and so on. Continuity Equation As previously explained, the zeroth moment (or continuity) equation is obtained by setting g = 1, which makes the average of the quantities contained in ( 4 11 ) equal to [53] ( 4 12 ) The term S represents the chemical reactions (gains and losses of particles due to collisions) [54] To evaluate the collision integral, some assumptions about the velocity PAGE 71 71 distribution function are necessary such as its deviation is small, which allows the Boltzmann equation to be solve d by the method of successive approximations. In this case, the dist ribution function is represented as a series in the powers of the parameters (f = f (0) + f (1) + f (2) affecting the particles and the concentration and temperature gradients) [53] The first term of the series (f (0) ) is the equilibrium (Maxwellian) di stribution, the second term (f ( 1 ) ) includes a linear combination of the parameters, the second term (f ( 2 ) ) a quadratic combination, and so on. If a homogeneous electric field E is the source of disequilibrium, the field determines the only isolated direction, therefore it is assumed that the 0z axis is parallel to it. This means that f may depend only on the velocity v and on the angle due to anisotropy, the angle between the directions v and E T he distribution function is expanded using a series, typically the orthogonal Legendre polynomials ( 4 13 ) T he function f n depend s only on the velocity. The Bol t zmann solver BOLSIG + used to calculate the reaction rates in this study, uses a two term expansion ( 4 14 ) The final zeroth moment (or continuity) equation for a species (j) is PAGE 72 72 ( 4 15 ) The species velocity u j introduced as an unknown, is solved using the first moment equation. Momentum Equation The first moment (or momentum) equation is obtained by setting g = v k which makes the average of the quantities contained in ( 4 11 ) equal to [53] ( 4 16 ) In ( 4 16 ) is a unit vector in the direction of the k axis, while the set of P kl values is the momentum flux density tensor [53] In t he momentum flux density tensor p is the normal scalar pressure and is the set of values of the viscous stress tensor due to the deviation of the distribution from spherical symmetry, where w is the random velocity vector. After substituting ( 4 16 ) in ( 4 11 ) t he first moment equation is defined as ( 4 17 ) PAGE 73 73 The first term in ( 4 17 ) is expanded and transform with the help of the zeroth moment equation ( 4 18 ) wh ile the last term in ( 4 17 ) can be decomposed as ( 4 19 ) Substituting ( 4 18 ) and ( 4 19 ) into the first moment equation and finding identical equations for the other two components of the vector u the vector equation for the directed velocity for a specie s (j) is [53] ( 4 20 ) In th is equation, the viscous effects due to can be neglected if the characteristic length scale of the bulk velocity change is much larger than the mean free path, and the random thermal motion is much faster than the bulk fluid velocity [53] These conditions are satis fied since the velocity distribution function is near Maxwelli an (small anisotropy ) Furthermore, in the absence of magnetic field, the pressure along all directions can be related to the temperature using the state equation ( 4 21 ) PAGE 74 74 where k B is the Boltzmann constant (1.602 10 19 J/eV) and T j is the particle temperature (eV) which requires an additional moment equation to be calculated [54] In a weakly ionized plasma, the charged neutral particle collision frequency is much larger than the collision frequency between charged particles [54] If collisions between charged particles are neglected and the distribution function over the random velocities can be considered isotropic then [53] ( 4 22 ) where ( ) is the charged neutral particle collision frequency for a specie s (j ) and U n is the neutral velocity vector (approximated as the bulk gas velocity vector ) Energy Equation The second moment (or energy) equation is obtained by setting g = mv 2 /2 which makes the average of the quantities contained in ( 4 11 ) equal to [53] ( 4 23 ) After substituting ( 4 23 ) in ( 4 11 ) the second moment equation is obtained as PAGE 75 75 ( 4 24 ) where is the sum of the energies due to random and directed motion and Q is the energy flux vector [53] ( 4 25 ) which depends on the heat flux density vector ( q ), the pressure tensor ( ), and other previously defined variables. The pressure tensor is simplified (the viscous stress tensor is neglected) by assuming that the particle velocity distribution function is near isotropic (small anisotropy) [53] T he heavy ion and metastable atoms species are as sumed to be in thermal equilibrium with the neutral gas. The feasibility of this assumption depends on the number of collisions between ions (or metastable atoms) and the neutral particles in a weakly ionized plasma [54] In th e case s shown in this study, the ions rapidly equilibrate with the gas since the ion neutral mean free path is orders of magnitude smaller than the device length scale me aning that many collisions happen between an ion and multiple neutral particles before the ion reaches one of the wall s [55] For the electrons, a simplified energy equation must be solve d This approach is called the Local M ean Energy Approximation (LMEA) [54] For the electrons, a common approximation is made T he kinetic energy due to random motion is considered much PAGE 76 76 larger than the kinetic energy due to directed motion, which makes the total ki netic energy approximately equal to the mean energy ( ) ( 4 26 ) After applying all the assumption s previously explained, the energy flux vector becomes ( 4 27 ) The second moment energy equation for electrons is ( 4 28 ) In the present model, the fourth moment, known as the electron heat flux ( ), is approximated using a closure model instead of solving a separate conservation equation, similarly to Houba [54] ( 4 29 ) Drift Diffusion Approximation The drift diffusion approximation is used to drastically reduce the computational cost of a plasma numerical simulation by eliminat ing the need to solve the full momentum equation separately since the velocity is calculated as a direct function of PAGE 77 77 the other variables. Several assumptions are necessary for this approximation. As previously mentioned, the viscous terms ( ) are assumed negligible in the momentum equation and the pressure along all directions can be related to the temperature using the state equation ( ) The second term ( ) in the momentum equation (quadratic with respect to u ) is neglected since it is much less than the term proportional to the pressure gradient ( ) when the anisotropy of the velo city distribution is small (L is the characteristic length over which the plasma parameters change substantially) [53] For a weakly ionized gas, t his assumption is true for electrons, which are not thermalized with the heavy p articles For ions, which are thermalized with the gas, the term ( ) is much smaller than ( ) T he first term of the momentum equation ( 4 20 ) is considered negligible when the characteristic time of variation in plasma parameter ( ) greatly exceeds the time between co llisions [53] For electrons, due to their very large electron neutral collision frequency ( ) the first term is typically negligible [43] If the frequency of the applied electric field is not large in comparison to the ion neutral collision frequency t he first term of the momentum equation is negligible even for heavier species [56] For the cases simulated in this study, the applied electric f ield varies slowly enough in the RF cases to neglect this term (the ion neutral collision frequency is about twenty times higher than the applied RF frequency of 13.56 MHz) In PAGE 78 78 the absence of magnetic field, taking the previous assumptions, the following m omentum equation is obtained ( 4 30 ) After ( 4 30 ) is di vide d by the collision frequency and the reduced mass ( m jn ) the species flux is obtained as ( 4 31 ) To obtain the final drift diffusion form, the mobility, diffusion and thermal diffusion coefficients are defined as ( 4 32 ) There exists a defined relationship between these three coefficients [53] For particle velocity distributions close to Maxwellia n the mobility and diffusion coefficients are related as PAGE 79 79 ( 4 33 ) which is known as the Einstein relation and For electrons and and for ions and The thermal diffusion coefficient is equal to [53] ( 4 34 ) The coefficient ( ) has an accuracy of the order of unity, to find a more accurate value the motion and heat flux equations must be solved simultaneously [53] Using this estimate, the thermal diffusion can be written as ( 4 35 ) The drift diffusion approximation of the species flux is ( 4 36 ) For electrons, the third term in ( 4 36 ) called the thermodiffusion flux, is usually smaller than the diffusion flux and it is typically neglected [43] For all the cases shown in this study, the previous statement is also true for ions, therefore PAGE 80 80 ( 4 37 ) Using the drift diffusion approximation of the species flux, the continuity equation for each species can be solved without the corresponding momentum equation s ( 4 38 ) Likewise, the drift diffusion approximation can be used in the electron energy equation. The energy equation can be cast in a form similar to the continuity equation with the help of the mean energy density ( ) as described by Hagelaar and Pitchford [57] ( 4 39 ) where the electron energy mobility and diffusion coefficients are related to the electron transport properties as [57] ( 4 40 ) PAGE 81 81 that induce forces on the charged particles in a plasma. In a linear medium, where the electric polarization and magnetization are proportional to the electric and magnetic fie [58] ( 4 41 ) ( 4 42 ) ( 4 43 ) ( 4 44 ) In these equations, ( ) and ( ) are the permittivity and permeability of the material. For an argon plasma, ( ) and ( ) are the permittivity and permeability of free space The net charge ( charge ) and net current ( J net ) are functions o f the number densities and fluxes of the charged particles, respectively, and are obtained from the mass conservation equation, defined as [54] ( 4 45 ) PAGE 82 82 ( 4 46 ) For all the cases in this study, no external magnetic field B is applied. The induced magnetic field generated by the time derivative of equation ( 4 43 ) is characterized by the magnetic Reynolds number defined as [54] ( 4 47 ) where U and L are typical velocity and length scales of the flow The variable ( ) is the electrical conductivity of the material in question, which in the case of weakly ionized plasmas is equal to ( 4 48 ) If the magnetic Reynolds number is small ( ) the induced magnetic field is negligible due to the smoothening effect of the magnetic diffusivity ( ). For all the cases shown in this dissertation, the maximum magnetic Reynolds number, based on the electron velocity ( U=u e ) is of an order of 10 2 Therefore, the induced magnetic field is neglected. In the absence of magnetic fields, equation ( 4 42 ) is satisfied automatically. ( 4 43 ) is s atisfied by defining an electric potential ( ) such that ( ) since PAGE 83 83 ( 4 49 ) ( 4 44 ) is satisfied by solving the system of continuity equations ( 4 38 ) E ach species continuity equation is multiplied by the species charge ( ) before they are added together ( 4 50 ) w here as long as the net charge is conserved through the constraint [54] ( 4 41 ) is the only equation that needs to be solved in addition to the conservation laws, completing the coupling between the species number density and the electric field ( 4 51 ) The governing equations defined by ( 4 38 ) ( 4 39 ) and ( 4 51 ) are written in the non dimensional form using reference number density of 10 16 (m 3 ), an electric potential of 100 (V), a length scale of 0.01 (m), an electron temperature of 10 (eV) and a time scale of 10 10 (s). PAGE 84 84 Transport Properties The solution of the system of equations defined by ( 4 38 ) ( 4 39 ) and ( 4 51 ) requires models for the transport coefficients ( j and D j ) and the reaction rates ( S j ). In the LMEA model the electron transport coefficients and reaction rates are calculated as functions of the mean electron energy ( ), which is computed from the electron energy equation given that ( ) The coefficients are tabu lated using an electron Boltzmann equation solver called BOLSIG + [57] The ion mobility coefficient is obta ined from Rafatov et al [59] as ( 4 52 ) where the electric field ( E ) is in (V/m) and the neutral number density ( n n ) is in (m 3 ), giving a mobility in (m 2 /Vs). The diffusion coefficient is obtained from Einstein relation ( ) The mobility is zero for all the metastable species in this model and the metastable diffusion is the same as in Deconinck [11] ( 4 53 ) PAGE 85 85 w here ( ) is the collision frequency between metastable species ( j ) and the background neutral species ( n n ), ( ) is an approximate hard sphere momentum transfer collision cross section based on the Lennard Jones in teraction potentials [60] and ( ) is the average relative molecular speed. The electron energy transport coefficients ( and ) may be calculated using BOLSIG + software [54] Otherwise, the approximations ( 4 40 ) may be used to decrease the nonlinearity of the numerical simulation. For an argon plasma with ions, electrons, and three types of metastable species, t he neutral number density ( n n ) is ca lculated by solving the state equation ( 4 54 ) where t he subscript s ( i ) and ( e ) represent ions and electron, while the subscripts ( S P and D ) represent three types of metastable atoms named after their energy levels [61] The reaction model includes elastic collisions and collisions that cause ionization, excitation, transition, and spontaneous emission The elementary react ions considered for low pressure, argon chemistry are given in Table 4 1 [57] The cross sections for these reactions were obtained from Biagi [62] and the kinetic model section of Petrov and Ferreira [61] marked with superscripts B and P, respectively PAGE 86 86 The source terms ( S j ) in the continuity equations are determined by the reactions occurring in the discharge described in Table 4 1 where the source terms for ions and electrons are identical d ue to the particle conservation, Table 4 1 Relevant Reactions. All reaction rates calculated using B OLSIG + unless otherwise stated [57] Index Reaction Type (eV) Coefficient 1 E lastic collision 0 Boltz. B 2 D irect ionization 15.7 Boltz. B 3 E xcitation 11.55 Boltz. B 4 E xcitation 13.0 Boltz. B 5 E xcitation 14.0 Boltz. B 6 S tepwise Ionization 4.07 Boltz. P 7 S tepwise Ionization 2.52 Boltz. P 8 S tepwise Ionization 1.66 Boltz. P 9 T rans. b/w Excited States 1.51 Boltz. P 10 T rans. b/w Excited States 0.90 Boltz. P 11 P enning Ionization 5. 0 10 16 m 3 /s P 12 P enning Ionization 5.0 10 16 m 3 /s P 13 P enning Ionization 5.0 10 16 m 3 /s P 14 T rans. b/w Excited States 5.0 10 44 m 6 /s P 15 T rans. b/w Excited States 5.0 10 17 m 3 /s P 16 R adiation 7.0 0 10 8 1/s P 17 R adiation 5.96 10 8 1/s P 18 R adiation 3.76 10 8 1/s P 19 R adiation 1.46 10 8 1/s P B C ross sections for these reactions were obtained from Biagi [62] P C ross sections for these reactions were obtained from Petrov and Ferreira [61] PAGE 87 87 ( 4 55 ) ( 4 56 ) ( 4 57 ) ( 4 58 ) where the subscript ( rev ) in some of the reaction rates ( R ) is for the reverse reaction In order to complete the model, the source term for the electron energy equation is defined as [11] ( 4 59 ) where N is the number density of the heavy particle in the given reaction, is the electron atomic elastic collision and is the energy loss (or gain) due to inelastic collision s during the reaction ( r ) These two terms in the equation represent the changes in electron energy du e to elastic collisions and inelastic processes, respectively. Boundary C onditions In order to solve the drift diffusion model a set of boundary conditions is necessary. The imposed boundary conditions are different depending on the behavior of the syst em near an anode electrode, cathode electrode, dielectric wall, or an open PAGE 88 88 surface exposed to the gas. A more extensive discussion of the typical boundary conditions for fluid models of gas discharge is presented by Hagelaar et al. [63] Open surface b oundary If the boundary is an open surface exposed to the gas, then the flux components due to the mobility and the diffusion normal to the boundary are equal to zero in ( 4 37 ) and ( 4 39 ) Similarly, the electric field (or potential gradient) normal to the boundary is ass umed to be zero because the given boundary is sufficiently far away from the area of high number density of charged particles (near the cathode) For the 2 D cases, the inlet boundary and the bottom of the computational domain have ( 4 60 ) where ( ) is the outward normal unit ve ctor to the boundary in question. This condition is imposed at the inlet due to numeral stability reasons. T he exit plane of the channel has ( 4 61 ) When the solution reach es steady state, the number densities are so small that th e s e boundary conditions become the same as ( 4 60 ) The electric field normal to the boundary is ( 4 62 ) PAGE 89 89 Anode b oundary The flux of particles towards a surface can be derived using kinetic theory as ( 4 63 ) where the coefficient ( ) is set equal to zero if the electric field normal to the wall multiplied by the sign of the charge of the species is in the direction pointing away from the wall, and it is set equal to one otherwise [54 ] The variable is the average thermal velocity of species ( j ). The diffusion of the charged species in the sheath region (near the wall) is smaller than the drift flux and it is usually neglected [64 ] For the electrons, ( ) for all the cases in this study due to the short distance betw een electrodes, which makes all the cases operate in an obstructed glow mode (previously explained in the section labeled as Qualitative C haracteristics of DC Glow Discharge ). The energy flux towards a surface has the form ( 4 64 ) The expression for the energy flux is not equal to the product of the wall flux ( ) with the mean energy ( ) because particles with higher energy hit the wall more frequently and contribute more power to the wall [44] PAGE 90 90 The electric potential ( ) at the anode is fixed to an applied value (Dirichlet boundary condition or first type ), which may vary in time for the RF cases Cathode b oundary The cathode boundary conditions are identical to the anode boundary conditions exc ept for an additional term in the electron particle flux and a corresponding additional term in the electron energy flux. These terms are due to the secondary electron emission, characterized by the secondary emission coefficient ( sec ) which depends on the cathode material and gas [54] The cathode boundary condition for electrons is ( 4 65 ) Dielectric w all The particle and electron energy fluxes for a dielectric wall are identical to the ones at the anode. At the interface between two materials (in the RGEJ case dielectric above and plasma below the interface) with different permittivity constants, the electric field is defined as [65] ( 4 66 ) PAGE 91 91 where ( ) is the c harge on the dielectric surface calculated by assuming that the charged particles striking the wall recombine instantly and no diffusion of charge along the wall occurs [54] The surface charge is calculated as [11] ( 4 67 ) For all cases shown, the electrodes are exposed instead of embedded inside the dielectric material. The boundary condition used in the RGEJ cases is ( 4 68 ) where it is assumed that the contribution of the dielectric material to the dielectric boundary equation is negligible due to several thrusters being stack together (similar to the design of FMMR) and separated by a relatively thin dielectric wall. Plasma Electrostatic Force and E lectrothermal Heating Source The volumetric plasma induced electrostatic force and the volumetric electro thermal heating source are calculated as [11] [66] ( 4 69 ) ( 4 70 ) PAGE 92 92 In equation ( 4 70 ) the term on the RHS is the ion Joule heating, since only the heavy particles thermalize with the bulk gas. The ion Joule heating term usually has a thermalization factor that accounts for the fraction of energy that is locally e quilibrated with the gas. In a xenon discharge at 100 (Torr), it was assumed by Boeuf et al [67] that only 25% of the energy was deposited in the neutral gas and the maximum gas temperature (~460 K) of their simulation matched relatively well with the gas temperature (~500 K) obtained from experiments. In this case, the ion Joule heating rapidly equil ibrates with the gas since the ion neutral mean free path is orders of magnitude smaller than the device length scale, a conservative approach is taken to obtain the upper bound for the gas heating, and the thermalization factor is assumed to be equal to one [66] Rarefied Gas Equations In the rarefied gas simulation, density based compressible flow equations were used with the assumption of ideal gas using argon as the working fluid. The continuity, momentum, and energy equations are given by Raju [68] with the additions of the terms ( F x F y ) in the right hand side of the momentum (in the x and y ) and energy equations, respectively. The c ontinuity equation compressible m omentum equations in 2 D (in the x and y ) with ideal gas assumption, and e nergy equation are ( 4 71 ) PAGE 93 93 ( 4 72 ) ( 4 73 ) ( 4 74 ) In these equations, ( ) is the density, ( u ) and ( v ) are the x and y components of velocity, ( T ) is the gas temperature, R is the gas constant (R=208. 1 J/ kg K) for a rgon, ( c v ) is the specific heat at constant volume , where the specific heat ratio ( ) is equal to 1.667 T he thermal conductivity ( ) and the viscosity ( ) are functions of the gas temperature [69] ( 4 75 ) ( 4 76 ) w here ( ) is in (Pa s), ( ) is in (W/m K), and (80 K T K ) PAGE 94 94 The governing equations are written in the non dimensional form using a velocity of 100 (m/s), a length scale of 0.01 (m), a pressure of 100 (Pa), a temperature of 300 K, and density from the ideal gas law using the s e parameters Boundary c onditions for r arefied g as : The Knudsen number ( Kn ) is defi ned as the ratio of the molecular mean free path length ( ) to a representative physical length scale, which in this case is the height ( H ) of the inlet of the micro channel and it is much smaller than the micro channel length The Knudsen number and mean free path are defined as [70] [71] ( 4 77 ) ( 4 78 ) where d and n are the atomic diameter and gas number density. Kn is used to determine which numerical model ing approach is more appropriate: statistical mechanics or continuum mechanics. The figure below shows the equations use to model different flow regimes. Figure 4 1 Gas flow regimes depending on Kn [68] PAGE 95 95 As Kn 0, the flow is assumed sufficiently continuous, while for Kn > 10, the flow is assumed free molecule. For 10 3 < Kn < 10 the flow is neither sufficiently continuum nor completely molecular [68] For this range, the flow is further divided into two subcategories: slip flow regime 10 3 < Kn < 10 1 and transitional regime for 10 1 < Kn < 10 as explained by Raju [68] In the RGEJ cases, Kn is in the slip regime by design. The boundary conditions for the rarefied gas simulation s of the RGEJ are fixed stagnation density and temperature at the inlet ( =P 0 /RT 0 T 0 =300 K, where P 0 =133.3 Pa) Isentropic flow assumption is used to calculate the static density and temperature at the inlet plane. At the walls there is no penetration, the normal velocity is equal to zero. At the outlet, static pressure is assumed to be ( P Out =0.05 Torr ) if Mach number is subsonic and ( =P Out /RT Out ), else is extrapolated from internal nodes. Boundary conditions for a rarefied gas are used for tangential velocity and temperature at the wall face, as described by Maxwell [72] and Smoluchowski [73] si milar to Raju [68] For example, for the top wall face inside the channel, the tangential velocity boundary condition is ( 4 79 ) The corresponding temperature boundary condition for an isothermal wall is ( 4 80 ) PAGE 96 96 For these equations, ( r ) is the distance from the wall normal to it, for the top surface ( r = H y ). If a case is adiabatic, or is used as a temperature boundary condition at the wall. In these equations, the new variables introduced are the Prandtl number ( Pr ), the specific heat ratio ( ), the tangential momentum accommodation coefficient (0 v 1), and the thermal accommodation coefficient (0 T Specular reflection happens when the gas molecules are reflected from the wall at an angle equal to the incident angle exerting no shear stress on the wall. Di ffuse reflection happens when the channel surface is rough and the gas molecules are reflected at random angles. The accommodation coefficients indicate the fraction of the molecules reflected diffusively from the walls [68] The tangential momentum accommodation coefficient ( [74] ) and the thermal accommodation coefficient ( [75] ) found in the formulas of slip flow boundary conditions are selected based on average values for argon interacting with different materials. Using the definition of mean free path ( 4 78 ) and th e expression of the viscosity of dilute gases [71] ( 4 81 ) obtained using the Chapman Enskog method where m is the molecular mass and a is the atomic radius, the boundary condition formulas can be modified. PAGE 97 97 The boundary conditions can be solved for the fluxes as [68] ( 4 82 ) ( 4 83 ) which allows the boundary conditions to be implemented into the simulation through the weak statement of the Galerkin m ethod. The rest of the boundary conditions needed are zero flux normal to the edges of the domain [ or ]. Heat Transfer Governing Equations For micro propulsion devices, heat loss via radiation and conduction are major concerns. As stated in chapter 1 FMMR was shown to have significant heat loss, a pro blem that ultimately prevented its application. In order to estimate the heat loss through radiation and conduction in the RGEJ the temperature distributions inside the channel walls and the external walls are used as boundary conditions or inputs to perform a thermal analysis These temperature distributions are obtain ed from the numerical simulation of the loosely couple gas modules of the numerical code The thermal analysis is not couple with the gas modules. Thermal Radiation Heat Transfer In the thermal radiation analysis a simplified model is used with the assumption of grey diffuse surface s with a non participating medium, allowing the radiation properties to be independent of wavelength and direction. To model the problem the PAGE 98 98 equation for radiative exchange between gray, diffuse surfaces was used and written in the following way ( 4 84 ) In this equation, are the emissivities of the given surface, are the unknown heat flux es are the known emissive powers, are the external irradiations. Additionally, is the Kronecker delta is the Stefan Boltzmann constant, T is the temperature of the surfaces, and are the view factors for N number of surfaces. View Factor Formulas The view factor of a plate with itself or any other plate on the same face is equal to zero since the faces are all flat surfaces. The view factors are calculate d using equations from Modest [76] Figure 4 2 View Factors. A) P arallel plates, B) perpendicular plates, C) any two parallel plates, and D) any two perpendicular plates [76] PAGE 99 99 The view factors are numerically calculated by first coding the simple formulas for el, directly opposed rectangles see Figure 4 2 (A), ( 4 85 ) t an angle of 90 to each other see Figure 4 2 (C), ( 4 86 ) These two previous simple formulas are subroutines of the more complicated functions see Figure 4 2 (B), PAGE 100 100 ( 4 87 ) where see Figure 4 2 (D), ( 4 88 ) where Conduction Heat Transfer To calculate the conduction heat transfer through the wall of the thruster, the heat transfer equation was used ( 4 89 ) In this equation, c p is the specific heat capacity, is the density, and is the thermal conductivity of the material, while is the volumetric heat source. For the case examined, only the solution at steady state is important and there is no volumetric heat source in the equation, therefore, only the L aplace equation is solve d as ( 4 90 ) PAGE 101 101 Conduction Heat Transfer Boundary Conditions The boundary conditions for the conduction simulation can be either a fixed temperature ( Di richlet or first type) or heat flux via radiation Dirichlet is used everywhere the temperature is known. Heat flux via radiation boundary is used in places where there is heat transfer from the thruster wall s to the internal components of the small satellite assuming the temperature of the components of a small satellite is known and is quasi statically constant. This boundary is written as [77] ( 4 91 ) for multiple layers of thermal radiation insulation dividing the space between the thruster and the other components of the small satellite separated from the thruster by a small distance The multiple layers are a layer of aluminum surrounding the dielect ric material ( ) and a thin layer of aluminum ( ) in between the thruster and all other components. The components of the small satellite are assumed to be a black body ( ) for simplicity. PAGE 102 102 CHAPTER 5 NUMERICAL MODEL In this study an existing, previously described in literature, in house modular Multi scale Ionized Gas (MIG) flow solver platform dev eloped by Balagangadhar and Roy [78] [79] was used. The MIG flow solver platform ha s been utilized for many different applications using different finite element modules including electric propulsion, micro flows, nanoscale flows, fluid dynamics, and plasma phys ics [80] [81] [82] [83] For this dissertation, two modules of the MIG code have been utilized : a n existing finite element based rarefied gas module (RGM) and a recently developed finite difference based ionized gas module (IGM) They are loosely coupled in the following sequence. The RGM runs until convergence, producing the gas density ( ), components of velo city (u, v), and temperature (T). These variables from the RGM are passed to the IGM. Then, the IGM runs until convergence, producing the volumetric plasma induced electrostatic force components ( F x F y ) and the volumetric electro thermal heating source ( q ) needed by the new run of the RGM. The information obtained by each module is exchanged in this fashion until the steady state solution is obtained. Figure 5 1 shows the process. In Figure 5 1 the first step initializes the variables of both modules, including the initial voltage (V i ) value, based on the plasma breakdown voltage of the gas (~150 V for argon). Each module uses the L 2 norm of the solution to test for convergence. After both modules have reached convergence at a given voltage, the voltage is increased by a voltage difference ( ) and the process is repeated, until the maximum voltage (V max ) desired is achieved. For al l cases in this study, is 50 V. After the IGM reach es local convergence a test for global convergence is performed. PAGE 103 103 Figure 5 1 Schematic of the new MIG code with rarefied gas and ionized gas modules. The RGM and the IGM are run for five time steps even if their L 2 norms are less than their tolerances. If the L 2 norm of each module is less than its tolerance for all five consecutive time steps for three consecutive iterative loops between the RGM and the IGM, then glo bal convergence conditions are satisfied. The L 2 norms for the RGM and the IGM, using the given non dimensionalization explained in the following sections, are 10 4 and 10 10 for time steps of 10 8 and 5 10 13 seconds for each module, respectively. In t his chapter an overview of the numerical schemes and relevant information for each model will be provided The Ionized Gas Module The system of governing equations for the plasma the system of equations defined by ( 4 38 ) ( 4 39 ) and ( 4 51 ) is solved using the semi implicit, finite difference scheme module of the Multiscale Ionized Gas (MIG) code [78] [79] The equations ar e solved decoupled starting with the steady Poisson equation, then the ion continuity PAGE 104 104 equation, and the three metastable atoms continuity equations. The electron continuity and electron energy density equations are solved coupled with each other. Poisson e quation is approximated by combining the continuity equations of ions and electrons to predict the charge at present time step as ( 5 1 ) The Poisson equation is solved using second order central difference scheme where the right hand side is treated as a source ( ) since it depends on previous time step information. ( 5 2 ) In this equation, t he superscripts ( i, j ) and ( n ) refer to the grid spacing and time discretization The continuity equations for ions, electrons, and metastabl e atoms are solved using second order central difference o n a staggered mesh and first order implicit Euler metho d time discretization schemes for numerical stability ( 5 3 ) where the fluxes are solved using the S c harfetter Gummel flux discretization scheme and the sources ( ) are treated explicitly [84] The S c harfetter Gummel flux PAGE 105 105 discretization was devised to solve the drift diffusion equations for semiconductor applications [54] An upwind discretization of the flux is embedded in the method, causing monotonic solutions for large values of electric field. The S c harfetter Gummel flux discret ization assumes that the flux between two grid points is constant and the potential is linear. The value of the flux is found at the center between these two grid points. The ( k ) subscript in these equations refers to ions, electrons, or three types of met astable atoms The discretization of one of the fluxes needed is given as ( 5 4 ) In the previous equation u is the x component of the gas velocity since ( ) The electron energy density equation is solved in a similar manner, simultaneously with the electron continuity equation using a Newton Raphson non linear solver PAGE 106 106 ( 5 5 ) The fluxes are discretized in a similar fashion as in equation ( 5 4 ) giving the following form ( 5 6 ) The two equations can be written into a matrix system of the form ( 5 7 ) ( 5 8 ) The electron continuity equation in matrix form is constructed by [ K 11 ], { n e }, and { F e }. [ K 11 ] is the matrix that contains all t he factors that will be multiplied by vector { n e } to discretize the fluxes and the part of the time derivative that is unknown. { n e } is the vector to be solved that stores the values of n e { F e } is the vector that contains the source of the electron continui ty equation as well as the known part of the time derivative. Similarly, PAGE 107 107 the electron energy equation is form by [ [ K 21 ] [ K 22 ] ] { { n e } { n } } T and { F }. Together, the two equations form the new equation ( 5 8 ) The Newton Raphson non linear solver for this equation can be written as ( 5 9 ) where { n } has the values of n e and n at iteration p and { n }={ n } p+1 { n } p This system of equations is solved using a GMRES numerical matrix solver [85] The Rarefied Gas Module The rarefied gas is simulated using equations ( 4 71 ) to ( 4 74 ) modeled using finite element methods and loosely coupled with the ionized gas module as shown in Figure 5 1 The finite element method (FEM) first gained popularity in the field of [86] After the development of weighted residual criteria, it has found new applications in the fields of fluid mechanics and heat transfer [87] FEM is a popular numerical technique used for solving partial differential equations subdomains called elements. The solution within the element is constructed from the basis functions. This method has several advantages for solving transport problems such as simplicity for dealing with different meshes (e.g. triangular or quadrilateral) and order of elements (e.g. linear or quadratic), which makes it ideal for complex geometries [88] In addition, complex Neumann (flux) or Robin (convection) boundary conditions are rel atively easy to implement. PAGE 108 108 The numerical simulation of the rarefied gas code is performed using an existing, finite element based module in the Multi scale Ionized Gas (MIG) flow solver platform developed by Roy et al. [78] [79] This module of MIG utilizes the Galerkin weak statement combined with the Newton Raphson nonlinear solver [89] Bilinear elements are used in the RGC module for the numerical analysis of the RGEJ thruster. In order to provide stability to the solution, the streamline up winding (SU) artificial diffusion method in 2 D is used (explained in sub section : Verification of the RGM using plane Poiseuille flow for a simpler case) [90] The MIG flow solver platform had been utilized for many d ifferent applications, including electric propulsion, micro flows, nanoscale flows, fluid dynamics, and plasma physics [80] [81] [82] [83] Most recently, the finite difference Galerkin Weak Statement The equation system ( 4 71 ) to ( 4 74 ) for the rarefied gas can be written in a concise form as ( 5 10 ) ( 5 11 ) PAGE 109 109 ( 5 12 ) where q is the state variable, f is the kinetic flux vector, f v the dissipative flux vector and s is the source term. In this set of equations the specific heat at constant pressure ( c p ) is equal to the specific heat at constant volume ( c v ) minus the gas constant ( R ) The pressure can be set in terms of density and temperature using the ideal gas law and based on Stoke s hypothesis. Any real world smooth problem can be approximated as a Taylor or power series of known functions x j The Galerkin weak statement can be used to construct an approximation to the solution as a series of known spatial function multiplied by a set of unknown expansion coefficients [89] An approximation for the proble m at hand can be constructed as ( 5 13 ) where a i are unknown coefficients and are known fuctions of x j The Galerkin weak statement approach requires the approximation error to vanish in an overall integrated sense. The weak statement, the mathematical expression for minimizing the weighted residual over the domain is ( 5 14 ) PAGE 110 110 in this equation, defines the domain for the problem statement and w is the weight function set which is made identical to the to the corresponding trial function set for the approximation of the state variables. This guarantees that the associated approximation error is a minimum since the error is orthogonal to all basis functions spanning the space of possible Galerkin solutions. The final form of the wea k statement formulation becomes ( 5 15 ) When using the Galerkin weak statement, the differentiability requirement for the For example, the second order terms reduce to first order in the set of equations. Finite Element Basis Function As described by B aker et al. [89] the finite element base s are a set of polynomials defined uniformly on every subdivision ( finite element ) of the solution domain as created by placing nodes as desired for accuracy, and hence constructing the domain discretization denoted as The introduction of the discretization of is a central conc ept of finite element analysis, which greatly simplifies the construction of a wide range of highly suitable trial function sets as well as evaluation of the integrals in the weak statement. The finite element bases are defined as the set of functions ass ociated with the trail function ( ) that span a single generic element Chebyshev, Lagrange, or Hermite interpolation polynomials can be use d as the finite element basis ( N k ) to the k PAGE 111 111 degree depending on the problem statement for one, two, or three dimensions. In this study, Lagrange interpolation polynomials are used. The discrete approximation of the spatially discretized domain (the mesh) yields a union of elements as ( 5 16 ) Similarly, the integrated variables can be represented as the union of spatially and temporally discretized elements ( 5 17 ) ( 5 18 ) The spatially discretized two dimensional quadrilateral finite element basis definition yields ( 5 19 ) where the intrinsic coordinates spanning the quadrilateral constitute the tensor product system [68] Bilinear bases are chosen for this study (bilinear quadrilateral elements). PAGE 112 112 Bilinear Basis For the bilinear case ( k = 1) the basis { N 1 } must involve polynomials expressible in the global coordinate system x i with powers no higher than unity. The geometry of rectangular and straight sided q uadrilateral finite element domains is completely defined by the coordinates of the four intersections (vertices) of the boundary segment generators. Therefore, four expansion coefficients must be involved in the global coordinate equivalent of equation ( 5 19 ) The polynomial that involves x and y to powers no higher than u nity is ( 5 20 ) and the a i ( ) are the to be determined expansion coefficients [89] The fourth term is called the bilinear term since xy remains linear in both arguments. To establish the elements of which is a column matrix of order four in order to match the four entries that must occur in Q el the requirement is to equate equations ( 5 20 ) and ( 5 19 ) A B Figure 5 2 Generic rectangular parallelepiped element A) Element in g lobal coordinate system, and B) in transformed local coordinate system. PAGE 113 113 With a little insight, one can immediately write down the elements of by inspection of Figure 5 2 as [89] ( 5 21 ) The bilinear bases are not restricted to a rectangula r parallel epiped finite element. T hey can be used for any generic straight sided quadrilateral element with no particular orientation and no sides parallel to each other or to a global coordinate axis. The global element geometry bears no impact whatsoever on the fi nite element basis description. The main new difference is the emergence of a nontrivial coordinate transformation between ( x, y ) and ( 1 2 ). For all case s the basis provi de the following transformation ( 5 22 ) Therefore, the semi discrete finite element formulation for the prob lem statement is repre sented as the weak statement ( ) ( 5 23 ) PAGE 114 114 ( 5 24 ) Where use of the Green Gauss divergence theorem exposes the indicated boundary fluxes on The S el symbolizes is coefficients into the global arrays by doing a non overlapping summatio n over all the elements. The and i n the equation differential elements on and The d ifferentiation of the basis function depends on both t he global and local coordinates ( 5 25 ) The inverse coordinate transformation, from i to x j is required to evaluate the formation of the assembly matrices containing the convection and diffusion information. The transformation, called the Jacobian matrix is the invers e of the forward transformation ( 5 26 ) The weak statement naturally yields the surface integrals after the application of the Gre en Gauss theorem in equation ( 5 25 ) which contains the unknown boundary PAGE 115 115 fluxes wherever Neumann boundary conditions are enforced. The zero gradient boundary conditions are automatically enforced when the surface integral is removed. The numerical integration method of choice for finite element matrix evaluation is Gaussian quadrature, a variable accuracy procedure that exactly integrates polynomials on regula r regions. Specifically, a P th order Gaussian quadrature rule integrates exactly a polynomial of degree 2P 1 on in transform space, that is, the unit square as explained on page 219 of Baker et al. [89] For a linear problem statement and bilinear elements the product of all terms forming each element of the element matrix can be no more than cubic in i Thus, P = 2 Gaussian quadrature rule is appropriate since 2P 1 = 3 (cubic polynomial) in this case. For this study, the problem statement is nonlinear and a ( P = 3 ) Gaussian quadrature is used. Newton Raphson Scheme The semi discretized weak statement ( ) of equation ( 5 23 ) always yields an ordinary differential equation (OD E) system of the following form ( 5 27 ) where Q is the time dependent finite element nodal vector, M = S el ( M el ) is the mass matrix associated with element level interpolation, and carries the element convection information, the diffusion matrix, and all known data. A implicit time integration procedure is used (for all cases show n in this study = 1 .0 meaning that the implicit first order Euler time stepping method is employed) and the terminal ODE is usually solved using the Newton Raphson scheme: PAGE 116 116 ( 5 28 ) ( 5 29 ) ( 5 30 ) ( 5 31 ) ( 5 32 ) ( 5 33 ) w here the subscript n represents the time level and p the iteration index of the Newton Raphson scheme. Equation ( 5 28 ) uses the substitution of equation ( 5 27 ) solved for Equation ( 5 29 ) is found by moving all the terms of ( 5 28 ) to the left hand side but F is onl y equal to zero at steady state. Equation ( 5 30 ) is the Newton algorithm for ( 5 29 ) where the Newton Jacobian definition ( J ) given as ( 5 31 ) is obtained as a direct result of ( 5 30 ) Equation s ( 5 32 ) and ( 5 33 ) show the Newton Raphson scheme for the FEM and how the solution at the new sub iteration is recovered The solution is declared convergent when the L 2 norm s of the residual and solution for each of the state variables becomes smaller than a chosen convergence PAGE 117 117 criterion The criterion for the L 2 norms of the solution is given in the introduction of this chapter and the criterion for the L 2 norms of the r esidual is not enforced since it is not as restrictive as the L 2 norms of the solution for the cases studied. Decoupled Heat Transfer Codes Conduction Numerical Code The Laplace equation ( 4 90 ) is used to solve the heat transfer through conduction using the finite element module of the MIG platform where ( 5 34 ) ( 5 35 ) similarly to the rarefied gas module in the section labeled as The Rarefied Gas Module In this equation, the time derivative vanishes at steady state. Thermal Radiation Heat Transfer Numerical Code In order to model the thermal radiation heat los s equation ( 4 84 ) can be rewritten i n terms of vectors and matrices as ( 5 36 ) PAGE 118 118 ( 5 37 ) and solved by matrix inversion as ( 5 38 ) For the open surfaces of the rectangular prism form by the channel of the micro thruster an artificial surface is used to close the openings. Since any radiation leaving the cavity will not come back, barr ing any reflection from other surfaces nearby, the artificial surface must be black. In addition, the artificial surface is not emitting. These criteria are satisfied by making the artificial surface a black body at zero Kelvin [76] The external irradiation is assumed negligible to obtain the highest possible heat loss The matrix inversion in ( 5 38 ) is solved using MATLAB [91] PAGE 119 119 CHAPTER 6 BENCHMARKING AND VALIDATION I onized G as M odule (IGM) Verification and Grid Convergence Study of the IGM In order to v erify the plasma module the veracity of the model was tested using the method of manufactured solutions with the same procedure as Houba [54] A solution to a 1 D drift diffusion equation with variable transport coefficients is developed as described in this subsection. The test equation is ( 6 1 ) where V is the drift velocity equivalent to ZE x and D is the diffusion coefficient. Both are defined to vary in space using a polynomial expression ( 6 2 ) The constant C 8 is set to zero for convenience to ensure V (x=0) is equal to zero, which is useful for Neumann boundary conditions Two types of boundaries are tested, Neumann and Dirichlet boundaries, to test the discretization scheme. The boundaries are given as PAGE 120 120 ( 6 3 ) where A and B are constants. A simple form for n that satisfies these boundary conditions is ( 6 4 ) The solution is manufacture by solving for a source term S that will produce e quation ( 6 4 ) as the exact steady state solution to the governing equation ( 6 1 ) The source term obtained after substituting the solution for n in the governing equations is [54] ( 6 5 ) Two cases were tested for the verification. In these cases the magnitude of V is different by a factor of a thousand to test the behavior of the scheme in a region of high electric field, where the drift flux dominates over the diffusion flux. The coefficients for each case are given in Table 6 1 The accuracy of the solution was calculated using the L 1 norm defined as ( 6 6 ) where N is the number of nodes in the mesh. PAGE 121 121 Table 6 1 Coefficients for the grid refinement study [54] Coefficients Case 1 Case 2 A 3 3 B 4 4 C 1 1 1 C 2 1 1 C 3 1 1 C 4 1 1 C 5 2 2000 C 6 2 2000 C 7 3 3000 C 8 0 0 The following figure shows the analysis of the convergence of the method. Figure 6 1 Order of accuracy for Scharfetter Gummel flux discretization scheme for Case 1 and Case 2 The Scharfetter Gummel flux discretization scheme behaves as expected similar to the study done by Houba [54] If the drift velocity is low (low electric field), the diffusion flux dominates and the scheme is second order. The logarithm of the L 1 Norm increases as the logarithm of the distance between node increases with a slope of ~2.03 in a linear manner If the drift v elocity is high enough for the drift flux to dominate, PAGE 122 122 then the scheme decreases to first order accurate with a slope of ~1.16 An important finding obtained from the convergence analysis is that the order of accuracy of the Neumann boundary affects the co nvergence of the entire s imulation If the boundary discretization scheme is first order, the solution is first order even when the diffusion flux dominates. The Scharfetter Gummel flux discretization scheme is used in this dissertation for its computation al speed advantage and numerical stability. T he Neumann boundaries are discretized using a second order scheme. The continuity equations for all particles and the electron energy equation used in this dissertation have the same form as the equation tested in this grid convergence study; therefore by verifying the discretization scheme of a single equation with a generic form, the plasma module is verified. The Poisson equation is always second order and its discretization was also tested using the method of manufactured solution, but its analysis was not included since it reaches maximum convergence with a smaller number of nodes than the drift diffusion equations. Validation of the IGM Additionally, a simulation in 1 D of a parallel plate, capacitively co upled, low pressure, symmetric RF discharge driven at 13.56 MHz was performed and compared with Godyak et al [92] Although the RGEJ cases in this study use a DC applied potential difference, the RF discharge experiment of Godyak et al [92] is employed for this validation since the numerical models used for RF and DC glow discharges are almost identical except for the applied alternating voltage in RF discharges. In previous studies, Sitaraman a nd Raja [66] and Deconinck [11] have successfully used the LMEA PAGE 123 123 model to investigate an RF discharge thruster and a DC discharge thruster, respectively. Godyak et al [92] measured the discharge electrical characteristics (voltage, current, etc.) using argon (99.998% purity) at low pressure inside a glass cylinder with inner diameter of 14.3 (cm), cross section area = 160 (cm 2 ), and a discharge gap formed by the two para llel plate aluminum electrodes of 6.7 (cm). The discharge as low as 3.0 (mTorr) could be studied without overlapping electrode sheaths, and small enough that it could be considered as a 1 D discharge. For the cases shown, the pressure was 1.0 (Torr) and temperature of 300 K was assumed. The computational grid for all cases tested was composed of 671 nodes equally spaced. The secondary electron emission ( sec ) depends sens itively on surface conditions, morphology, impurities, and contamination; the commonly assumed value for ( sec ) for pure aluminum in a DC discharge is 0.1 [41] Since in Godyak et al [92] the purit y of the aluminum or its surface condition are not given and the discharge is not DC, several values of ( sec ) were tested to validate the code. The discharge voltage amplitude is applied to each electrode with equal magnitude but opposite phase. The peak to peak voltage is the total voltage of the simulation at peak value. Figure 6 2 shows the comparison between the experimental values and three different numerical simulations performed with different ( sec ) PAGE 124 124 Figure 6 2 Comparison between numerical model and experiment [92] The cases with ( sec ) equal to 0.01 matched the experimental results with the least error as shown in the following table Table 6 2 Comparison between experiment [92] and numerical model with ( sec =0.01 ) for (300 V < V Pk Pk < 800 V). Current Amplitude Experiment Data Voltage Amplitude Numerical Model ( sec =0.01 ) Voltage Amplitude Percent Error I A (A) V A (A) V A (A) Error (%) 1.229 169.2 167.5 1.03 1.726 233.0 225.0 3.44 2.366 317.6 300.0 5.53 2.795 353.5 350.0 1.00 3.214 388.1 400.0 3.07 A linear interpolation was used to interpolate between the experimental data points and to calculate the voltage amplitude at the given current amplitude obtained from the numerical model. The maximum percent error of the voltage amplitude in the PAGE 125 125 range of interest for the peak to peak voltage (300 to 800 V) is 5.53% for ( sec ) equal to 0.01, a reasonable error by drift diffusion model standards. Drift diffusion models typically have relatively large errors due to inaccuracies in the input coefficients as we ll variation in the reduced mobility when using different collision cross section libraries, in a Boltzmann solver such as B OLSIG +, for the electron mean energy (7 10 eV) region is of ~5% [44]. Rarefied G as M odule (RGM) Verification of the RGM using plane Poiseuille flow A verification of the RGM is performed using a simple case of plane Poiseu ille flow. Since the fluid dynamic solver used in the RGM is a multi physics code that could be used to solve a variety of problems, a plane Poiseuille flow is used to test the implementation of the same derivatives and artificial diffusion method that are used to solve the rarefied gas equations The incompres sible Navier Stokes equations are used for the plane Poiseuille flow case, which can be written in non dimensional form as [54] ( 6 7 ) ( 6 8 ) In an incompressible flow, the pressure waves m u st propagate at an infinite speed. This type of behavior is not a trivial task to numerically model. In order to solve PAGE 126 126 this problem, the continuity equation is modified by adding a time derivative of the pressure multiplied by the inverse of a compressibility factor ( ) which introduces a finite speed for the pressure disturbances to propagate through the domain. This approach is called the artificial compressibility method [93] and it produces a modified continuity equation given as ( 6 9 ) At steady state, the additional term vanishes and the original continuity equation is recovered. Addi tionally, the equations require an artificial diffusion mothed for numerical stability purposes. As explained in the section labeled as The Rarefied Gas Module the streamline up winding (SU) artificial diffusion method in 2 D is used [90] For th e case s in this dissertation where only rectangular, bilinear elements are used and the velocity components ( u 1 = u u 2 = v ) are defined in the same direction as the axe s ( x, y ), respectively, the artificial diffusion tensor ( ) is ( 6 10 ) where u is the velocity vector and PAGE 127 127 ( 6 11 ) The modified incompressible Navier Stokes equations are ( 6 12 ) ( 6 13 ) These modified equations approached the original equations at steady state with grid refinement. The system of equations are solved with the same approach explained in the section labeled as The Rarefied Gas Module where the equations can be written in a concise form as ( 6 14 ) ( 6 15 ) s imilarly to the rarefied gas equations PAGE 128 128 The plane Poiseuille flow case studied consist s of a channel form by two parallel plates with a pressure difference between its inlet and outlet. The following figure shows its geometry and the boundary conditions used in the numerical simulation Figure 6 3 Plane Poiseuille flow in a (L=10, H=1) channel with a (Re = 100) using a (2111) nodes mesh. All variables are non dimensional and only half of the domain is numerically simulated due to symmetry. For this simulation, Re and were set equal to 100. Three different meshes were tested: (2111), (4121), and (8141) nodes The test revealed the solution was mesh independent. Plane Poiseuille flow has an analytical solution assuming the flow is steady state and fully developed. The analyti cal solution is PAGE 129 129 ( 6 16 ) The comparison between the analytical and numerical solution s shows near perfect agreement For the x component of velocity ( u ), the L 1 norm is ~ 3.1 4 10 8 and the L norm is ~ 5.7 010 8 where the L 1 norm and the L norm are defined as ( 6 17 ) respectively, and N is the number of nodes of the mesh. The discrepancy between the analytical and the numerical solutions for the pressure ( P ) and the y component of velocity ( v ) are smaller than for u therefore they are not shown. Benchmarking /Validation of the RGM In this section, a benchmarking of the rarefied gas module is done for subsonic gas flows through a micro channels using results from Chen et al [94] which were validated within 1.15% accuracy with experimental results of Pong et al [95] The model assumes the gas flows through two parallel plates of length ( L = 3000 m) and width ( W PAGE 130 130 = 4 0 m) separated by a distance ( H = 1.2 m) The end effects are neglected and only the two dimensional geometry stretching in the x and y directions is considered. Table 6 3 Microchannel Properties of Fluid for Subsonic Gas Flows [94] Fluid P in /P out P out T in T w R (kPa) (K) (K) ( Pa s) (W/m K) (J/kg K) N 2 2.701 100.8 314 314 1.85 10 5 0.0259 1.4 296.7 In the previous table, T in is the inlet gas temperature and T w is the isothermal wall temperature. P in /P out is the ratio of inlet pressure vs. outlet pressure A ll other parameters and boundary conditions are given by Chen et al [94] The inlet and outlet pressures produce Knudsen numbers of 0.0217 and 0.0585, respectively. The flow in the microchannel is in the slip flow regime, which is the regime of inte rest for this study. The following table shows the grid dependence test done by Chen et al [94] Table 6 4 G rid dependence test of the center line u velocity (m/s) at different x locations [94] Grid x = 500 ( m ) x = 1000 (m) x = 1500 (m) x = 2000 (m) x = 2500 (m) 15007 0.4765853 0.5225323 0.5845432 0.6745498 0.8216915 300013 0.4760927 0.5222115 0.5845377 0.6749501 0.8229470 600023 0.4759963 0.5222144 0.5845423 0.6749584 0.8229740 The maximum discrepancy between Chen et al [94] and our results, in the u velocity at the centerline, occurs at x = 2500 m This maximum discrepancy of 1.6 % Table 6 5 Grid dependence t est done using MIG of the center line u velocity (m/s) at different x locations. Grid x = 500 (m) x = 1000 (m) x = 1500 (m) x = 2000 (m) x = 2500 (m) 15007 0.47601 0.52127 0.58127 0.66794 0.80982 300013 0.47606 0.52127 0.58126 0.66798 0.80984 600023 0.47620 0.52111 0.58124 0.66794 0.80988 PAGE 131 131 can be attributed to using different values of thermal conductivity and different numerical schemes. Chen et al [94] did not provide the thermal conductivity and explicit finite difference method was used to solve the governing equations. The pressure discrepancy along the centerline is similar to the u velocity discrepancy. Exact values for the pressure along the cen terline are not given in Chen et al [94] only a figure is provided, therefore an exact comparison of the pressure is not performed. The temperature remains near constant throughout the domain (~T w =314 K), and the v velocity i s close to zero. Overall, the rarefied gas module matched closely with results in literature. A comparison between the experimental results for the pressure at the centerline from Pong et al. [95] given in Chen et al. [94] are provided in the following table. Table 6 6 Comparison of centerline pressure ratios at different x locations [95] Pong et al. MIG (600023 Nodes) x/L P/P out P/P out Rel ative Error (%) 0.102841 2.559144 2.571570 0.485558 0.351046 2.243178 2.231138 0.536748 0.601134 1.831036 1.837706 0.364236 0.858915 1.376518 1.343269 2.415414 The maximum difference occurs closer to the exit of the channel, while Chen et al. [94] overshoots the solution by 1.15% the MIG code solution is lower than the experimental value by 2.4% This discrepancy could be due to the different viscosity 1 8540724 10 5 (Pa s) a nd thermal conductivity ( 0.02675 W/m K) at gas temperatures of ~ 314 .01 K found at the given data point in the experimental solution. PAGE 132 132 CHAPTER 7 CASES STUDIED Geometry and Grid The micro thruster was designed with a long (20 mm), narrow (3 mm) slot to prevent the possibility of catastrophically following figure. Figure 7 1 Geometry and design of the RGEJ thruster using several stack slots similarly to the FMMR thruster The absence of an expansion nozzle at the ex it of the channel is due to predicted low Reynolds numbers (< 100). In the limit of continuum isentropic flow through a large throat area [4] ( 7 1 ) where W and H are width and height. PAGE 133 133 The Reynolds number gives a measure of nozzle efficiency in terms of viscous [4] ( 7 2 ) the Reynolds number must remain constant or increase as the device is miniatur ized [4] Since small satellites require lower thrust and cannot operate at high enough plenum pressures, the operational Reynolds number for micro nozzles may decrease to values as low as 100, and as heat is added, the flow e xperiences a further decrease in Reynolds number. Micro thrusters with low throat Reynolds number (~100) do not experience any gains from an expansion nozzle [96] The low plenum pressure operation condition is chosen to scale the thrust and for the additional benefit of reduced propellant storage pressure, therefore easing the propellant tank mass and valve leakage requirements [4] Figure 7 2 Domain region of the RGEJ numerically simulated PAGE 134 134 Figure 7 2 shows the domain region numerically simulated (light blue). The IGM module models only the region inside the channel, (0 mm< x <20 mm), since the charged particle number densities are negligible at the exit plane. The mesh inside the channel has 401 31 nodes. For the cases tested, the anode is between (1 mm< x <2 mm) and the cathode (8 mm< x <19 mm). The RG M module models the channel region and the plume. The plume has 201 201 nodes and is 10 10 (mm). All cells are rectangular and have constant and y It is only necessary to solve one half of the domain due to symmetry. Cold gas thruster and constant thermal heating source thrusters results and comparison In order to understand the effect of gas heating in the RGEJ, results for adiabatic wall condition cases with different constant volumetric thermal heating source values = 1 .0, 2.0 3.0 and 4. 0 (MW/m 3 ) were obtained and compared with a cold gas thruster of the same design. The cold gas thruster simulation with ( = 0 .0 MW/m 3 ) is called the base case Figure 7 3 Regions of applied Only half of the domain is shown. The value of is constant over the given region. The f our regions in the domain where was applied are shown and label ed 1 4. PAGE 135 135 The performance parameters in this study are calculated at t he exit plane of the thruster using the fol lowing set of equations. ( 7 3 ) The following table shows the performance parameters of the different cases ferent regions shown in Figure 7 3 The total thermal heating source ( ) Four different were tested. These values are 150, 300, 450 and 600 (mW). The mass flow rate, given in SCCM, is calculated by the program and depends on the pressure difference between the inlet and outlet. PAGE 136 136 Table 7 1 Values of total thermal heating source ( ), mass flow rate ( ) thrust ( F Thrust ), shear force ( ) specific impulse ( I sp ), specific impulse increase ( I sp Inc.) compared to the base case exit plane Mach number ( M Exit ) at the centerline and exit Knudsen number ( Kn Exit ) at the centerline are displayed. The unit ( SCCM ) means cubic centimeter per minute at STP b Case Region of Applied as shown in Figure 7 3 F Thrust I sp a I sp Inc. M Exit Kn Exit (mW) (SCCM) (mN) (mN) (s) (%) Base case None 0 180.3 2.288 1.955 44.11 1.192 0.046 Case Q1 150 1 150 149.9 2.069 2.112 48.00 8.82 1.138 0.063 Case Q1 300 1 300 120.3 1.841 2.298 53.21 20.64 1.056 0.089 Case Q1 450 1 450 91.0 1.600 2.516 61.15 38.65 0.924 0.133 Case Q1 600 1 600 62.1 1.370 2.746 76.65 73.79 0.722 0.208 Case Q2 150 2 150 154.6 2.145 2.046 48.24 9.38 1.143 0.061 Case Q2 300 2 300 129.4 1.991 2.167 53.50 21.29 1.074 0.085 Case Q2 450 2 450 103.7 1.816 2.323 60.89 38.05 0.969 0.121 Case Q2 600 2 600 69.2 1.563 2.566 78.45 77.87 0.750 0.203 Case Q3 150 3 150 161.6 2.243 1.966 48.28 9.45 1.152 0.059 Case Q3 300 3 300 143.5 2.189 1.998 53.06 20.30 1.102 0.077 Case Q3 450 3 450 125.8 2.125 2.050 58.74 33.18 1.039 0.100 Case Q3 600 3 600 107.8 2.047 2.124 66.03 49.71 0.954 0.131 Case Q4 150 4 150 169.5 2.360 1.885 48.41 9.76 1.119 0.056 Case Q4 300 4 300 158.9 2.421 1.830 52.96 20.06 1.071 0.068 Case Q4 450 4 450 148.7 2.472 1.789 57.80 31.05 1.030 0.082 Case Q4 600 4 600 138.5 2.512 1.756 63.04 42.93 0.991 0.098 a The values shown were obtained using adiabatic conditions, not taking into account the heat loss through the walls due to conduction and radiation. b STP means standard temperature (273.15 K) and pressure (100000 Pa). Adding thermal energy to the gas increases the I sp for all cases. A decrease in the mass flow rate requirement occurs for all cases as is increased. In the cases where the thermal heating source is applied away from the exit plane, with regions of applied label ed 1, 2, or 3, the thrust also decreases due to an increase in the overall shear stress at the wall. Only in the cases where the heat source is placed near the exit PAGE 137 137 plane, at region 4, there is a gradual increase in thrust with increasing These set of behaviors for each region where was applied can be observed in the following figure. Figure 7 4 Graphical representation of different parameters. A) M ass flow rate, B) thrust, C) shear stress at the wall, D) specific impulse, E) thrust effectiveness and F) thruster total efficiency vs. total thermal heating source for each different region in Figure 7 3 using values from Table 7 1 In Figure 7 4 (A D), w hen the heat source is located near the exit plane, in region 4, we observe a linear decrease of mass flow rate require ment with a negative slope of 69.5 (SCCM/W) and a near linear increase in thrust with positive slope of 0.37 (mN/W) PAGE 138 138 This behavior in the thrust is due to an inversely proportional decrease in the shear force at the wall. When the heat source is located closer to the inlet, in region 2, we observe a decrease in the mass flow rate requirement with a slope of 1 82.0 (SCCM/W), and a decrease in thrust with a slope of 1.19 (mN/W). The I sp for all cases increases linearly while is less than or equal to 300 (mW), independently of the location where is applied. For greater values of with applied in regions 1 and 2, the I sp increases with a quadratic trend and the thruster is no longer choked ( M < 1), for cases with applied in regions 3 and 4 the linear positive trend is maintained. For the cases with equal to 600 (mW), the I sp increases the closer is applied to the inlet, with the exception of case Q1 600 with applied in region 1, due to the interaction between the applied and the inlet boundary condition. Case Q1 600 characteristics are most likely due to numerical effects. In Figure 7 4 (E F), the thrust effectiveness and total efficiency were calculated using the following equatio ns ( 7 4 ) ( 7 5 ) since the mass flow rate and thrust decrease for some cases with added a longer width ( W new ) is assumed for the given case to match the mass flow rate of the base PAGE 139 13 9 case The base case thrust is subtracted from the thrust of the given case to take into account only the additional thrust at a given mas s flow rate produced by adding thermal energy. For the cases shown in Table 7 1 P W is assumed to be equal to In Figure 7 4 (E), the thrust effectiveness decreases linearly for the cases with applied in either region 3 or 4 F or the cases with applied in region 1 or 2, the slope of the thrust effectiveness is negative and linear while the exit Mach number is (> 0.9), but for the cases with the highest with (M < 0.9), the sl ope of the thrust effectiveness does not follow the same negative, linear trend For the cases with the highest for all four regions t he value of the thrust effectiveness seems to correlate with the values of the exit temperature, which is the highest (~ 632 K) for Case Q2 600 The shift in the slope trend of the thrust effectiveness for cases with the highest applied in regions 1 or 2, indicates that for a thruster operating with low Mach num bers (M < 0.75) throughout the domain, the thrust effectiveness increases primarily with increasing exit temperature. In Figure 7 4 (F), t he total efficiency of the thruster increases with increasing and similarly to the thrust effectiveness, for the cases with the highest and low exit Mach numbers (M < 0.75), it depends strongly on the exit temp erature. The relatively low values of total efficiency are due to the under expansion of the flow since the exit pressure (~ 30 Pa) is higher than the external ambient pressure (~ 0.05 Pa) and only a fraction of the thermal energy is converted to kinetic e nergy. This is typical of non optimized micro electro thermal thrusters. The following figure show s two comparison s: between the base case and case Q2 600 and between the base case and case Q4 600 PAGE 140 140 Figure 7 5 Comparisons between base case and case Q2 600 and base case and case Q4 600 A, B) Density, C, D) x component of velocity E, F) temperature, and G, H) pressure The two comparison s between the base case and case Q2 600 and between the base case and case Q4 600 shown in Figure 7 5 (A H), revealed that adding thermal energy to the flow increases the temperature of the gas around the region whe re is added, and downstream of this region. An inversely proportional decrease in density happens where the temperature is increased in the domain and a distortion in the profile of pressure occurs in comparison to the base case in both cases, case Q2 600 and case Q4 600 The x component of velocity (or tangential velocity) displays a zone of PAGE 141 141 lower values than the base case before the thermal energy is added for both cases; followed by an acceleration zone after the thermal energ y is added. Figure 7 6 Comparison of base case case Q2 600 and case Q4 600 A) D ensity, B) temperature, C) x component of velocity, and D) pressure at the centerline. E ) shear stress at the wall, F ) Mach number at the centerline G ) thermal creep component of the shear stress at the wall, H ) mass flux at the wall, and I ) temperature at the wall Figure 7 6 show s a comparison of the base case case Q2 600 and case Q4 600 at two different cross sections : the wall and the centerline. If the thermal creep term is neglected in the slip flow boundary condition (for simplicity of the analysis), the shear PAGE 142 142 stress is proportional to u and T 0.5 at the wall. Since the mass flow rate requirement is smaller for cases with added thermal energy than for the base case the mass flux ( ) is smaller in most of the domain along the wall for cases with added as shown in Figure 7 6 (H) Alternatively, the shear stress shown in Figure 7 6 (E), increases with wall tempe rature shown in Figure 7 6 ( I ) These two competing contributions, ( ) and T cause the shear stress to decrease before the region where is added and increase right after. The thermal creep see Figure 7 6 (G), plays a minor role that increases this effect since molecules creep from cold towards hot regions [97] The thermal creep contribution to the shear stress is negative while T is increasing and positive while T is decreasing along the tangential direction at the wall. The shear force in case Q2 600 at the wall increased by about ~30% when compared to the base case For case Q4 600 the shear force is ~10% lower than the base case For both cases, the shear force causes an inversely proportional change in the thrust of similar percentage. The total shear force experienced by the fluid and the thrust produced by the device depends on the fraction of the wall area that is exposed to the higher temperatures. This ob servation shows that if maximizing the thrust in the device was the most desirable objective, the should be applied closer to the exit plane to minimize viscous losses. For case Q2 600 the most drastic change in comparison to t he base case due to the addition of thermal energy occurs in its density and temperature, shown in Figure 7 6 ( A ,B ) At the exit plane, T increases by a factor of (~3) and decreases inversely proportional. The increase in T at the exit plane, as observed in Figure 7 6 ( B ), increases the speed of sound. The Mach number at the exit plane, Figure 7 6 ( F ), decreases from 1.0 for the base case to 0.75 for case Q2 600 despite u increasing by 13% as shown PAGE 143 143 in Figure 7 6 (C) The addition of thermal energy has caused the flow to change from choked flow ( M Throat ~1) to subsonic ( M Throat < 1), meaning that viscous losses play an important role in this particular pressure regime for this channel design and prevent the exit tangential velocity from increasing with ( T 0.5 ) in a directly proportional manner. The pressure profile, Figure 7 6 ( D ), has a gradual slope until heat is added at ( x = 6 mm) and decreases with a steeper slope from this point on. Exit plane P is virtually unaffected, only 4% higher. At the exit plane, P and u change only a small amount in comparison to ( ) for case Q2 600 Given that ~ ( ) exit F Thrust ~ ( 2 + P ) exit and I sp ~ ( u + P/ ( ) ) exit the mass flow rate and thrust decrease as decreases T he I sp increases due to the ( P/ ( ) ) exit term increase, since ( P/ ( ) ) exit (1/ ) exit if P exit r emain s unchanged when compared to the base case which in these cases it does and the change in and u exit is small when compared with exit For case Q4 600 at the exit plane T increases by a factor of (~2), decreases by 40%, P and u increase by 33% and 23%, respectively, in comparison to the base case. Although the increase in P and u are beneficial to increase the I sp the I sp percent increase for this case is not as significant as for case Q2 600 The decrease in at the exit plane plays the most important role in increasing the I sp and it is inverse ly p roportional to the increase in T The base case has an exit Reynolds number of ~33, proving our initial assumption that an expansion nozzle would be counterproductive to increase the specific impulse since (Re < 100). All other cases have smaller exit Reynolds number due to heat addition. Increasing the tangential momentum accommodation coefficient in the slip flow boundary conditions from 0.89 to 1.0 increased the I sp for case Q2 450 by PAGE 144 144 4%. The thermal creep (transpiration) effects can affect the varia tion of pressure caused by tangential temperature gradients [97] The thermal creep was only significant for cases with very high temperatures (~ 900 K), such as case Q2 600 where the shear stress is affected the most in the region between (5 mm < x < 10 mm) by to the thermal creep as shown in Figure 7 6 (G) The higher order slip flow boundary condition presented b y Xue and Fan [98] as described in Zhang et al. [99] were tested for cases with low temperature, e.g. case Q2 450 and compared with cases using equation ( 4 82 ) and no effect in the I sp was detected, indicating that the first order approximation equation used in this study, equation ( 4 82 ) is sufficient to model the flow. The Kn at the centerline for the base case is 1.9 10 2 at the inlet and 4.5 10 2 at the exit plane, while for case Q2 600 the Kn is 1.8 10 2 at the inlet, an d as high as 2.0 10 1 at the exit plane. The Kn e xit for case Q2 600 is the highest of all cases and higher than the typical recommended range of values for slip flow regime (0.001 < Kn < 0.1). Maurer et al [100] estimated the upper limit of the slip flow regime as Kn = 0.3 0.1 [45], where Kn is based on the channel height as in th is dissertation For RGEJ, only the cases with ( Kn < 0.1 or Kn ~ 0.1) will be presented and studied. Based on th ese results region 2 wa s picked as the best location to apply the thermal heating source in order to increase the I sp The closer is applied to the inlet, the higher the T exit see case Q2 600 vs. case Q4 600 in Figure 7 6 ( B I ), which proportionally relates to the I sp Region 1 was not picked due to its interaction with the inlet boundary condition. The decrease in thrust experience by the cases with applied in region 2 could be counteracted by using a greater number of slots or a wider slot in the device if a given thrust is necessary. PAGE 145 145 RGEJ thruster performance The following table shows the performance parameters of several cases with plasma aided technology at different voltages with geometry, mesh, and boundary conditions as defined in Figure 7 2 Table 7 2 Values of: voltage (V), mass flow rate ( ) thrust ( F Thrust ), shear force ( ) specific impulse ( I sp ), specific impulse percent increase ( I sp Inc.), current ( I ), total thermal heating source ( ), total electrical power ( P W ), and fraction of the total electrical power converted into thermal heating source ( /P W ). Case V F Thrust I sp a I sp Inc. I P W / P W V SCCM mN mN s % mA mW mW % Base case 0 180.3 2.288 1.955 44.1 0.000 0 0 Case 450V 450 136.5 2.039 2.113 51.9 17.7 0.538 240 242 99.4 Case 550V 550 123.7 1.964 2.187 55.2 25.1 0.570 308 313 98.3 Case 650V 650 113.6 1.898 2.253 58.1 31.7 0.565 360 367 98.0 Case 750V 750 105.6 1.843 2.308 60.7 37.6 0.541 400 406 98.5 a The values shown were obtained using adiabatic conditions, not taking into account the heat loss through the walls due to conduction and radiation. In general, the plots of F Thrust , and I sp as functions of follow similar trends as the previous cases in Table 7 1 with constant thermal heating source applied in region 2. The current, thermal heating source and power vs. voltage are discussed in the section labeled as RGEJ thruster discharge characteristics In Table 7 1 the vs. I sp have positive slopes, 31.5 (s/W), approximately linear for the range of between 0 300 (mW) and independent of the location w here is applied. The plasma aided cases of RGEJ, shown in Table 7 2 have a significantly higher vs. I sp positive slope, 54.8 (s/W), which indicates that concentrating the total PAGE 146 146 value of in a smaller region of the domain is beneficial to increase the I sp for adiabatic cases. Case 750V has an I sp of 60.7 (s), a 37.6% improve ment over the base case The I sp of RGEJ, (60.7 s), operating at 750 (V) is 16% higher than the I sp of a highly optimized argon propellant cold gas thruster, (52 s), and 35% higher than the I sp of the argon propellant FMMR, (45 s) [4] The increase in I sp for case 750V over cold gas thrusters is achieved with only 406 (mW) per centimeter of width of the device and 98.5% of the total electrical power is converted to heating of the neutral gas This percentage is higher than the range (81 95%) predicted by Houba and Roy [101] for a device operating at an assumed constant temperature (300 K) and pressure (0.6 Torr) using air as the working fluid. The remainder of the discharge power goes into the electrons, which lose energy in inelastic collisions due to the various ionization, attachment, and excitation reactions [101] The thrust effectiveness of the device and total efficiency are shown in the followin g figure. Figure 7 7 Parameters of the RGEJ thruster for the given voltage regime. A) Thrust effectiveness and B) total efficiency. In Figure 7 7 (A), the thrust effectiveness decreases by 2.25% with a parabolic trend over the range of voltages tested (450 750 V) For higher voltages, the thrust PAGE 147 147 effectiveness is expected to increase b ased on the trend observed in Figure 7 4 (E) for In Figure 7 7 (B), the efficiency of the RGEJ thruster increases linearly with increasing voltage at a rate of 17.4 (%/kV) and it doubles over the range of voltages tested. The thrust effectiveness and tota l efficiency obtained for case 750V are 1240 ( N/W) and 10 %, respectively. For comparison, t he MPT operating at 100 Torr plenum pressure and wall temperatures of 300 K has a mass flow rate of 0.14 mg/s, a thrust of 67.4 N for the cold gas thruster case, and a thrust of 100 N for the 750 (V) and 650 (mW) plasma aided case which produces an I sp of ~ 74 (s) [10] [11] These values make the MPT thrust effectiveness and total efficiency equal to 50 ( N/W) and 0.6% respectively. The thrust effectiveness of the MPT is comparable to other electro thermal thrusters such as the VAT with a thrust effectiveness of 10 N/W [15] The thrust effectiven ess of the RGEJ is two orders of magnitude higher than other competing technologies. Figure 7 8 Left plots show a c omparison of the rarefied gas results between the base case ( on the bottom half ) and case 750V ( on the top half ). Right plots show a c omparison of the rarefied gas centerline results for the base case and all the plasma cases. A,B ) D ensity, C,D) tangential velocity E,F ) temperature and G,H) pressure Top wall anode (blue) and the cathode (red). PAGE 148 148 Figure 7 8 show s a comparison of the results of the rarefied gas simulations. Figure 7 8 ( E ) show s how the plasma locally heats the gas in case 750V to temperatures as high as (~640 K) near the corner of the cathode electrode where high electric field causes the volumetric electro thermal heating source to increase sharply. The localized heating cause s the gas density to decrease in the areas of high temperature, see Figure 7 8 ( A ) By comparison, the base case has a temperature profile that decreases along the x axis direction as we approach the exit plane due to the flow expansion and a similar density profile In Figure 7 8 ( F ), the exit temperature at the centerline is higher for the higher voltage cases with case 750V having a temperature of (~400 K), twice th e value of the base case (~200 K). The speed of sound doubles for case 750V The pressure profile, see Figure 7 8 ( G H ), is affected by the addition of thermal energy to the flow. The base case shows a near linear decrease in pressure for (x < 11 mm) with ( ) ~3.5 (Pa/mm), followed by a parabolic decrease. Four important effects are encountered in micro flows: rarefaction compressibility viscous heating and thermal creep [97] Out of those, compressibility and rarefaction are competing effects. The curvature in the pressure distribution found in channels with compressible flows is due to compressibility effects, the higher the Mach number the greater this effect becomes, the curvature increases as the inlet to outlet pressure ratio is increased [45]. Rarefaction decreases the curvature in the pressure distribut ion, which becomes increasingly linear as the free molecular flow regime is approached with increasing Kn numbers [97] In the pressure profile in Figure 7 8 ( H ), rarefaction is dominant in the base case for (x < 11 mm Kn = 0.03 ) and the compressibility effects become increasingly important in the rest of the domain due to increasing Mach number. For all PAGE 149 149 the plasma aided cas es, case 450V 750V the pressure profile is very similar independent of voltage (thermal energy input), but the pressure profile has two distinct linear regions, one before and one after the area where thermal energy is deposited. The effect of rarefaction dominates for (x < 16 mm Kn = 0.077 ) due to the higher Kn numbers found in the plasma aided cases in comparison to the base case For case 750V the two distinct linear regions have ( ) ~ 1.5 (Pa/mm) for (x < 7 mm) and 5.5 (Pa/mm ) for (7 mm < x < 16 mm), respectively. These two different regions are caused by the difference in temperature which affects the shear stress at the wall in the same manner as the cases in Table 7 1 The tangential velocity plots, in Figure 7 8 ( C D ), show the effect of having two distinct ( ) regions in case 750V The base case has a constant acceleration of the flow along the channel, but case 750V has approximately constant velocity before thermal energy is added due to the low ( ) for (x < 7 mm) follow by an acceleration region where ( ) is higher. The acceleration region is similar to a shorter ch annel operating with the same inlet to outlet pressure ratio at a higher inlet design in the future for optimization purposes. The addition of thermal energy increases the tangential velocity at the exit plane centerline from 311 (m/s) in the base case to 362 (m/s) in case 750V and de creases the density from 0.73 ( g/m 3 ) to 0.36 ( g/m 3 ), increasing the I sp of the thruster. Due to low Reynolds numbers (< 30) in all cases and d ominant viscous terms, no shock discontinuities are found. As the Reynolds number decreases with increasing PAGE 150 150 temperature (thermal energy input), the viscous losses increase causing a degradation of the thrust, which could be counteracted by extending the wi dth of the device. Figure 7 9 Shear stress at the wall for all cases in Table 7 2 The thermal creep a rarefaction effect, plays a significant role close to the cathode electrode corner at (~ 8 mm) where the magnitude of ( ) is the highest, causing a sharp discontinuity in the profile of shear stress at the wall. In this location, the shear stress abruptly increases due to the thermal creep but the thermal creep effect is negligible in the rest of the domain due to low ( ) values In the cases presented, the viscous heating effect is not apparent since for pressure driven compressible flows the expansion cooling negates it [97] Figure 7 10 show a comparison of the results of the ionized gas simulations for the given DC voltages. The applied volt age was varied from 450 to 750 V this range is within the operating conditi ons of the validation cases. In Figure 7 10 ( A, C, E,G, I, K ) display a comparison of the contours and ( B, D, F, H, J, L ) display a comparison of the centerline cross section, respectively, for the electrons number density, positive ions number density, the three types of metastable atoms number densit ies and electron energy density. PAGE 151 151 Figure 7 10 Comparison of plasma discharge results for case 450V and case 750V (on the left), and comparison of plasma discharge centerline results for all cases (on the right). A B ) E lectron number density in m 3 C, D) ion number density in m 3 E, F, G, H, I, J ) t hree types of metastable atom number densit ies in m 3 and K, L) electron energy density in J m 3 The top wall is shown to illustrate the position of the anode electrode (blue) and the cathode electrode (red). For case 750V the maximum values reach 4.03 10 16 (m 3 ), 4.23 10 16 (m 3 ), 1 18 10 1 5 (m 3 ), 1 09 10 1 5 (m 3 ), 1 55 10 1 4 (m 3 ), and 0.045 (J/m 3 ), respectively. The metastable atoms number densities are approximately one to two orders of magnitude lower than the charged particles and their contribution to ionization is very small for the given pressure operation regime. For case 4 50V the maximum number density values are very similar as in case 7 50V but the peak values shif t towards the inlet. PAGE 152 152 The electron temperature ( T e ), calculated using the electron energy and number densities, is (~4 eV) in the plasma column for all cases. In the cathode fall, where electron number density decreases to very small values ( <10 13 m 3 ) by comparison to the peak value (~10 16 m 3 ), the ( T e ) is over predicted and can increase exponentially in this region of vanishing electron densities due to a numerical artefact of the fluid model [10] This behavior of ( T e ) does no t affect the accuracy of the simulations for the other variables since the electron energy content is negligible in this part of the domain [10] Figure 7 11 Comparison of plasma discharge resul ts for case 450V and case 750V (on the left), and comparison of plasma discharge centerline results for all cases (on the right). A B ) E lectric potential in V, C, D) tangential component of the electric field in V/mm E, F ) charge separation in mC/ m 3 K, L) tangential component of the volumetric plasma induced electrostatic force in mN/ m 3 and I, J) the electro thermal volumetric heating source in W/ m 3 The top wall is shown to illustrate the position of the anode electrode (blue) and the cathode ele ctrode (red). PAGE 153 153 In Figure 7 11 (I, J ) show the volumetric electro thermal heating source concentrated in the cathode fall, reaching a maximum value of ~ 4 10 2 (W/ c m 3 ) near the corner of the cathode for case 750V Although this value looks large, it is concentrated in a small region of the domain and the integrated value of the total electro thermal heating source ( ) is ( 1 00) mW In Figure 7 11 (G, H ) the induced volumetric plasma induced electrostatic force is shown, which proportionally increases with the applied voltage, but its magnitude in the tangential direction is too insig nificant to contribute to the thrust of the device. An interesting finding in this comparison of cases at different voltages is that the cathode fall region is increased with increasing voltage while compressing the region of the plasma column. When case 4 50V and case 750V are compared in Figure 7 10 ( A, C, E, G, I, K ) and Figure 7 11 (A, E ) we can observe a shift towards the inlet of the charged particle number de nsities, electric potential peaks and charge separation. This phenomenon is observe d also in the corresponding centerline cross section plots across the range of voltages tested. I n Figure 7 10 (B, D), t he peak values of the number densities of charged particles remain relatively constant, ~ 4 10 16 (m 3 ) for cases with voltages of 550 to 750 V The peak value for case 450V is ~ 3 10 16 (m 3 ). For constant gas pressure and temperature cases, the number density of charged particles would drastically increase with increasing voltage, but for the plasma cases studied the discharge seems to be self limiting. For the cases sh own in Table 7 2 a t a given applied voltage between 550 and 650 V the current reaches a maximum and decreases with increasing voltage from this point on instead of increasing along with the applied voltage T he behavior of the PAGE 154 154 current is examined in the next section This phenomenon is caused by the decrease in the number density of neutrals as the gas temperature increases due to the localized h eating. T he normal ion flux ( ) the main contributor to the discharge current at the cathode, is higher for cases with higher voltages near the corner of the cathode electrode, with values ranging from 2.9 10 20 to 2.29 10 20 (m 2 s 1 ) for the given voltage range between 450 and 750 V. Conversely, is lower for the cases with higher voltages over most of the remainder of the cathode (8.45 mm< x < 13 mm ) and negligible for ( x >13 mm ). This means that cases wit h higher voltages have a higher current density near the corner of the cathode electrode, but operate at a lower discharge current since most of the area of the cathode electrode is exposed to lower values of The reduction of over most of the cathode (8.45 mm< x < 13 mm ) is due to the effects of the ion mobility and diffusion. T he ion mobility ( ) increases with a decrease in neutral number density or with a decrease in electri c field magnitude ( | E | ). For the range of voltages studied, the electric field magnitude remains surprisingly similar for (8.45 mm< x < 13 mm ) therefore, the ion mobility increases at the wall with increasing applied voltage due to the decrease in neutral number density which is due to the increase in gas temperature The ion diffusion ( ) is further increased by the increase in ion temperature which is assumed to be the temperature of the gas due to rapid thermalization of the io ns with other heavy particles. Since the ion diffusion is more sensitive to the effects of gas heating than the mobility, the diffusion term of the ion flux causes the normal ion flux ( ) and the ion number density ( ) to decrease over the given area of the PAGE 155 155 cathode (8.45 mm< x < 13 mm ) for cases with higher voltages. This behavior causes the discharge to be current self limiting and prevents the plasma from covering a greater area of the cathode electrode as the discharge voltage increases Figure 7 12 Case 750V A) Sum of the photon emitting reaction sources, and B) the ionization source for ions and electrons in their respective continuity equation, where ( S i = S e ). In Figure 7 12 (A), the sum of the photon emitting reaction sources is shown for case 750V This value gives the number density of all three t ypes of metastable atoms per time destroyed in processes that generate photons. Therefore, it is the production source term for photons and it reveals that most photons are produced in the center of the domain before the cathode fall in the same region whe re the peak values for metastable atoms occur. Figure 7 12 (A) provides a qualitative comparison between the numerical code results in this dissertation and the experimental results obtained for the RFET [14] working at low pressures (1.5 Torr) where the discharge is brightest at the center of the domain. In Figure 7 12 (B), m ost of the ionization happens close to edge of the plasma column before reaching the cathode fall and it is concentrated around the centerline PAGE 156 156 similar to the RFET [14] when it is worki ng at low pressures (1.5 Torr) The ions flow from this region of high electric potential, but a low electric field, to the walls where they recombined. In contrast, the electrons flow from the cathode to the rest of the walls, most of them flow into the anode and the rest into the dielectric sections of the walls to balance the ion current. The number of ions leaving the thruster through the exit plane is negligible; therefore their contribution to the thr ust is neglected. Some power is always lost through inelastic collisions to solid walls and the outflow, but the heating of the walls due to neutralization is not taken into account since there will always be some amount of thermal heating loss through con duction and radiation. O ur assumption of adiabatic walls is just an approximation to simplify the numerical simulation of a well insulated thruster. RGEJ thruster discharge characteristics The performance plasma discharge characteristics of RGEJ are presented in this section. The discharge current ( I ) is obtained by integrating the species current over the length of the electrode and the plasma force effectiveness ( ) is obtained using the tangential component of the plasma induced electrostatic force ( ) and the total electrical power ( P w ) [101] ( 7 6 ) ( 7 7 ) PAGE 157 157 ( 7 8 ) ( 7 9 ) w here t he tangential volumetric plasma induced electrostatic force ( ) is obtained from equation ( 4 69 ) and W = 1 cm. The following figure shows the performance plasma discharge characteristics of the RGEJ as functions of voltage Figure 7 13 Parameters of the RGEJ thruster for the given voltage regime. A ) T otal thermal heating source, total electrical power and discharge current as func tions of the applied voltage, B) specific impulse vs. voltage, and C) tangential component of the plasma induced electrostatic force and its effectiveness as functions of the applied voltage In Figure 7 13 (A ), the power consumed by the thruster varies from 242 to 406 (mW) over the range of voltages studied following a near linear trend, increasing with voltage. The current voltage distribution shows a peak at ~ 550 (V), the current decreases for higher voltag es. Deconinck et al [10] showed that current decreases in a thruster with constant gas wall temperatures when the wall temperature is increased at the cathode. This phenomenon is due to the decrease in ion flux at the cathode ca used PAGE 158 158 by the increase in gas temperature with increasing voltage as explained previously in the section labeled as RGEJ thruster performance This decre ase in current may be beneficial to decrease erosion of the cathode electrode permitting the thruster to operate at higher voltage and power. From 450 to 550 (V), the positive differential resistivity in the I V characteristics indicates that the thruster operates in the abnormal glow regime, for higher voltages the thruster is transitioning to the normal glow regime, where only a portion of the cathode is covered by the plasma. The I sp increases linearly with increasing voltage with a slope of 0.0292 (s/V ). For the voltage regime tested, the I sp has not saturated, meaning that it could be further improved by increasing the voltage. The t angential component of the plasma induced electrostatic force ( ) and the total t hermal heating source ( ) increase almost linearly with increasing voltage. Although the direct contribution of the force to the total thrust is negligible, the plasma force term near the cathode corner is about a tenth of the pressure gradient term in the Navier Stokes equations. F x and F y locally affect the flow by creating a low pressure region that pulls to flow towards the wall upstream of the corner of the cathode and accelerates the flow right after in the tangential direction. The highest plasma force produced in the tangential dire ction ( ) is 38 ( ), by comparison the plasma force produced by Houba and Roy [101] is 5 ( ) at 600 (V) and 0.6 (Torr) using air as the working fluid. The effectiveness of the plasma force increases with increasing voltage approaching a saturation point. PAGE 159 159 Thermal analysis A simple thermal analysis of the thruster was performed by taking the temperature distribution at the channel walls and exit walls for case 750V to calculate the hypothetic al heat loss through conduction and radiation. The analysis is performed with go verning, boundary equations and numerically modeled described in sections labeled Heat Transfer Governing Equations and Decoupled Heat Transfer Codes respectively The heat loss via conduction is a minimum for thruste that are stacked next to each other, similarly to FMMR In the following conduction analysis, the heat s slot. Figure 7 14 Condu ction heat loss analysis for Case 750V assuming 1 cm of width, with the temperature distribution given along the wall as T =f(x) In Figure 7 14 the conduction analysis of t he thruster is performed using a (401 101 ) mesh, assuming 1 0 mm thick, insulating walls made of silica aerogel with a thermal conductivity of 4.2 (mW/m K) for pressures (< 10 Torr). The thruster is covered PAGE 160 160 with an aluminum layer and an additional aluminum thin layer separates the thruster from the surroundings. For a well insulated satellite, a MEMS thruster system external average temperature is 285 K [102] The exit plate is th e external wall of the small satellite At steady state operation conditions, the rarefied gas module predicts an A thin gap separates the thruster walls from the e xit plate to prevent heat loss via conduction through that side of the channel walls. Neglecting the thermal resistivity of the electrode and dielectric materials due to their thinness, and neglecting any heat loss due to radiati on through the exit plate due to their assumed low emissivity, the heat loss through conduction of the internal walls at T =f(x) is 0.23 (mW) per centimeter of width of the thruster. Figure 7 15 Radiation heat loss analysis for Case 750V assuming 1 cm of width, with the temperature distribution of given along the wall as T =f(x). In Figure 7 15 the radiation analysis was performed u sing 400 plates in each internal wall, assuming each is a diffuse grey surface. The inlet plane is considered a diffuse grey surface with the emissivity of the electrodes since the plenum chamber will PAGE 161 161 contain micro machined pillars or a porous metal m ateri al to heat the propellant. The radiation analysis has shown that if micro machined pillars or a porous metal m aterial is not used in the inlet then a significant amount of thermal energy would scape the walls into the plenum, which would heat the gas in a non localize manner. The background radiation temperature of outer space is assumed to be 0.0 K. The emissivities of the dielectric and electrodes were assumed to be = 0.24 ( aluminum oxide ) and = 0.02 ( aluminum ), respectively [103] The heat loss due to radiation is 30.1 (mW). The radiation heat loss depends strongly on the emissivity of the internal walls. The adiabatic assumption in the cases presented is intended to provi de an upper heat loss of a thruster working in outer space is 30.34 (mW), which is 7.5 % of the input power in C ase 750V depending on the internal temperature of the small satellite the insulation layer of the thruster, and the emissivity of the internal walls. In laboratory conditions, with an environment temperature of 300 K, the heat loss through conduction would be negligible in a well insulated thruster and the h eat loss through radiation would be 23. 9 (mW), which is 5.9% of the input power, making the assumption of adiabatic walls a reasonable approximation for a comparisons with experiments It is important to indicate that 1.5% of the total power in C ase 750V i s not transferred to the gas and it is mostly los t to the walls in the form of heat (98.5% of the total power is transferred to the gas) This loss is expected to decrease the discrepancy between the adiabatic assumption cases and future experiment results even further. PAGE 162 162 CHAPTER 8 SUMMARY Summary of the Results The RGEJ d evice shows encouraging results. For a single channel with a width of 1 cm operating at the highest voltage tested of 750 V, the RGEJ thruster requires a current of 0.541 mA and a total electrical power of 406 mW to heat the flow t o temperatures as high as ~640 K. Under the given operation conditions, the RGEJ produces a thrust of 1.843 mN at a mass flow rate of 105.6 SCCM (3.095 mg/s), which produces a specific impulse (I sp ) of 60.7 s. This I sp is an improvement of 37.6% over a col d gas thrusters with the same geometry and working parameters, and it is an improvement of 16% and 35% over optimized argon propellant cold gas thrusters and argon propellant FMMR respectively Based on the I sp vs. voltage characteristics, the thruster co uld operate at an even higher voltage to further increase the I sp The thrust effectiveness and total efficiency obtained for the highest voltage tested are 1240 ( N/W) and 10%, respectively. For comparison, t he MPT operating at 100 Torr plenum pressure an d wall temperatures of 300 K has a mass flow rate of 0.14 mg/s, a thrust of 67.4 N for the cold gas thruster case, and a thrust of 100 N for the 750 V and 650 mW plasma aided case which produces an I sp of ~74 s [10] [11] These 50 ( N/W) and 0.6%, respectively. The Vacuum Arc Thruster (VAT) has a thrust effectiveness of 10 N/W [15] The thrust effectiveness of the RGEJ is two ord ers of magnitude higher than the s e similar plasma aided technologies. For the configuration of electrodes studied and range of voltages, the majority of the energy is co nverted into gas heating (~98%). This perce ntage is higher than the PAGE 163 163 range (81 95%) predicted by Houba and Roy [101] for a device operating at an assumed constant temperature (300 K) and pressure (0.6 Torr) using air as the working fluid. T he plasma force plays a negligible role in the injection of momentum and the RGEJ operates as an electro thermal thruster (heating the gas) instead of operating as an electrostatic thruster (accelerating the flow by ion collisions) High heat loss is typical in MEMS scale micro t hrusters [4] For example, FMMR required twice (200% of) the power use d to heat the gas to operate due to heat losses. In contrast, a simple thermal analysis estimated that the proposed RGEJ thruster would require 7.5 % more in put power to maintain the same performance due to heat loss in outer space, this number could be decreased by reducing heat loss through conduction emissivity. T he peak values of the number densities of charged parti cles remain relatively constant ( ~ 4 10 16 m 3 ) for cases with voltages of 550 to 750 V For constant gas pressure and temperature cases, the number density of charged particles would drastically increase with increasing voltage, but for the plasma cases studied, the discharge is self limiting due to the increase in gas temperature as the applied voltage is increased. T he current voltage distribution shows a peak at ~ 550 V and the current decreases in cases with higher voltages. From 450 to 550 (V), the positive differential resistivity in the I V characteristics indicates that the thruster operates in the abnormal glow regime, for higher voltages the thruster is transitioning to the normal glow reg ime, where only a portion of the cathode is covered by the plasma. This decrease in current may be beneficial to decrease erosion of the cathode electrode due to sputtering, PAGE 164 164 permitting the thruster to operate at higher voltage and power. The sputtering rat e at the front corner of the cathode was calculated to be 3.68 10 8 cm/s for copper electrodes. The cathode corner region experiencing peak ion current density will be eroded 1 mm every 755 hours of continuous thruster operation. This sputtering rate can b e decrease d by using stronger materials in the corner of the cathode (e.g. tungsten). The RGEJ material compatibility between the propellant and the surface when operat ing with argon is not optimum because argon causes high sputtering, but if RGEJ can ope rate with water propellant, the long aluminum cathode electrodes could be silver plated to help mitigate the formation of hydrogen peroxide, serving as a catalyst to decompose it. Due to the low plenum pressure and the long slot design, RGEJ is expected to have very low or negligible valve leakage problems and no passage clogging complications that could result in a single point failure. The system reliability and durability has to be determined, but the only real concern would be the electrode erosion due to sputtering, which is the most damaging in the front corner of the cathode electrode. Due to the RGEJ design simplicity, with no moving parts, other reliability and durability problems are avoided. The RGEJ can even op erate as a cold gas thruster The RG EJ produces thrust in the (mN) necessary for slew maneuvers at very low power The minimum impulse bit of the RGEJ should be less (~ 1/2) than for the Moog cold gas thruster ( model 58 125A [19] ) with a minimum impulse bit of 100 N s, since the minimum impulse bit depends on the thrust ( Moog: F Thrust maximum open/close response time. If liquid or solid propellant is used, the weight and storage density of the propulsion system could be greatly reduce d due to a reduction in the tank mass [4] T he PAGE 165 165 storage density of the propellant is important to minimize the volume required for propellant tanks. The integration complexity of the system is similar to FMMR [4] consisting of a propellant tank, drain/fill valve, filter assemblies, control valve, and the thruster itself. The RGEJ design can be developed using simple MEMS fabrication techniques and common materials, which results in low cost batch f abrication. The efficiency of the RGEJ is low (10%) but higher than other electro thermal thrusters using plasma aided technology (e.g. MPT). The working principle of the RGEJ has been tested and the device has shown promising results that indicate its p otential for real world applications as a micro propulsion thruster for small satellites Further investigation is need ed to optimize and improve the thruster. The Numerical Model Loosely coupling the finite element based rarefied gas module and the finite difference based ionized gas module in the MIG framework resulted in a stable approach to solve internal, slip flow problems with glow discharges. These type of problems are highl y unstable since the plasma and gas interact strongly with each other, but their time scales are widely different by (10 5 s). The approach used in this study circumnavigates this problem by increasing the voltage in small increments and by having severa l conver gence criteria as explained in the introduction of Chapter 5 Using finite difference to solve the ionized gas is faster than using finite volume for rectangular geometries. While using the finite element based rarefied gas module will allow the analysis of more complex geometries in the future. The MIG numerical platform is computationally cost effective and it was able to solve each individual case PAGE 166 166 in a single processor under a period of a month using a mesh with 52801 nodes in the RGM to model the channel and the plume and a mesh with 12431 nodes in the IGM to model the channel Future Directions to Continue Improving RGEJ Argon was selected as the working fluid because it is the noble gas of choice for benchmarking pla sma numerical codes due to the few reactions necessary to model the glow discharge at low pressure. However, thermal thruster developers prefer to use gases with a lower molecular weight and a higher gas constant. Our goal is to develop a phase change thru ster concept that uses liquid or solid propellant, instead of argon, to avoid heavy storage tank and valve leakage problems. The low minimum required operating pressure is selected by design to eventually develop this class of thruster. In future simulatio ns, a different gas with better suited properties for the RGEJ thruster should be use d and the geometry of the thruster shoul d be optimized. A suitable chemistry model should be used to simulate the thruster using the new propellants, which inevitably will require a greater number of species equations. As part of the proposed future work, the parallelization of the ionized gas module would be a necessity in order to perform these new simulations. The numerical model could be improved by coupling a thermal a nalysis module with the existing rarefied and ionized gas modules; it could also be improved by building and testing a prototype of the thruster for validation purposes. PAGE 167 167 LIST OF REFERENCES [1] P. Romero et al., "Optimal station keeping for geostationary satellites with electric propulsion systems under eclipse constraints," Mathematics in Industry, vol. 12, p. 260, 2008. [2] D. M. Goebel and I. Katz, Fundarmentals of Electric Propulsion: Ion and Hall Thrusters, Hoboken, NJ: John Wiley & Sons, 2008. [3] D. J. Barnhart, T. Vladimirova and M. N. Sweeting, "Very Small Satellite Design for Distributed Space Missions," Journal of Spacecrafts and Rockets, vol. 44, no. 6, pp. 1294 1306, 2007. [4] M. M. Micci and A. D. Ketsdever, Micropropulsion for Small Spacecraft, Reston, VA: AIAA, 2000. [5] M. Leomanni, A. Garulli, A. Giannitrapani and F. Scortecci, "Propulsion options for very low Earth orbit microsatellites," Acta Astronautica, vol. 133, pp. 444 454, 2017. [6] W. Marshall and C. Boshuizen, "Planet Labs' remote sensing system," in Proceedings of the 27th AIAA/USU Conference on Small Satellites Logan, U.S.A., 2013. [7] A. Ketsdever, A. Green, E. P. Muntz and S. Vargo, "Fabrication and Testing of the Free Molecule Micro Reistojet: Initial Results," in 36th Joint Propulsion Conference & Exhibit Huntsville, AL, 2000. [8] Soni, "Characterization of Plasma Actuator Base d Microthruster Concepts for High Altitude Aircrafts and Cubsats," University of Florida, Gainesville, FL, 2014. [9] S. Roy, "Method and Apparatus of Small Satellite Propulsion". US Patent US20120204618, 6 December 2012. [10] T. Deconinck, S. Mahadevan and L. L. Raja, "Computational simulation of coupled nonequilibrium discharge and compressible flow phenomena in a microplasma thruster," Journal of Applied Physics, vol. 106, p. 063305, 2009. [11] T. D. Deconinck, "Simulation Studies of Direct Current Microdischarges for Electric Propulsion," Ph.D. Dissertation, Austin, Tx: The University of Texas at Austin, (2008). PAGE 168 168 [12] A. Greig, C. Charles and R. Boswell, "Simulation of main plasma parameters o f a cylindrical asymmetric capacitively coupled plasma micro thruster using computational fluid dynamics," Frontiers in Physics vol. 2, no. 80, pp. 1 9, 2015. [13] T. Takahashi, Y. Takao, K. Eriguchi and K. Ono, "Numerical and experimental study of mi crowave excited microplasma and micronozzle flow for a microplasma thruster," Physics of Plasmas, vol. 16, no. 8, pp. 1 14, 2009. [14] C. Charles and R. Boswell, "Measurement and modelling of a radiofrequency micro thruster," Plasma Sources Science and Technology, vol. 21, no. 2, pp. 1 4, 2012. [15] J. Mueller, R. Hofer and J. Ziemer, "Survey of propulsion technologies applicable to cubesats," in 57th JANNAF Propulsion Meeting Pasadena, CA, 2010. [16] J. Mueller, "Thruster Options for Microspacecraft: A Review and Evaluation of Existing Hardware and Emerging technologies.," in Proceedings of the 33rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference Seattle, WA, 1997. [17] S. W. Janson, H. Helvajian, W W. Hansen and J. Lodmell, "Microthrusters for Nanosatellites," in The Second International Conference on Integrated Micro Nanotechnology for Space Applications Pasadena, CA, 1999. [18] D. B. Scharfe and A. Ketsdever, "A Review of High Thrust, High De lta V Options for Microsatellite Missions," in 45th AIAA Joint Propulsion Conference & Exhibit Denver, CO, 2009. [19] J. Mueller, J. K. Ziemer, R. R. Hofer, R. E. Wirz and T. O'Donnell, "A Survey of Micro Thrust Propulsion Options for Microspacecraft a nd Formation Flying Missions," in In: 5th Annual CubeSat Developers Workshop San Luis Obispo, CA, 2008. [20] J. M. Parker, D. Thunnissen, J. Blandino and G. Ganapathi, "The Prelimininary Design and Status of a Hydrazine MilliNewton Thruster Development ," in 35th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Los Angeles, CA, 1999. [21] A. London, A. Epstein and J. Kerrebrock, "A High Pressure Bipropellant Micro Rocket Engine," in 36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Huntsville, AL, 2000. [22] B. Marcu, G. Prueger, A. Epstein and S. Jacobson, "The Commoditization of Space Propulsion: Modular Propulsion Based on MEMS Technology," in 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference Tucson, AZ, 2005. PAGE 169 169 [23] A. D. Ketsdever, R. H. Lee and T. C. Lilly, "Performance testing of a microfabricated propulsion system for nanosatellites applications," Journal of Micromechanics and Microengineering, vol. 15, no. 12, pp. 2254 2263, 2005. [24] R. Lee, A. Bauer, M. Killingsworth, T. Lilly, J. Duncan and A. Ketsdever, "Performance Characterization of a Free Molecule Micro Resistojet Utilizing Water Propellant," in 43th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Cincinnati, OH, 2007. [25] F. Rysanek, J. Hartmann, J. Schein and R. Binder, "MicroVacuum Arc Thruster Design for a CubeSat Class Satellite," in 16th Anual AIAA/ISU Conference on Small Satellites Logan, UT, 2002. [26] J. Schein, M. Krishnan, R. Shotwell and J. Ziemer, "Vacuum Ar c Thruster for Optical Communications Missions," in 39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Huntsville, AL, 2003. [27] R. Burton, J. G. Eden, L. Raja, S. J. Park, J. Laystrom Woodard and M. de Chadenedes, "Thrust and Efficiency Performance of the Microcavity Discharge Thruster," University of Illinois at Urbana Champaign, Arlington, VA: Air Force Office of Scientific Research, 2011. [28] R. Wirz, M. Gale, J. Mueller and C. Marrese, "Miniature Ion Thrusers for Precision Formation Flying," i n 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit Ft. Lauderdale, FL, 2004. [29] 2.5 A New Optimized Microthruster of Giessen University," in 31st International Electr ic Propulsion Conference University of Michigan, MI, 2009. [30] V. Khayms and M. Martinez Sanchez, "Fifty Wall Hall Thruster for Microsatellites," in Micropropulsion for Small Spacecraft (Progress in Astronautics and Aeronautics vol. 187, Reston,VA, AIAA, 2000, p. Chap. 9. [31] K. Polzin, Y. Raites, J. Gayoso and N. Fisch, "Comparisons in Performance of Electromagnet and Permanent Magnet Cylindrical Hall Thrusters," in 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit Nashville, TN, 20 10. [32] K. Polzin, T. Markusic, B. Stanojev, A. Dehoyos, Y. Raitses, A. Smirnov and N. Fisch, "Performance of a Low Power Cylindrical Hall Thruster," Journal of Propulsion and Power, vol. 23, no. 4, pp. 886 888, 2007. PAGE 170 170 [33] A. Smirnov, Y. Raitses and N. Fisch, "Performance Studies of Miniaturized Cylindrical and Annular Hall Thrusters," in 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Indianapolis,IN, 2002. [34] N. Demmons, V. Hruby, D. Spence, T. Roy, E. Ehrbar and J. Zwahlen, "ST7 DRS Miss ion Colloid Thruster Development," in 44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Hartford, CT, 2008. [35] V. e. a. Hruby, "ST7 DRS Colloid Thruster Systerm Development and Performance Summary," in 44th AIAA/ASME/SAE/ASEE Joint Propulsion Confe rence Hartford, CT, 2008. [36] D. Simon, B. Land and J. Emhoff, "Experimental Evaluation fo a Micro Liquid Pulsed Plasma Thruster Concept," in 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit Sacramento, CA, 2006. [37] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, Second Edition, New YorK: Plenum Press., 1984. [38] I. Langmuir, "Oscillations in Ionized Gases," National Academy of Science, vol. 14, no. 8, pp. 627 637, 1928. [39] M. Mitchner and C H. Kruger, Partially Ionized Gases, New York: Wiley, 1973. [40] A. von Engel, Ionized Gases, Pittsford, NY: American Institute of Physics, 1997. [41] M. A. Lieberman and A. J. Lichtenberg, Principles of Plasma Discharge and Materials Processing, Ne w York: John Wiley & Sons, 1994. [42] M. J. Druyvesteyn and F. M. Penning, "The mechanism of electrical discharges in gases of low pressure," Reviews of Modern Physics, vol. 12, no. 2, pp. 87 174, 1940. [43] Y. P. Raizer, Gas Discharge Physics, New Y ork: Springer, 1991. [44] J. R. Roth, Industrial Plasma Engineering, Volume 1: Principles, Bristol, United Kingdom: Institude of Physics Publishing, 1995. [45] M. A. Laughton and D. D. Warne, Electrical Engineer's Reference Book 16th ed., Oxford: Newnes, 2002. [46] A. Zumbilev, "Carbonitriding of Materials in Low Temperature Plasma," in Heat Treatment Conventional and Novel Applications InTech, 2012, p. 420. PAGE 171 171 [47] R. K. Marcus, Glow Discharge Spectroscopies, New York: Plenum Press, 1993. [48] M. J. Madou, Fundamentals of Microfabrication: the Science of Miniaturization, 2nd Edition, New York: CRC Press, 2002. [49] H. Conrads and M. Schmidt, "Plasma gen eration and plasma sources," Plasma Sources Science and Technology, vol. 9, no. 4, p. 441, 2000. [50] A. A. Fridman, Plasma Chemistry, Cambridge University Press, 2008. [51] A. M. Howatson, An Introduction to Gas Discharges, 2nd edn., Oxford: Pergamon Press, 1976. [52] R. Mason and P. Melanie, "Sputtering in a glow discharge ion source pressure dependence: theory and experiment," Journal of Physics D: Applied Physics, vol 27, no. 11, pp. 2363 2371, 1994. [53] V. E. Golant, A. P. Zhilinsky and I. E. Sakharov, Fundamentals of Plasma Physics, New York: John Wiley & Sons, 1980. [54] T. Houba, "Parallel 3 D Numerical Simulation of Dielectric Barrier Discharge Plasma Actuators," University of Florida, Gainesville, FL, 2016. [55] A. Blanco and S. Roy, "Rarefied gas electro jet (RGEJ) micro thruster for space propulsion," Journal of Physics D: Applied Physics, vol. 50, no. 455201, pp. 1 15, 2017. [56] N. Balcon, "At mospheric pressure Radio Frequency discharges, diagnostic and numerical modeling," University of Toulouse, Toulouse, France, 2007. [57] G. J. M. Hagelar and L. C. Pitchford, "Solving the Boltzmann equation to obtain electron transport coefficients and ra te coefficients for fluid model," Plasma Sources Science and Technology, vol. 14, pp. 722 733, 2005. [58] J. C. Maxwell, "A dynamical theory of the electromagnetic field," Philosophical Transactions of the Royal Society, vol. 155, pp. 459 512, 1865. [59] E. A. Rafatov, D. Akbar and S. Bilikmen, "Physics Letters A," Modelling of non uniform DC driven glow discharge in argon gas, vol. 367, pp. 114 119, 2007. [60] R. Kee, G. Dixon Lewis, J. Warnatz, M. Coltrin and J. Miller, "A fortran computer code package for the evaluation of gas phase multicomponent transport properties," Sandia, 1995. PAGE 172 172 [61] G. M. Petrov and C. M. Ferreira, "Numerical modeling of the constriction of the dc positive column in rare gases," Physical Review E, vol. 59, no. 3, 1998. [62] S. F. Biagi, "Biagi v7.1 (Magboltz version 7.1)," June 2004. [Online]. Available: www.lxcat.net. [Accessed 30 January 2015]. [63] G. J. M. Hagelaar, F. J. de Hoog and G. M. W. Kroesen, "Boundary conditions in fluid models of gas discharges," Physi cal Review E, vol. 62, no. 1, pp. 1452 1454, 2000. [64] S. T. Surzhikov and J. S. Shang, "Two component plasma model for two dimensional glow discharge in magnetic field," Journal of Computational Physics, vol. 199, pp. 437 464, 2004. [65] D. J. Griffiths, Introduction to Electrodynamics, Third Edition, Upper Saddle River, NJ: Prentice Hall, 1999. [66] H. Sitaraman and L. L. Raja, "Simulation studies of RF excited micro cavity discharge for micro propulsion applications," Journal of Physi cs D: Applied Physics, vol. 45, no. 18, pp. 1 11, 2012. [67] J. P. Boeuf, L. C. Pitchford and K. H. Schoenbach, "Predicted properties of microhollow cathode discharges in xenon," Applied Physics Letters, vol. 86, no. 7, pp. 1 3, 2005. [68] R. Raju, "Hydrodynamic Model for Investigation of Gas Flows Through Micro Geometries and Nanopores," M.S. Dissertation, Mechanical Engineer Dept., Flint, MI: Kettering University, (2003). [69] H. Hanley, "The Viscosity and Thermal Conductivity Coefficie nts of Dilute Argon, Krypton, and Xenon," J. Phys. Chem. Ref. Data, vol. 2, pp. 619 642, (1973). [70] R. J. Goldston and P. H. Rutherford, Introduction to Plasma Phyiscs, Bristol and Phyladelphia: Plasma Physics Lab., Princeton University, Intitute of P hysics Publishing, (1995), p. 150. [71] A. R. Choudhuri, The Physics of Fluids and Plasmas: an Introduction for Astrophysicists, Cambridge, United Kingdom: Cambridge University, Cambridge University Press., (1998), p. 19. [72] J. C. Maxwell, "On Stresses in Rarefied Gases Arising from Inequalities of Temperature," Philosophical Transactions of the Royal Society Part1, vol. 170, pp. 231 256, (1879). PAGE 173 173 [73] Annalen der Physik und Chemi, vol. 64, pp. 101 130, (1898). [74] A. Agrawal and S. V. Prabhu, "Survey on measurement of tangential momentum accommodation coefficient," Journal of Vacuum Science and Technology A, vol. 26, no. 4, pp. 634 645, 2008. [75] W. M. Tro tt et al., "Experimental Measurement of Thermal Accommodation Coefficients for Microscale Gas Phase Heat Transfer," AIAA, no. 2007 4049, 2007. [76] M. F. Modest, Radiative Heat Transfer 2nd Ed., Amsterdam: Academic Press, 1993. [77] F. P. Incropera, D. P. Dewitt, T. L. Bergman and A. S. Lavine, Fundamentals of Heat and Mass Transfer 6th Ed., Hoboken, NJ: John Wiley & Sons, 2007. [78] D. Balagangadhar and S. Roy, "Design Sensitivity Analysis and Optimization of Steady Fluid Thermal Systems," Computer Methods in Applied Mechanics and Engineering, vol. 190, pp. 5465 5479, (2001). [79] S. Roy, "Combining Galerkin Matrix Perturbation with Taylor Weak Statement Algorithms," Computer Method in Applied Mechanics and Engineering, vol. 184, pp. 87 9 8, (2000). [80] S. Roy and B. P. Pandey, "Development of a Finite Element Based Hall Thruster Model," Journal of Propulsion and Power, vol. 19, pp. 964 971, 2003. [81] S. Roy, R. Raju, H. Chuang, B. Cruden and M. Meyyappan, "Modeling gas flow through microchannels and nanopores," Journal of Applied Phyisics, vol. 93, pp. 4870 4879, (2003). [82] S. Roy and D. Gaintonde, "Force interaction of high pressure glow discharge with fluid flow for active separtion control," Physics of Plasmas, vol. 13, no. 023503, 2006. [83] H. Kumar and S. Roy, "Multidimensional hydrodynamic plasma wall model for collisional plasma discharges with and without magnetic field effects," Physics of Plasmas, vol. 12, no. 093508, 2005. [84] D. Scharfetter and H Gummel, "Large Signal Analysis of a Silicon Read Diode Oscillator," IEEE Trans. Elec. Dev., vol. 16, pp. 64 77, (1969). PAGE 174 174 [85] Y. Saad and H. S. Martin, "GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systerms," SIAM J. Sci. Stat. Comput., vol. 7, no. 3, pp. 856 869, 1986. [86] J. Oden and E. R. A. Oliveira, Lectures on Finite Element Methods in Continuum Mechanics, (A Collection of Lectures), Huntsville: University of Alabama in Huntsville Press, 1976. [87] A. J. Baker, Finite Element Computational Fluid Mechanics, New York: Hemisphere Publishing Corporation, 1983. [88] C. Wang, Numerical Modeling of Microscale Plasma Actuators, Gainesville: University of Florida, 2009. [89] A. J. Baker and D. W. Peppe r, Finite Elements 1 2 3, McGraw Hill, 1991. [90] J. Donea and A. Huerta, Finite Element Methods for Flow Problems, Chichester, England: John Wiley & Sons, 2003. [91] The MathWorks, Inc., MATLAB, Natick, Massachusetts, United States, 2017. [92] V. A. Godyak, R. B. Piajak and B. M. Alexandrovich, "Electrical Characteristics of Parallel Plate RF Discharges in Argon," IEEE Transactions on Plasma Science, vol. 19, pp. 660 676, 1991. [93] A. J. Chorin, "A numerical method for solving incompressible viscous flow problems," Journal of Computational Physics, vol. 2, pp. 12 26, 1967. [94] C. S. Chen, S. M. Lee and J. D. Sheu, "Numerical Analysis of Gas flow in Microchannels," Numerical Heat Transfer, Part A, vol. 33, pp. 749 762, 1998. [95] K. C. Pong, C. Ho, J. Liu and Y. Tai, "Non Linear Pressure Distribution in Uniform Microchannels," ASME PUBLICATIONS FED, vol. 197, p. 51, 1994. [96] A. D. Ketsdever, M. T. Clabough, S. F. Gimelshein and A. A. Alexeenko, "Experimental and numerical determinat ion of micropropulsion device efficiencies at low Reynolds number," AIAA Journal, vol. 43, no. 1, pp. 633 641, 2005. [97] G. E. Karniadakis, A. Beskok and N. Aluru, Microflows and Nanoflows: Fundamentals and Simulation. Vol. 29, Springer Science & Business Media, 2006. PAGE 175 175 [98] H. Xue and Q. Fan, "A high order modification on the analytic solution of microchannel gasesous flow s.," in In: Proceeding of ASME fluids engineering division summer meeting Boston, USA, 2000. [99] W. M. Zhang, G. Meng and X. Wei, "A review on slip models for gas microflows," Microfluids and nanofluidics, vol. 13, no. 6, pp. 845 882, 2012. [100] J. Maurer, P. Tabeling, P. Joseph and H. Willaime, "Second order slip laws in microchannels for helium and nitrogen," Physics of Fluids, vol. 15, no. 9, pp. 2613 2621, 2003. [101] T. Houba and S. Roy, "Numerical study of low pressure air plasma in an ac tuated channel," Journal of Applied Physics, vol. 118, no. 23, pp. 1 9, 2015. [102] S. Kang and H. Oh, "On Orbit Thermal Design and Validation of 1 U Standardized CubeSat of STEP Cube Lab," International Journal of Aerospace Engineering, vol. 2016, pp. 1 17, 2016. [103] J. Henninger, "Solar Absorptance and Thermal Emittance of Some Common Spacecraft Thermal Control Coatings," National Aeronautics and Space Administration, Washington DC, 1984. [104] T. Trudel, S. Bilen and M. Micci, "Design and Perfo rmance testing of a 1 cm Miniature Radio Frequency Ion Thruster," in 31st International Electric Propulsion Conference University of Michigan, MI, 2009. PAGE 176 176 BIOGRAPHICAL SKETCH Ariel Blanco conducted his undergraduate studies at the University of Florida, where he graduated with a dual Bachelor degree in m echan ical and a erospace e ngineering i n December 2008 with the honors of Cum Laude. After graduation, Ariel Blanco decided to continue his education at the same university He obtained a Master of Science degree in Aerospace Engineering i n May 2010. During the last semester of his Master of Science degree, a fter careful consideration, he realize d that he wanted to study electric propulsion in order to co and joint the Applied Physics research group (APRG) under the tutelage of Dr. Subrata Roy and started his doctoral research in the field of micro electric propulsion. He received his Doctor of Philosophy d egree from the University of Florida in the fall of 201 7. |