LIBRARY Michigan State University This is to certify that the dissertation entitled TRANSPORT PHENOMENA OF FLOW THROUGH HELIUM AND NITROGEN PLASMAS IN MICROWAVE ELECTROTHERMAL THRUSTERS presented by Scott Stanley Haraburda has been accepted towards fulfillment of the requirements for Ph.D. degreeinghfimimljngineering Major professor Date ‘70 ”E. / «1%0700/ MSU is an Affirmative Action/Equal Opportunity Institution 0- 12771 PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 cJCIRCIDdoDuepss-pJS OF FLOW THROUGH HELIUM AND NITROGEN _. WAVB EECTROTHERMAL THRUSTERS By Scott Stanley Haraburda A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY lpepartment of Chemical Engineering .1901 A g - ”Dimer Martin "I, Flaw.- .c; ABSTRACT TRANSPORT PHENONIENA OF FLOW THROUGH HELIUM AND NITROGEN PLASMAS IN MICROWAVE ELECTROTHERMAL THRUSTERS By Scott Stanley Haraburda Electric rocket thrusters have effectively been demonstrated for uses in deep space and platform station keeping applications. However, the operational thruster lifetime can significantly decrease as the electrodes erode in the presence of the propellant. The Microwave Electrothermal Thruster (MET) would be an alternative propulsion system that would eliminate the electrode altogether. In this type of thruster, the electric power would be transferred from a microwave frequency power source, via electromagnetic energy, to the electrons in the plasma sustained in the propellant. The thrust from the engine would be generated as the heated propellant expands through a nozzle. Diagnostic methods, such as spectroscopic, calorimetric, and photographic methods using the TM011 and TM012 modes in the microwave resonant cavity, have been used to study the plasma. Using these experimental results, we have expanded our understanding of plasma phenomena and of designing an operational MET system. As a result, a theoretical and computational based model was designed to model the plasma, fluid, and radiation transport phenomena within this system using a helium and nitrogen mixture based propellant. Additionally, a literature search was conducted to initially develop potential non—propulsive applications of microwave generated plasma systems. 7" ~ ., _ ~ W il‘yibo'uise, Jessica Allyson, and Christine Frances, who have endured ". a... c , WMMW me while I was writing this dissertation. . wWea was - My, for his support of my research at West Point. Appreciation is given to hacienda: in the chemical engineering department, at Bayer Corporation, at General Electric Plastics, and the United States Army for their eating assistance in helping me complete this dissertation. This research was supported in part by fully-funded schooling from the United States Army under the provisions of Army Regulation 621-1 and by grants from the National Mounties and Space Administration — Lewis Research Center. (Iii-EYE} ?~ _ -__ ._-n ‘1‘ l a-I'“ .n-rr ._ :_. J,,.. l .— m'W—n h n D . n.- that TABLE OF CONTENTS LIST OF TABLES ........................................................................... xi LIST OF FIGURES ........................................................................... xiii NOMENCLATURE ........................................................................ xviii CHAPTER 1 Introduction ................................................................. 1 1.1 Problem Description ............................................. 1 1.2 Problem Significance ............................................. 3 1.3 Research Objectives ............................................... 4 1.4 Problem Solving Approach ....................................... 6 1.5 Experimental System ............................................. 9 1.5.1 Microwave Cavity ....................................... 9 1.5.2 Plasma Containment ..................................... 12 1.5.3 Flow System .............................................. 12 1.5.4 Microwave Power ........................................ 12 1.5.5 Temperature Probes ............................. ' ......... 15 1.5.6 Spectroscopy .............................................. 15 1.6 Experimental Results ............................................. 17 CHAPTERZ BACKGROUND ........................................................... 23 2.1 Plasma Properties and Applications ............................ 23 2.2 Electrothermal Propulsion ....................................... 24 2.3 Microwave Induction ............................................. 26 2.4 Research Direction ................................................ 27 CHAPTER 3 FLOW THROUGH A MICROWAVE GENERATED PLASMA 29 3.1 3.2 3.3 3.4 3.5 Basis of Model ..................................................... 29 Model Development Using Separate Transport Processes .. 31 Development of Computational Technique ................... 32 Examination of Parameters ...................................... 32 Past Research Review ............................................ 33 3.5.1 Experimental .............................................. 33 3.5.2 Theoretical ................................................ 34 3.5.2.1 Chapman ................................. 34 3.5.2.2 Morin .................................... 34 3.5.2.3 Haraburda ............................... 35 CHAPTER 4 MODEL DEVELOPMENT .............................................. 36 4.1 Overview ........................................................... 36 4.2 Realbody Radiation (Section 1) ................................. 38 4.3 Outer Chamber (Section 2) ...................................... 44 4.4 Coolant Chamber (Section 3) ................................... 47 4.5 Discharge Section (Section 4) .................................. 49 4.6 Validity of Assumptions ......................................... 55 4.6.1 Realbody Radiation ...................................... 55 4.6.2 Outer Chamber ........................................... 55 4.6.3 Coolant Chamber ......................................... 55 4.6.4 Discharge Section ........................................ 56 4.6.5 NASA Program Simulation ............................. 57 4.6.6 Summary ................................................... 57 4.7 Summary of Model Equations ................................... 57 CHAPTER 5 MODEL SIMULATIONS ................................................ 59 5.1 General .............................................................. 59 5.2 Realbody Radiation ................................................ 59 5.3 Outer Chamber ..................................................... 62 5.4 Coolant Chamber .................................................. 77 5.5 Discharge Section ................................................. 86 5.6 NASA Program Simulations ..................................... 97 5.6.1 Pressure Changes ......................................... 103 5.6.2 Energy Changes .......................................... 103 5.6.3 Nitrogen Mixtures ....................................... 107 5.7 Comparison with Experimental Results ....................... 109 CHAPTER 6 SCALE-UP ANALYSIS ................................................. 111 6.1 Introduction ........................................................ 1 1 1 6.2 Scale—up Issues .................................................... 1 1 1 6.2.1 Operability ................................................ 11 1 6.2.2 Maintainability ........................................... 112 6.2.3 Controllability ............................................ 112 6.2.4 Cost ........................................................ 112 6.2.5 Schedule ................................................... 113 6.2.6 Performance ............................................... 113 6.2.7 Public Acceptance ........................................ 114 6.2.8 Six Sigma .................................................. 114 6.2.8.1 DMAIC ........................................... 117 6.2.8.1.1 Design .......................... 118 6.2.8.1.2 Measure ........................ 118 6.2.8.1.3 Analyze ........................ 119 6.2.8.1.4 Improve ........................ 119 6.2.8.1.5 Control ......................... 120 6.2.8.2 DFSS ............................................. 120 6.2.8.2.1 Define ........................... 122 6.2.8.2.2 Scope ........................... 122 6.2.8.2.3 Analyze ......................... 122 6.2.8.2.4 Design .......................... 123 6.2.8.2.5 Implement ...................... 123 6.2.8.2.6 Control .......................... 123 6.3 Summary ........................................................... 124 CHAPTER 7 NON-PROPULSIVE APPLICATIONS ................................ 125 7.1 Detoxification of Hazardous Materials ........................ 125 7.2 Surface Treatment of Commercial Materials ................. 126 7.3 Novel Methods in Chemical Reaction Procedures ........... 128 CHAPTER 8 CONCLUSIONS .......................................................... 129 CHAPTER 9 RECOMIVIENDATIONS ................................................ 132 9.1 Model Development .............................................. 132 9.2 Material Development ............................................ 133 9.3 Advanced Computer Simulations ............................... 134 9.4 Flight Simulations ................................................. 134 9.5 Alternative Applications ......................................... 135 9.6 Scale-up ............................................................ 136 9.7 Summarized List of Recommendations ........................ 136 REFERENCES ............................................................................... 138 APPENDIX A: FORTRAN PROGRAMS ............................................... 151 A.1 Gauss Elimination ................................................ 152 A2 Curve-Fitting ....................................................... 156 A3 Outer Chamber ..................................................... 161 A.4 Coolant Chamber .................................................. 165 A5 Discharge Section ................................................. 169 A6 Statistical Mechanics ............................................. 175 A.7 Electromagnetic Field ............................................. 188 A8 Chemical Kinetics ................................................. 190 APPENDIX B: ATOMIC ENERGY LEVELS ......................................... 196 B] Helium .............................................................. 197 B2 Nitrogen ............................................................ 199 APPENDIX C: A TWO-DIMEN SION AL KINETICS PROGRAM SIMULATION ................................................................................ 203 vii APPENDIX D: PLASMA TRANSPORT PHENOMENA ............................. 220 DJ Collisional Processes ............................................. 220 D.1.1 Types ...................................................... 220 D.1.2 Cross-Section ............................................. 223 D.1.3 Coulomb Forces .......................................... 223 D2 Charged Particle Motion in BM. Field ......................... 225 D.2.l E.M. Field (TM012) ...................................... 225 D.2.2 Equations of Motion ..................................... 226 D.2.3 Power Absorption ........................................ 228 D3 Distribution Function ............................................. 229 D.3.1 Statistical Mechanics .................................... 229 D.3.l.1 Partition Functions ..................... 230 D.3.1.2 Chemical Equilibrium Reactions 233 D.3.1.3 Species Mole Fraction ................ 233 D.3.1.4 Average Molecular Weight ........... 235 D.3.1.5 Compressibility Factor ................ 235 D.3. l .6 Plasma Density ......................... 236 D.3.1.7 Energy / Enthalpy ...................... 236 D.3.1.8 Entropy .................................. 237 D.3. 1.9 Chemical Potential ..................... 238 D.3.1.10 Heat Capacity ........................... 239 D.3.2 Boltzmann Equation ..................................... 240 D.3.3 Conservation Equations ................................. 242 D.3.3.1 Continuity Equation ................... 242 D.3.3.2 Momentum Equation .................. 242 D.3.3.3 Energy Equation ........................ 243 D.3.4 Collisional Processes .................................... 243 D.3.4.1 Neutral-Neutral ......................... 243 D.3.4.2 Neutral-Ion .............................. 245 D.3.4.3 N eutral—Electron ........................ 246 D.3.4.4 Charged Particles ....................... 246 D.3.4.5 Excited Species ......................... 247 D.3.4.6 Collision Integral and Rates ........... 247 D.3.5 Transport Coefficient .................................... 248 D.3.5.1 Electrical Conductivity ................ 249 D.3.5.2 Thermal Conductivity ................. 250 D.3.5.3 Mobility ................................. 250 D.3.5.4 Viscosity ................................ 251 D.3.5.5 Diffusion Coefficient ................. 251 D4 Chemical Kinetics ................................................ 256 D.4.1 Reaction Rate ............................................. 256 D.4.2 Reaction Time to Equilibrium .......................... 259 viii APPENDIX E: FLUID TRANSPORT PHENOMENA ................................ 262 El Flow Functions .................................................... 262 E2 Conservation Laws ................................................ 264 E.2.1 Continuity ................................................. 264 E22 Motion ..................................................... 265 E.2.3 Energy ...................................................... 266 E3 Compressible Fluid Flow Variables ............................ 267 E4 Method of Characteristics ........................................ 268 APPENDIX F: RADIATION TRANSPORT PHENOMENA ....................... 272 El Blackbody .......................................................... 273 F2 Graybody ........................................................... 279 APPENDIX G: COMPUTATIONAL METHODS .................................... 283 G.1 Algebraic Sets of Equations ...................................... 283 G.1.1 Linear Equations ......................................... 283 G.1.2 Non-linear Equations .................................... 285 6.2 Data Curve-Fitting ................................................ 287 GB Classification of Partial Differential Equations ............... 288 G4 Galerkin Method .................................................. 289 G5 Finite-Difference .................................................. 291 G6 Finite-Element ..................................................... 293 G7 Analysis of NASA's TDK Computer Program ................ 294 0.7.] Assumptions .............................................. 295 G.7.2 Thermodynamic Data .................................... 295 6.7.3 Nozzle Geometry ......................................... 296 APPENDIX H: COMPUTATIONAL PARAMETERS .............................. 298 H.1 Thermodynamic Properties ...................................... 298 H.1.l Mole Fraction ............................................. 298 H.1.2 Molecular Weight ........................................ 298 H.1.3 Compressibility ........................................... 301 H.1.4 Plasma Density ............................................ 301 H.15 Electron Density .......................................... 304 H.1.5 Enthalpy ................................................... 304 H.1.6 Entropy ..................................................... 304 H.1.7 Heat Capacity ............................................. 304 H2 Transport Coefficients ............................................ 309 H.2.1 Electrical Conductivity .......... . ....................... 311 H.2.2 Thermal Conductivity ................................... 313 H.2.3 Viscosity .................................................. 316 H.2.4 Diffusion .................................................. 316 H.3 Chemical Kinetics ................................................ 316 H.3.l Reaction Rates ............................................ 316 H.3.2 Reaction Time to Equilibrium .......................... 321 ‘4, ' . ' . [,3pefimental Values 324 Researcthnstraints 324 TWCoefficimts 330 I H WNW 90W? in :. I Experimcmr‘ i': Rein-2:54.51 : . " ' Cw: (33:2 "- T ‘Polymnm: Car-”rim: - " «11:. . .:-.;...‘ 5.3 5.4 6.1 7.1 El 6.] G2 H.l Electron Temperature Range Erato ' v . ' Ex" ' tal Power Distributions .............................. 19 Experimental Plasma Dimensions .............................. 20 Research Parameters ............................................. 30 Radiation Media Description ................................... 39 Research Model Equations ...................................... 58 Emissivity Values for Selected Materials ..................... 61 Outer Chamber Simulation List ................................ 65 Coolant Chamber Simulation List ............................. 81 Simulation vs. Experimental Values .......................... 1'10 Sigma Significance ............................................... 115 Current Plasma Detoxification Systems ....................... 126 Radiation Model Parameters .................................... 272 Thermodynamic Coefficients (NASA Program) .............. 296 Nozzle Geometry Parameters .................................... 297 Chebyschev Polynomial Coefficients .......................... 309 Coefficients for Different Polynomial Orders ................. 311 Polynomial Coefficients for Transport Coeffients ............ 313 q 4 v y V! ‘5 ye IA (i cf. (Cp/RTo) .. Mn Conductivitymholcm) ...... ............ ' i ' imminent-104m sec/m2) Wmmiumamomavityuethrg/cmsewq ..... .5514 lipflirDifl'usion Coefficient (cm2/sec) '\ .zirn'. ‘ “ i’ 15'." .- R;K.tri’h‘ -~ ROSES . ' 0mm t. lttlr t . “ a s.‘ Chaim 1‘3 Wt 327 328 329 329 330 331 331 332 Wm . 10 Wye Cavity ................................................. 11 fl» Plasma Containment Tubes ...... 13 M ‘ ' Microwave Power Source ........................................ 14 w Spectrompy System .. ...... 16 15"“ Calorimetry System ................................................ 18 ' 11? Thruster ............................................................. 26 i . 23* Discharge Properties ......... . . . .. 28 | 4-1' Microwave Plasma Model Overview ........................... 37 43 Plasma Surface Temperature .................................... 39 43" Radiation Environment Media .................................. 40 43‘ Radiation Emission Sketch ....................................... 41 ‘95“ Radiation Energy Balance ........................................ 43 iii Radiation Plasma Surface Temperature ........................ 60 ‘. 5:35 OuterChamber Sketch 63 - ii‘i 34M ammrmmm (Outer Chamber) ...................... 66 fl va- »~ . . ‘ 67 “Mmumcmmsxa grid) as r ‘ I . “VET“:‘M: he Mama-d) Dunn-enn- ”s - ‘ 7 l 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 5.20 5.21 5.22 5.23 5.24 5.25 5.26 5.27 5.28 5.29 530 Temperature Gradient (Outer Chamber, 9x9 grid) ............ 71 Temperature Gradient (Outer Chamber, 11,11 grid) ......... 72 Temperature Gradient (Outer Chamber, 27x3 grid) 73 Temperature Gradient (Outer Chamber, 3x27 grid) . . . 74 Temperature Gradient (Outer Chamber, 9x9p1 grid) ........ 75 Temperature Gradient (Outer Chamber, 9x9p2 grid) ........ 76 Temperature Gradient (Coolant Chamber Sketch) ........... 79 Boundary Temperature (Coolant Chamber ................... 81 Temperature Gradient (Coolant Chamber, 3x3 grid) ........ 82 Temperature Gradient (Coolant Chamber, 9x9 grid) ........ 83 Temperature Gradient (Coolant Chamber, 9x9p1 grid) 84 Temperature Gradient (Coolant Chamber, 9x9p2 grid) 85 Temperature Gradient (Discharge Section, 400 torr) ........ 89 Temperature Gradient (Discharge Section, 600 torr) ........ 90 Temperature Gradient (Discharge Section, 800 torr) 91 Temperature Gradient (Discharge Section, 600 torr, mix) .. 92 Velocity Gradient (Discharge Section, 400 torr) .............. 93 Velocity Gradient (Discharge Section, 600 torr) .............. 94 Velocity Gradient (Discharge Section, 800 torr) .............. 95 Velocity Gradient (Discharge Section, 600 torr, mix) 96 Electron Density Gradient (Disch. Sect, 400 torr) ........... 98 Electron Density Gradient (Disch. Sect, 400 torr) ........... 99 Electron Density Gradient (Disch. Sect, 600 torr) ........... 100 xiv D“ \H De. I.| .- 5.31 5.32 5.33 5.34 5.35 5.36 5.37 5.38 5.39 5.40 5.41 5.42 6.1 6.2 6.3 7.1 9.1 D.l D2 03 D4 E1 E2 Electron Density Gradient (Disch. Sect, 800 torr) ........... 101 Electron Density Gradient (Disch. Sect, 600 torr, mix) 102 Specific Impulse Pressure Plot .................................. 103 Helium Mole Fraction Gradient (5 Energy Levels) ........... 104 Temperature Nozzle Gradient .................................... 105 Pressure Nozzle Gradient ......................................... 105 Mach Number Nozzle Gradient ................................. 106 Specific Impulse Nozzle Gradient .............................. 106 Helium Mole Fraction Gradient (5 Mixture Levels) ........ 107 Electron Mole Fraction Gradient .............................. 108 Specific Impulse Mixture Plot ................................... 108 Mach Number Mixture Plot ..................................... 109 Six Sigma Process Improvement Graph ....................... 116 DMAIC Process ................................................... 117 DFSS Process ...................................................... 121 Plasma Surface Application Sketches .......................... 127 Discharge Chamber Cross—Section .............................. 135 Collisions ........................................................... 22] Classical Potential Plot ........................................... 225 Collision Path ...................................................... 249 Ambipoiar Diffusion .............................................. 254 Characteristics ..................................................... 269 Characteristic Nozzle ............................................. 271 XV 'n P: P. ’.—'._‘ F.l F2 'F.3 E4 F.5 G.1 H1 H.2 H3 H4 H5 H6 H7 l-I.8 H.9 H. 10 H.11 H. 12 H13 H. 14 H.15 H.16 H.17 Radiation Heat Transfer Mode] ................................. 273 Blackbody Radiation Energy .................................... 275 Radiation Model Schematic ..................................... 278 Graybody Radiation Emissivity ................................. 281 Graybody Radiation Energy ..................................... 281 ODE Nozzle Geometry .......................................... 297 Helium Mole Fraction Plot (0.1 ATM) ........................ 299 Helium-Nitrogen Mole Fraction Plot (0.1 ATM, 50% mix) 299 Nitrogen Mole Fraction Plot (0.1 ATM) ....................... 300 Helium Mole Fraction Plot (1 ATM) ........................... 300 Molecular Weight Plot (Helium) ................................ 301 Compressibility Plot (Helium) ................................... 302 Compressibility Plot (Helium-Nitrogen mix) ................. 302 Plasma Density Plot (Helium) ................................... 303 Plasma Density Plot (Helium-Nitrogen mix) .................. 303 Electron Density Logarithmic Plot ............................. 305 Electron Density Logarithmic Plot (Helium-Nitrogen) ...... 305 Plasma Enthalpy Plot (Helium) ................................. 306 Plasma Enthalpy Plot (Helium-Nitrogen mix) ................. 306 Plasma Entropy Plot (Helium) ................................... 307 Plasma Entropy Plot (Helium-Nitrogen mix) ................. 307 Heat Capacity Plot (Helium) .................................... 308 Heat Capacity Plot (Helium-Nitrogen mix) ................... 308 xvi ”ii n ‘- - E18 H.19 H21 H22 H.23 H24 H25 H26 H27 H.28 H29 H30 H31 H32 H33 H34 H35 11.36 H.37 Electrical Conductivity Plot (Chebyschev Polynomial) 310 Electrical Conductivity Plot (Different Order Approx.) ..... 310 Electrical Conductivity Plot (0.1 ATM) ...................... 312 Electrical Conductivity Plot (1 ATM) ......................... 312 Non-Reacting Thermal Conductivity Plot (0.1 ATM) ........ 314 Non-Reacting Thermal Conductivity Plot (1 ATM) .......... 314 Reacting Thermal Conductivity Plot (0.1 ATM) .............. 315 Reacting Thermal Conductivity Plot (1 ATM) ................ 315 Viscosity Plot (0.1 ATM) ........................................ 317 Viscosity Plot (1 ATM) .......................................... 317 Diffusion Constant Plot (0.1 ATM) ............................. 318 , Diffusion Constant Plot (1 ATM) ............................... 318 Reaction Rate vs. e' Density Plot (0.4 ATM, 8000 Kelvin) .319 Reaction Rate vs. e‘ Density Plot (0.4 ATM, 10000 Kelvin) 319 Reaction Rate vs. e‘ Density Plot (0.4 ATM, 12000 Kelvin) 320 Reaction Rate vs. e' Density Plot (1 ATM, 10000 Kelvin) 320 Electron Density vs. Time Plot (0.4 ATM, 8000 Kelvin) ...322 Electron Density vs. Time Plot (0.4 ATM, 10000 Kelvin) ...322 Electron Density vs. Time Plot (0.4 ATM, 12000 Kelvin) ...323 Electron Density vs. Time Plot (1 ATM, 10000 Kelvin) .....323 xvii . Ebb-«mi Lt #313002 1 t J . . “ : - ' i'i ' ,. -- . . l ‘ A :‘.~ . me Molecule $11866sz Sound Speed We Field Impact Parameter WCap‘acity (at constanct pressure or volume) Diameter; Diffusion Constant DeSign For Six Sigma Design, Measure, Analyze, Improve, Control Design Of Experiment Defect Per Million Energy Electron Force Function Gravity Guage Repeatability & Reproducibility Enthalpy Planck’s Constant Hazard and Operability Study Electric Flux Bessel Function Riemann Invariants Siluare Root of Negative One Equilibrium Constant for Reaction "x" tam .1' --i an o v 1“! l/ MtConstam O ‘ Radius a ' mm ‘ ’1: . , Temperature SID-Stunts Time Vi ' Volume r Velocity X,‘ Mole Fraction of "i" Zr Compressibility Greek Characters a. Polarizability; Absorptivity 3“ M 31:: Propagation Constant 1‘ Deflection Angle ; 8 Energy Level; Emissivity T Y‘ Reduced Initial Velocity; Heat Capacity Ratio. 1?": Molar Flux. h“ MM‘Y ..'J kl‘ Thermal Conductivity; Scalar Value; Wavelength - ' 1:) DebyeLength. , Message Rotational Constant. i-'. v. ' i . ChemicalPotential;Mean it no 'uu o.“ 9:. am ..; ~..:‘Averase; Base. 1'2 Particles. mil“? .2 .‘r Effwfive.‘ Electronic " Cenu'ifugal. $351,? t CartesianCoordinates. til‘fl'ofv: 7.": ‘ c r32, Cylindrical Coordinates. 139.0, "‘ Spherical Coordinates. mm ...' ;' - TMModes. iti ith and jth Component. 33.1mm . ; Surface VP Constant Volume / Pressure “3 Rotational m "‘ “11'1" Translational & Himmler tyltwoml um fitting: r. t :1 g. . W’YJy-‘t; piggyizfit 3753“. . x.“ ;.' - _ on... . ' afulctfir‘way‘ fi . 7 - ".1- - ‘~‘_ ‘\ A} *- .tr. _-r.. -»‘\ twmmM.“ row . ‘Vr‘. L. " ‘ vb». ‘P‘F D ‘t -- , 1., . ID,“ 4 CHAPTER 1 INTRODUCTION 1.1 Problem Description. Plasmas provide a useful role in jet propulsion for space flights. They can be used for electric rocket thusters. One such thruster is the Microwave Electrotherrnal Thruster (MET). This type of electric thruster under research and development by NASA. Work was done at Michigan State University (MSU) to develop a better understanding of the transport mechanisms within the thruster. A combined joint effort of both the Electrical Engineering and Chemical Engineering Departments was done in this research effort at MSU. The research focus for the Electrical Engineering Department was on the development and optimization of the microwave cavity; whereas, the focus for the Chemical Engineering Department was on the transport mechanisms within the propellant and the cavity. A microwave plasma is very efficient for uses in jet propulsion. Production of these plasmas involves using plasma columns the size of conventional resonance cavities. These cavities are stable, reproducible, and quiescent. These plasmas develop as a result of surface wave propagation and are characterized by ion immobility. The major Variables involved with a microwave plasma system, such as one used for a Microwave Electrothermal Thruster (MET) system (described in more detail in Chapter 2) are: - I1 . ~- h.. "i. . “ o. v 0 Physical process conditions, such as the nature of the gas. 0 ‘Gas pressure. 0 Dimensions and material of the containment system. 0 Frequency and configuration of the electromagnetic field. 0 Power transferred to the plasma. In the electromagnetic environment, free electrons would be accelerated about the heavier and neutral molecules. These electrons collide with other elements or molecules of the gas to cause them to ionize as previously bound electrons are stripped off. In essence, kinetic energy would be transferred from the accelerated electrons to the gas. A cold propellant would receive microwave energy resulting in the production of a plasma. This plasma would give off radiation and heat energy. The excited species would flow away from the plasma while the could species would flow towards it. Downstream of the plasma, the excited propellant would be recombined with an increase kinetic energy. This thermalized propellant would exit through the nozzle as propulsion thrust. The research focus for my Master’s degree was on obtaining experimental data in the laboratory at MSU (Haraburda, 1990]. One set of data obtained in this research effort was for the energy distribution within the experimental microwave cavity. This energy distribution data included power absorption percentage of the propellant as a function of gaseous pressure and flow rate. Another set data was obtained for the dimensions of the plasma as a function of pressure and microwave power input. The last set of data was for . ‘ .-_, -7 2 Altar... the electron temperature of the plasma as a function of pressure. Additional data was also available from previous researchers in this area. The problem focused in this doctoral research was on the development of a theoretical model describing the transport mechanisms within the experimental microwave cavity. This theoretical model was based upon the empirical data received by other researchers and myself within the laboratory at Michigan State University. This model described the radiation energy transport within the cavity, the thermal energy and mass transport within each chamber of the cavity, and the reaction mechanisms within the discharge chamber of the cavity. The model simulations were for nitrogran and helium gaseous mixtures over a pressure range of between 400 - 1000 torr. 1.2 Problem Significance. These model equations are needed by NASA design engineers for building an operational Microwave Electrothermal Thruster. For example, it would not be desirable for a plasma to touch the walls of the discharge chamber because it would significantly increase the power lost from the propellant. These equations would accurately predict the plasma dimensions are various conditions. Also, some undesirable plasma reaction may result in losing the plasma. Thus, knowledge and prediction of the reactions within the plasma could be used for thruster design optimization and for proper selection of the propellant. ~~'J':" \. u..\-"‘ ’f3 '1 t ._ P n . a. t.’ 1 . I r I 1 l These equations could be used for non-propulsion applications. A chapter in this ' dissertation has been added to identify the following potential applications: 0 Detoxification of hazardous materials using the energy transport processes of plasmas. 0 Surface treatment of materials using plasmas to speed the growth of an oxide layer. 0 Novel methods such as constructing an operational microwave generated plasma chemical reactor. 1.3 Research Objectives. The following were the objectives of this doctoral research: 0 To describe the gases used within the experimental system. These were the simple monatomic (helium), diatomic (nitrogen), and mixtures thereof. 0 To describe the electromagnetic fields within the cavity system and to relate its effects upon the overall thruster system. 0 To develop the following transport modes with a goal of predicting the performance of an actual thruster. 1. Radiation heat transport model within the microwave cavity to explain the transport phenomena within the experiemental system. In addition to the nitrogen and helium mixtures, air at atmospheric pressure conditions was tie... n. 7’ v\ used for the outer chamber and the cooling chamber. The plasma, itself, l was modeled as a solid object. 2. Fluid transport model for the outer chamber, coolant chamber, and discharge chamber within the microwave cavity. The fluid used was nitrogen and helium mixtures at pressures near atmospheric. 3. Plasma transport model for the plasma reactions within the discharge chamber of the microwave cavity to explain the effects that condition changes have upon the reactions. 0 To identify the important system variables and their relationship with one another and within the thruster system. 0 To determine scale-up issues and concerns in the design and fabrication of a full- scale and operational MET. To accomplish these research objectives, the following were the research activities that were performed: 0 Conducted a literature review of parameters to be used in these model equations. 0 Developed a curve-fittin g procedure to use these parameters over a wide range of conditions (pressure and temperature). 0 Developed a statistical mechanics method for predicting the thermodynamic properties within the plasma. 0 Developed the computational simulations of the model equations. Conducted these simulations and verified the results with experimental data. fit. A ‘- ‘0 0 Used the NASA computer program to examine the nozzle performance applications using the information obtained from these model simulations. 0 Developed scale-up issues and concerns to design and fabricate a full-scale and operational MET. 0 Identified and described several non—propulsion applications which required the use of the same model equations and variables used to describe the thruster system. 1.4 Problem Solving Approach. First, the problem in this research was to define a steady-state analysis approach for each region within the system. Thies involved writing a macroscopic balance in each of these regions. Because of the complexity of the plasma region. a microscopic balance was also used. Finally, a relationship (such as correlation of the data) was done for each variable used in the model. Second, the basis of the model was developed. This involved the identification of the different types of transport processes (plasma, fluid flow, and radiation heat transport) that was predominant within the region . A computational technique was developed for solving the equations. This was done to conduct a microscopic analysis between the gases and the parameters in the model. \ '.. v.v o ‘Ils'n‘u ‘ i ..s .s... ‘K' " tantawasgathered This included the review of previous -- a» sahikoinvolvddlegatheringfimflmas ,, ,mmunufloruseintheequations. Thedaennodynamic ' ,4. parameters were obtained from various literature values (and . amount of-data available, selection of the appropriate data was “laminating of the database involved selecting the data (both empirical and ”Shawna obtained at similar conditions to that which was in the model We Typically, the most recent data was used. Fourth, the radiation heat transport model within the microwave cavity was done. The simple blackbody radiation model was not used. The realbody radiation model used thereflectivity of the cavity walls. During an experimental run for Haraburda’s Masters thesis, the condition of the microwave cavity wall was significant. A dull and dirty cavity wall reduced the efficiency of the energy transfer to the propellant. Thus, literature values were obtained for the emissivities of the various materials within the cavitiy. Plasma dimensional data and energy transfer data from Haraburda’s experimental research were used for boundary conditions in the simulations of the model equations. Fifth, the fluid transport model within each chamber of the microwave cavity was W Emeeonservation equations (continuity, momentum, and energy) were used to i i Wooden equations. The spatial dimensions and the boundary ' €95: . ,. . ‘7 - gym were used for the model aluminium. The Sixth, the plasma transport model, using a microscopic analysis, was developed. This model was used to characterized the reaction mechanisms within the discharge chamber. Statistical mechanics and curve-fitting methods were used for the necessary parameters. These parameters, which were primarily the thermodynamic properties of a plasma, were listed in Appendix H. Unlike the other model equations, this data was not obtained from the experimental system. Thus, the data obtained from the model similations would be considered predictions. Future experiments should be designed to verify these model equations. Seventh, the model equations using the computation methods were simulated and verified. The data from the model equation simulations were verified by comparing them to the experiemental data. The NASA program simulation was added to provide additional information on an application model. This applications model was used to describe the operating performance of an electrotherrnal rocket thruster. This was the interface between the experimental data and the applications. Finally, a discussion (literature search) identifying non-propulsion applications was done. Only a brief discussion was provided in this area as this step was not in the original problem scope. However, the variables developed in these models were still inlportant ones in plasma non-propulsion applications. .‘ i. s .4 v-~ ..ao...‘t v ~10... an». . u... .u‘ ‘u. ... ~¢~ . e 'u . i. . \ . . ‘ u 1-5 W The system used in this experiment was designed to conduct diagnostic measurements of three elements of plasma characteristics [Haraburda, 1990]. At the macroscopic level, the power distribution and plasma dimensions were detemrined using thermocouples and visual photography respectively. At the microscopic level, the electron temperatures were measured using an optical emission spectrometer. See Figure 1.1 for the overall set-up. 1.5.1 Microwave Cavity. An electromagnetic system was needed to generate a plasma. The microwave cavity body was made from a 17.8 cm inner diameter brass tube. As seen in Figure 1.2, the cavity contained a sliding short and a coupling probe (the two major mechanical moving parts of the cavity). The movement of this short allowed the cavity to have a length varying from 6 to 16 cm. The coupling probe acted as an antenna that transmitted the microwave power to the cavity. The sliding short and coupling probe were adjusted (or moved) to obtain the desired resonant mode. A resonant mode represents an eigenvalue of the solution to Maxwell’s equations. Two separate resonance modes were used in these experiments: TMon (L5 = 7.2 cm) and TMmz (L5 = 14.4 cm) [Chapman, 1986]. Additional features of this cavity included: two copper screen windows located at 90 degree angles from the coupling probe (which allowed photographic and spectral measurements), and two circular holes (in both the base and top plates) to allow propellant and cooling air flows through the cavity. G“ Gal Exhaust M Gal Pump v ' I - - - - ’ Thermocouple . : _ l-—----- Air Exh‘u't "‘ VIII]; 'IIIIl’ . ' 4 .. ..>. - _.9 _ _ _ -> Spectroscopy ” l " ' ; u I 7 Microwave ; 5 Cavity I , f Insulation ; I l I I Ill (I II 'I J h - - - a - -‘ ’ Thermocouple ‘ Water Outlet @G 66) Figure 1.1 Experimental Setup cob? LEGEND 1. Cavity Wall 7. Microwave Power 2. Sliding Short 8. Coupling Probe 3. Base Plate 9. Air Cooling Chamber 4. Plasma Discharge Fg. Gravity Force 5. Viewing Window Lp. Probe Length 6. Discharge Chamber Ls. Short Length Figure 1.2 Microwave Cavity 1.5.2 Plasma Containment. The plasma was generated in quartz tubes placed within the cavity (see Figure 1.3). The inner tube is 33 mm outer diameter and was used for the propellant flow. The outer tube was 50 mm outer diameter and was used for air cooling of the inner tube. Both tubes were about 2 V2 feet long and were epoxied to aluminum collars. These collars fed the gas and air to and from the cavity. For additional protection, water cooling was done on the collar downstream of the cavity. 1.5.3 Flow System. Flow of 99.99% pure nitrogen and helium was controlled using a back pressure regulator and a % inch valve in front of the vacuum pump. A Heise gauge with a range from 1—1600 torr was used to measure the pressure of the plasma chamber. Four sets of flow meters were used to measure the gas, water, and air flow. Thermocouples were used to measure the temperature of the air and water both entering and exiting the cavity. 1.5.4 Microwave Power. A Micro-N ow 420B1 (0—500 watt) microwave power oscillator was used to send up to 400 watts of power at a fixed frequency of 2.45 GHz to the cavity (see Figure 1.4). Although rated for 500 watts, energy was lost from the microwave cable, circulator, and bi-directional coaxial coupler. Connected to the microwave oscillator was a Ferrite 2620 circulator. This circulator provided at least 20 dB of isolation to each the incident and reflected power sensors. The circulator protected the magnetron in the oscillator from reflected signals and increased the accuracy of the power measurements. The reflected power was absorbed by the Terrnaline 8201 coaxial 12 Gas Outlet Water Inlet/Outlet Air Outlets Aluminum Input '1 I l Collar Air Cooling Passage _____ ____._ Thermocouple .__._.._ Quartz Tube (50 mm O.D.) Quartz Tube (33 mm O.D.) Plasma Gas Passage , , . , Aluminum Output Collar Gas Inlet Air Inlet Figure 1.3 Plasma Containment Tubes III—P 436A Power Meters Incl dent Rlflootod Power Power H-P 348“ Power Sensors Attenuators 20 dB Microwave Cavity Nerds. Microllne 3022 Bidirectional Coaxial ‘\\.\\\\\‘ Micro-Now 420131 --60 0 W) Microwave Power Oscillator Ferrite 2620 Tex-manna 8201 Ctrculator Coaxial Resistor (Dummy Load) Figure 1.4 Microwave Power Source ‘5' r" . 7‘ \ .w-u‘ ‘I“ "'~- v . 1 all \ s . 'I'a resistor. The incident and reflected powers were measured using Hewlett-Packard 8481A I power sensors and 435A power meters. 1.5.5 Temperature Probes. Type T thermocouples (copper constantan) with braided glass insulation were placed at the inlet and outlet for the water and air cooling (see Figure 1.1). An Omega 400B Digicator was used to measure the temperature at these four locations. 1.5.6 Spectroscopy. The radiation emitted by the plasma was measured using a McPherson Model 216.5 Half Meter Scanning Monochromator and photomultiplier detector. A high voltage of 900 volts was provided to the photomultiplier tube (PMT) using a Harrison (Hewlett-Packard) Model 6110A (DC) power supply. The output from the PMT was processed through a Keithly Model 616 digital electrometer. The processed output is sent to a Metrabyte data acquisition & control system and recorded on a Zenith 80286 personal computer (see Figure 1.5). The monochromator was positioned about 100 cm from the plasma. The emission radiation was focused on the monochromator using two 25 cm focal length glass lenses. This lens system concentrated the emission radiation on the entrance slit opening of the monochromator. To optimize intensity of the Spectroscopic emissions, the slit widths for this experiment were set at 100 microns for the entrance slit and 50 microns for the exit slit. The atomic spectra were taken using the 1200 grooves per millimeter grating (plate) with a range of 1050 - 10000 A. This groove E E lieu-ante 1‘01wa 516 DAS—ls Digital Electrometer Viewing Win dow Focusing Lenses H—P 6110A (DC) Power Supply Figure 1.5 Spectroscopy System mug allowed for a large range of wavelengths to be observed. The reciprocal linear Won was 16.6 A per millimeter. The focal length of the spectrometer was one half I' ‘1:ng 3.7 .1. kittiat‘.’ a ,_ .- . '. we > 1.6 Ex 'mental Results. The experimental results that were used for this research were obtained from Haraburda’s previous experimental research [Haraburda, 1990]. Measurements were taken for the helium and nitrogen plasmas for the following three parameters: 0 macroscopic power distributions 0 plasma dimensions 0 electron temperatures The experimental pressure range was from between 200 and 1000 torr. The gas flow rates varied from 0 to 2000 SCCM. The power input varied from 200 to 275 watts. The water cooling flow rate was 5.75 ml / sec. The air cooling flow rates were 2 SCFM for helium gas experiments and 3 SCFM for nitrogen gas experiments. Finally, the microwave resonance cavity was either the TM011 or the TM012 mode. The typical macroscopic power distribution was for the three flowing fluids: air, propellant, and the water. A simple calorimetry system was used to obtain the power distribution within the experiment. See Figure 1.6, which illustrates the power source, the radiation loss to the chamber wall and the fluid flows. The air was used to cool the quartz tube. The propellant was the helium or nitrogen gas. And, the water was used to cool the microwave cavity (brass chamber). Table 1.1 lists the approximate power distributions for these macroscopic power distribution experiments, which were used for the simulations in this doctoral research [Haraburda, 1990]. 17 PropellantT Air Figure 1.6 Calorimetry System 2’13 M ~ 4:; . £41,510 Wt . I .' v I—f‘ Tillie 1.1 Etperimental Power Distributions Helium Nitrogen Nitrogen The plasma dimensions were taken using a 35—mm camera, mounted on a tripod andpositioned about 2 centimeters from the microwave cavity wall (see Figure 1.2). Table 1.2 lists the plasma volumes for these plasmas using a 250 watt power source [Haraburda, 1990]. Two different measurements were obtained for the different regions of the plasma. The strong ionization regions were photographed as intense white color. The weak region were a different colored region, purple for helium and orange for nitrogen. The electron temperature was only done for the helium gas in a TM012 mode with no flow for the propellant and with a power supply of 220 watts. The spectroscopic F ‘ system was used to measure this temperature. Although the pressure varied from 400 to 3m tour, a small change in the temperature was seen. For this experiment, the electron Wattle was assumed to be that of the electronic temperature under the assumption »! , am eguilibrium occurred. As a result of this measmement and new Ii — HJ'I w nAhl-u. t‘ u... .. Xh-u. ,\ 1.2-u 0 Lu“. . \..- o 'l Strong Region (cubic cm) (cubic cm) 2.50 11.68 4.83 9.75 4.63 8.24 4.70 8.10 6.35 ' 13.15 5.10 10.29 4.85 8.53 4.77 8.01 m 1144 400 6.76 13.18 600 5.55 10.69 800 5.05 8.45 1000 4.77 8.37 Nitrogen 0 400 10.13 I 15.23 450 9.10 15.54 500 8.84 I 14.70 Nitrogen 102 400 11.74 16.39 450 9.91 15.83 500 10.14 16.34 Nitrogen 204 400 11.87 16.51 , 450 10.79 16.74 16.38 'M. 'V-\ . ‘N. V I ""w. As illustrated by the numbers within Figure 1.6, the experimental system can be broken into four sections: 1) radiation; 2) outer chamber; 3) coolant chamber; and 4) discharge chamber. The experimental results identified within this section are used for the simulations described in Chapter 4 of this dissertation. Figure 4.1 has a good schematic of the various regions that were modeled in this research. The conditions used for the simulations within Chapter 5 were listed in Table 3.1. Pure helium, pure nitrogen, and 25% (mole fraction) of nitrogen in helium were the propellant gases used in this research. Although several species are contained within Appendix A, only a few were used in the simulations because of the low temperature (less than 15,000 Kelvin). At higher temperatures, one would expect to see the other species (from electron ionization). As such, only the following species were used for these simulations: 21 MWM I;- . “1,.“ , n. \ ‘nw;1 “"55 . -r-l \~ s ‘ l . ' Wont, mm ranging from 250 to 4,000 ‘ , “the pressure was not held constant. In fact, the pressure varied from 0 l t', ms positions within the nozzle and discharge chamber. Additionally, ' em mixtures varied from pure helium to 90% mole fraction nitrogen. . ”Warrant 100% of the ions recombined within the discharge chamber (or Wm rttitused. In the simulations, many ions exited the system (which is another mm of power to the propellant). . As for the parameters used within the plasma, temperatures ranging from 300 to 50.000. Kelvin were used. The chemical kinetics and reaction rate simulations were done using collisional cross sections for ionization. Equilibrium conditions were used to calculate the associated values for the recombination reaction. ”fin-313%,? _ ”3% GS :‘nl. 2.2.1:; a. ,. reset-«1i terms. ranging from be! cumming. - a:- ' ’ . 'I -‘ 'r‘ 2 ' ‘ ‘ ‘. ‘iuwarswngi .,-, .v 5'3 3".“ 1* ._ A” .. a..", huh-u .5 '|> ... " V-l.l :~_.,~ - ..‘u, \ CHAPTER 2 BACKGROUND The use of plasmas in engineered products (such as rocket engines) is a relatively new area of technology. This chapter will outline its generic properties and general application. The ultimate goal of this research effort is for the use of plasmas in electrotherrnal propulsion. In addition to a brief discussion of this type of propulsion is the concept of using a microwave to form the plasma for use in this type of rocket thruster. Finally, this section ends with an outline discussion of the research effort involved in this dissertation. 2.1 Plasma Properties and Applications. Over a hundred years ago, a state of matter (other than the well-known solids, liquids, and gases) was observed. This state of matter was characterized by an enclosed electrically neutral collection of ions, electrons, atoms and molecules. Also, this state displayed relatively large intermolecular distances and large internal energy, resulting in a high degree of electrical conductivity. Because this substance did not display the characteristics of any of the three well-known states of matter, it was referred to as the "Fourth State of Matter," and later called the plasma state [CRC Press, 1988]. Plasmas may exist in several forms, ranging from hot classical plasmas found in “it magnetospheres of pulsars to the cold, dense degenerate quantum electron plasma of a 23 -.h white dwarf. Unlike the electrical insulation characteristics of a typical gas, plasmas . could be used as a useful electrical component, such as a good conductor. Because of these potential uses, plasmas have been artificially produced in the laboratory by such means as shock, spark discharge, nuclear reaction, chemical reaction (of large specific energy), and electromagnetic field bombardments [National Research Council, 1986]. Plasmas can have many applications, some of which include production of nuclear fuels, research and diagnostics in medicine, agriculture research, and environmental tracking of pollutants. Compared to conventional metal combinations, the use of plasma thermocouples can allow one to extract more thermoelectric power from nuclear reactors [Hellund, 1961]. As for military (and industrial) purposes, plasmas can be used for filtration systems in a toxic chemical environment [Carr, 1985]. Also, plasmas can provide useful sources for producing emission spectra from chemical analysis. Additionally, plasmas provide a unique and useful role in jet propulsion for space flight. 2.2 Electrothermal Propulsion. In general, there are three major types of rocket thrusters: chemical, nuclear, and electrical. Chemical rocket thrusters, in which energy is transferred to the working fluid through chemical reactions (combustion), are the most commonly used type of thruster, Nuclear rocket thrusters, in which energy is transferred through nuclear energy, are practically and politically difficult to use. Electrical rocket thrusters, in which energy is ”"3" 3‘ ‘ ir- transferred via heating coils or EM waves to the propellant fluid, are not practical in a -.h large force region, such as gravity from large celestial bodies [Dryden, 1964]. There are three basic types of electrical rocket thruster systems: I Electrothermal thrusters use electric energy to heat a conventional working fluid. I Electrostatic thrusters use ions or colloidal particles as the working fluid. - Electromagnetic (EM) thrusters use EM fields to accelerate the working fluid, usually in the plasma state. A proposed electrothermal propulsion system can use microwave induced plasmas. Although this system can use an electromagnetic wave, it would be classified as an electrothermal thruster because it would use a nozzle (not EM waves) to accelerate the propellant. Schematically shown in Figure 2.1 would be a version of this system [Hawley, 1989]. Microwave or millimeter power beamed to a spacecraft from an outside source (such as a space station or planetary base) could be focused onto a resonant cavity to sustain a plasma in the working fluid. The hot gas would expand through a nozzle to produce thrust. Alternatively in a self-contained situation, power from solar panels or nuclear reaction could be used to run a microwave frequency oscillator to sustain the PM [Haraburda, 1989]. Doomed lies-mu or Millimeter Fan Power Nozzle Energy Absorption lPropellunt Storage nun ' Velocity Gaseous Proponents Figure 2.1 Thruster 2.3 Microwave Induction. A microwave plasma could be very efficient for uses in jet propulsion. Production of these plasmas would involve using plasma columns the size of conventional resonance cavities. These cavities would be stable, reproducible, and quiescent. These plasmas would develop as a result of surface wave propagation and would be characterized by ion immobility. The major physical processes governing the (“Charge would be: (a) discharge conditions (such as the nature of the gas), (b) gas Milk) dimensions and material of the vessel, (d) frequency of the EM field, and (e) .5 ‘ ..... . 1' D I"! ,. Aroma} «L‘A- For microwave plasma electrothermal rocket thrusters, pressures near atmospheric and gas-temperatures near 2000 Kelvin were being investigated. In the EM environment, free electrons would be accelerated about the heavier and neutral molecules. These electrons would collide with other elements or molecules of the gas to cause them to ionize as previously bound electrons would be stripped off. In essence, kinetic energy would be transferred from the accelerated electrons to the gas. Figure 2.2 would illustrate the various discharge properties within a microwave system [Haraburda, 1990]. The cold prOpellant would receive microwave energy resulting in production of a plasma. This plasma would give off radiation and heat. The excited species would flow away from the plasma while the cold species would flow towards it. Finally, the plasma excited propellant would be recombined outside the plasma with increased kinetic energy. This thermalized propellant would exit through a nozzle as propulsion thrust. 2.4 Research Direction. The research in this dissertation has been separated into four areas and will be discussed in more detail in Chapter 3. 0 The first area would be the development of the model. The model involved the analysis of the separate heat transfer mechanisms (radiation, conduction, and convection of heat) and of the three major areas within the experimental system (outer chamber, coolant chamber, and discharge chamber). The second area involved the simulations of the model for the purpose of predicting and analyzing plasma behavior within a microwave discharge system. The third area was the development of the parameters needed for the model simulations. Finally, the last area was a literature review of potential non-propulsive applications using a microwave generated plasma. Thrust Microwave Figure 2.2 Discharge Properties 28 “’fl- l.'\ "I d..-»- r “tp~ +\\. . s..\. n. ..q CHAPTER 3 FLOW THROUGH A MICROWAVE GENERATED PLASMA Optimizing an electrothermal thruster system required a model characterizing the fluid within the thruster. Accurate models existed for characterizing the fluid both upstream and downstream of the plasma discharge region. The purpose of this dissertation was to link these two regions by developing a simple model to characterize the fluid within the plasma discharge region. As a motivation for this development, this model could be used in optimizing rocket propulsion such as using existing computer programs which NASA had for characterizing the fluid flow through a nozzle. This chapter outlines the basis behind the calculations of the model describing the fluid flow through a plasma. It also includes a description of the computational technique used in the calculation. Next, it will include a brief overview of the different types of parameters used in the calculations. Finally, a review of the past research effort conducted at Michigan State University in this field is done. 3.1 Basis of Model. This model development only considered using helium and nitrogen as a propellant. Unfortunately, these gases would not normally be considered as a pr0pellant (such as hydrogen and hydrazine) for rocket propulsion. However, the simplicity in using a monatomic gas allowed one to develop an easily understood theoretical model of fluid 29 ...i .Iv\'5" flow through a highly complex region (that of a plasma). The use of nitrogen in this model was used to demonstrate propellant contamination and the complexity of simulating & modeling polyatomic gases. This model was used for calculating the density, velocity, and temperature profiles of electrons, neutrals (i.e., N 2), and ions (i.e., He+, N2+). Table 3.1 listed the range of values used for the parameters within this research. Table 3.1 Research Parameters PARAMETERS VALUES E.M. Field TMO 12 Fluid Flow 0-1500 SCCM Power to Plasma ZOO-300 watts % Power to Cavity Wall 17% % Power to Air Coolant 38% Fluid Helium and Nitrogen Pressure 400-1000 torr Plasma Tube 33 mm O.D. Cavity Wall 17.8 cm I.D. 3O 3.2 Model Development Using Separate Transport Processes. Each of the separate regions within the microwave cavity system considered each of three transport processes : plasma transport processes, fluid flow transport processes, and radiation heat transport processes. As shown in Figure 2.2, this was a three dimensional problem. Plasma transport processes allowed one to develop individual particle actions within an electromagnetic and high temperature region. Statistical mechanics was used to augment this process in providing thermodynamic and transport coefficient parameters for the model equations. Fluid flow transport processes allowed one to develop particle and heat flow predictions using the conservation equations. These equations were linked to the plasma transport processes through the parameters calculated. However, these equations could have been linked rigorously through magnetohydrodynamic equations because the electrons would be highly dependent upon the electromagnetic field. However, because of their low relative mass and short residence time in a high velocity region, I neglected the electromagnetic effect on the ions in the plasma. 0 Radiation heat transport processes allowed one to develop a relationship and prediction for energy losses through electromagnetic means. Radiation losses accounted for the majority of the energy loss from the plasma propellant. This “a. r. x . La-I.I.. as»... ulna“: 1‘. . ll 1 x... ,1. process model was developed and simulated separately from the other two and linked together (after the simulations) as a one lump—sum value energy loss. 3.3 Development of Computational Technique. The model equations were highly non-linear and required sensitive methods for simulation with a plasma discharge region. Three types of numerical methods were considered for solving the transport equations, written as partial differential equations in computational space: Galerkin Method, Finite-Difference Method, and Finite-Element Method. These methods were considered because of their popularity in solving multi- dimensional partial differential equations. Error estimations concerning the numerical methods were analyzed. Computer programs were written to test the sensitivity (or error approximation) for each numerical method. Known solutions to various partial differential equations were used for this test. This included verification that each subroutine worked as intended. 3.4 Examination of Parameters. Three different sets of parameters were used for this model. The first set of parameters was the experimental research constraints, such as those listed in Table 3.1. The errors in these parameters were linked to experimental errors and not to the model development. The second set of parameters was the thermodynamic values used in the 32 at “‘9‘: I t v “his I Q.'C“~: 'h~.«~ in... .‘ V'.‘ «1" model equations. These parameters, used in the model simulations, were obtained through statistical mechanics (using empirically detemrined energy levels). The last set of parameters was the transport coefficients. These parameters were obtained through curve-fitting methods of previously obtained (known) data. These last two sets of parameters were compared to previously determined values and error ranges were discussed for each parameter. 3.5 Past Research Review. 3.5.1 Experimental. Besides my experimental work, R. Chapman conducted many experiments at MSU using the same microwave cavity system [Chapman, 1986]. His research focus was done on hydrogen gas at low pressures (0.5 - 10 torr) and at low microwave power (20 — 100 watts). Similar to my experimental work, he showed a 20% net power absorption to the exiting gas from the cavity system. Also, his calculated vibrational temperature range of 4000 - 17000 Kelvin was within the range of my gaseous (one-temperature) plasma system for helium. Similar to the results of my theoretical models, he measured an ionization percentage of between 0.001 and 0.1 % and demonstrated that the electron density increased with pressure and energy (temperature). Although my plasma system was based upon helium at higher pressures (near 760 torr) with higher microwave power (over 200 watts), the experimental data obtained by Chapman for hydrogen correlated quite well with that of mine. 33 V ‘ "J ‘ .. . {bet n...s_ ‘ a "' mm. 3.5.2 Theoretical. Two major research efforts at Michigan State University had been conducted by Chapman and Morin. 3.5.2.1 - Chapman. Besides his experimental work, Chapman provided a simple heat transfer model [Chapman, 1986]. He assumed a hard object (sphere) model with a plasma-wall interaction area. This model was useful for initial calculations of convective heat transfer. However, the actual phenomena of heat transport would be more complex. The heat transport within my work assumed this hard object only for radiative heat transfer and no wall boundary for convective heat transfer using a statistical mechanics based model of heating individual species (neutrals, ions, and electrons). 3.5.2.2 Morin. The majority of the theoretical work done at MSU concerning plasma systems was done by T. Morin [Morin, 1985]. Like that of Chapman, he focused his work solely upon a diatomic gas, primarily hydrogen. Thus, most of these equations used vibrational and translational degrees of freedom, which would not be present for a monatomic gas, such as helium. Morin’s main focus in his research work was that of modeling collision induced heating in a non equilibrium environment for a Weakly ionized plasma (less than 0.1% ionization). He combined the use of statistical mechanics and kinetic theory. Because of the equilibrium based nature of statistical mechanics, Morin spent much of his theoretical development on kinetic theory, which involved dynamic changes from one non equilibrium state to another. He used a Boltzmann based kinetic theory of gases. His simple chemical reaction models used both 34 95?: L‘ Ali-‘- 5 . , I” .Il;r b: bib ’,n 'r a..., r 5 i‘w‘j I.h\..g 1;- 5 .- ‘:\. ~\.s\,. N. ‘39 Plug Flow Reactor (PFR) and Continually Stirred Tank Reactor (CSTR) models. Some of the results of Morin’s research indicated that lower molecular weight molecules were superior as working fluids in collision induced heating, and that the concentration and temperature calculations in his PFR and CSTR models suggested a domination of the electron-molecule kinetics scheme. 3.5.2.3 Haraburda. The goal of this research work differed from Morin’s in that I took a simple equilibrium based theory to develop useful space- dependent parameters for popular transport equations. This research used the lowest molecular weight (as suggested by Morin) monatorrric gas (helium) for this first approach calculations. The kinetics involved in my work used a simple finite elemental section within the plasma using a Batch Reactor model to determine the residence time of the reaction to equilibrium. The kinetics calculation was only done to determine the difference of using equilibrium based and reaction based models in the plasma system. Thus, these kinetic calculations were not used for rate determining calculations. As such, the model equations within this research were for equilibrium conditions only. A detailed discussion of this difference was done in Appendix H3. 35 CHAPTER 4 MODEL DEVELOPMENT 4.1 Overview. Modeling of the experimental system was broken into four sections. The microwave cavity system was broken into three separate regions, which would be three of the four sections. The radiation heat transfer section (the remaining section) assumed a hard-sphere body for the plasma; whereas, the discharge chamber section did not make this assumption. Figure 4.1 depicts the four sections of this model as seen from inside the microwave discharge cavity. Each of these sections is described in more detail in this chapter with the assumptions and the resulting reduced equation (highlighted in black borders). From the reduced equation, the computational code was determined and provided. The parameters used within the code were determined from experimental data. These sections were not modeled (calculated) simultaneously. Instead, they were done sequentially using the results from one set of calculations from another section. The final model equation for each section would be identified by the equation enclosed within the bold outlined box. A comparison with experimental results would be provided with the results of the calculations. 36 - 33>) >...-...>.2..9 Cavity Wall / (1) / Real Body é Radiation . ca 5: ~ g / .53 / (2) (a) (4) ‘3‘ ; Outer Coolant Discharge é Chamber Chamber Chamber Figure 4.1 Microwave Plasma Model Overview The first section was that of the radiation heat transfer from the plasma to the cavity walls. The energy equation (F24), described in Appendix F, was used as the characteristic model for the realbody radiation within the cavity. This radiation was modeled through each of the three chambers of the microwave cavity. The second section was the outer chamber of the microwave cavity. The energy equation (E.l9) described in Appendix E was used as the characteristic model. Steady state conditions (B/Bt = 0) were assumed for this section. The third section was the coolant chamber. The energy equation (E.l9) described in Appendix E was used as the characteristic model. With the exception of the 37 coolant flow through the chamber, steady state conditions (ii/8t = 0) were assumed for this section. 0 The fourth section was the discharge chamber. The energy equation (13.19) and the momentum equation (E15), both described in Appendix E, were used as the characteristic model. With the exception of the propellant and energy flow through the chamber, steady state conditions (d/Bt = 0) were assumed for this section. 4.2 Realbody Radiation (Section 1). The body of the cavity was not considered a black body because it contained reflective metal (brass). Table 4.1 listed the emissivity for various conditions of brass and silver. The emissivity values in this table were average ones and not dependent upon temperature, wavelength, or direction [Siegel, 1981]. The condition of the cavity wall in my experiments was considered to be between dull and polished. It was my estimation to use the value of e = 0.2. Using this value in equation F .24, one could obtain a value for the plasma surface temperature. Figure 4.2 showed the plasma surface temperature as related to pressure for helium gas in the TM 012 mode (for the strong region). These calculations (using the equations developed in Annex F) assumed the radiation was emitted in a vacuum. This was clearly not the situation during the experiment. There existed at least five different media regions between the cavity wall and the plasma. As shown in Figure 4.3, numbered I through V with dimensions given, the media was identified in Table 4.1: 38 AXV IVN-‘FUNL-i—hhlih Table 4.1 Radiation Media Description REGION MEDIA DESCRIPTION I Air 1 atm, 300 K II Quartz 1.4mm thick, 300 K 111 Air 2 SCFM flow, 300-398 K IV Quartz 1.4mm thick, 300-700 K V Propellant Flowing, 300-1500 K (est.) Temperature (K) PLASMA SURFACE TEMPERATURE (TM 012 mode, Helium Gas) 1400 1 380— 1360- 134-0— 1320- 1 o . , , a - g - 503 o 4630 560 e60 760 300 9 0 Pressure (torr) Figure 4.2 Plasma Surface Temperature 39 nag"; .«r-mxrawhv .. . \..‘~ ..~ \ ..-., / _. / M as B / E s/ I II III IV V U, r; .2 t6 CL 0; fl 41:6.5 mnrg g] 25.0 mm \5 54 89.0 mm A; Figure 4.3 Radiation Environment Media These five regions absorbed, emitted, and scattered radiation. Thus, these regions magnified upon the complexities of the radiation heat transfer model. For example, fractions (fx) of the energy emitted by the plasma could be absorbed in each region (RP—9x )- This was graphically shown in Figure 4.4. Mathematically, this would look like the following: 40 hl‘ 2:1 - -~....ep >..~ ~>.~..v lap—)2 =f2 Ep Ep—>3 =f3 Ep Ep—>4 = f4 Ep Ep—95 =f5 Ep with the fractions adding up to l. Cavity Wall \\\\\\\\\:\\ fo=1 1 II III iv A\ T<1 /b A\ 1\ T 4\ o Figure 4.4 Radiation Emission Sketch 41 Plasma 4.1 ‘lu In. hub Each region could emit energy, thus producing six times as many relationships. A further analysis into this could be done by conducting an extensive study into determining the temperature gradient within each region. Nonetheless, a simple analysis was done by considering just the glass tubing media. An assumption was made that there was no effect by the media in regions I, III, and V. Regions II and IV were made of the same material (quartz glass). The emissivity and absorptivity (or) of this glass were not equal. The difference between the two would be the net amount of emitted energy, which was defined as: Cx :8x “1x 4'2 The following values for quartz glass were used [Touloukian, 1970]: aglass = 0.76 aglass = 0.03 A simple energy balance could be used (see Figure 4.5) resulting in the following energy relationship: and 42 - 4.. 2r \f um>lunv ‘t ,1 in» I II III IV V F4 '3 Cd 3 E3 E2 E1 5 >3 4 L L m +9 \ \ \ co "-4 F" g a. O \\\ 3 L\\\\\\\ Figure 4.5 Radiation Energy Balance through substitution of the above two equations, the following relationship could be seen: 4.5 E3 = C2 131 Therefore, the new energy balance at the surface of the cavity could be modified with the following being the model equation for the radiation section of this microwave plasma system: _ 4.6 4 'Qc =85(C2ApTg “AcTc )' 43 ;) mun-n - «bu r l . Ink “1 iv ‘1' "v x N“; 4.3 Outer Chamber (Section 2). The assumptions were made that there was angular symmetry, no net particle motion, and no more than 1 watt of heat conduction. This allowed the use of the steady state heat transfer equation. V2T=0 4.7 or (in cylindrical coordinates), 4.8 azr a r 3% + + = 0 a r2 r a r 8 22 Using the centered approximation Finite Difference method, the following approximations were made: it: ~ T(r + Ar,z)-T(r - Anz) 4-9 a r 2 Ar and, 4.10 32.1. .~. T(r + Ar,z)— 2T(r,z)+ T(r — Ar,z) a r2 (Ar)2 ,- I-‘l' 9’4 . u. vi “a . 4“ ‘ For ease in developing the computer algorithm, the following nomenclature was used: Tirl =T(r,z) 4.11 Ti+1,j=T(r+A r,z) 4.12 Ti,j+1=T(r,z+A z) 4.13 Now, the above heat transfer equation could be rewritten in a computer algorithm as: C1Ti+1,j +C2 Ti-1,j +C3 Ti,j+1+ C3 Ti,j-1-C4 Ti,j = 0 4-14 with the coefficients. being identified as: C1 = 1 + 1 4.15 (A r)2 2 r A r C2 1 1 4.16 =(A r)2 -21'A r 45 .1“... ‘.v|. 6 ‘u'rh‘. 5. ll‘ H.‘ e 4.17 C3 = 4.18 C4 = + _ For simplification in the computer algorithm, the two dimensional temperature was linearized as: Tid- = T(NZ[i -1]+ j) 4.19 With "NZ" being the number of nodes in the axial direction. The error estimation for this algorithm was second order. The error (E) was the difference between the actual temperature and the computed one. E = T(r,z)- Tm- 4.20 502(Ar,Az) The second order error estimation was calculated for the radial component usrng the following: 46 ru ‘1. h" I“; £ A1,3)63 4.21 02(Ar)s E v 25$§S89mm 6 3r3 s(Ar)3 Kr with "Kr" being a constant. The same estimation was done for the axial component. Thus, the overall error estimation was approximated to be: 02(Ar,Az)S (Ar)3 Kr + (Az)3 KZ 4°22 4.4 Coolant Chamber (Section 3 1. The assumptions were made that there was angular symmetry, no angular or radial motion, ideal fluid with a linear velocity profile, 125 watt heat transfer to the cooling air from discharge side wall, and negligible Viscosity effects. Additionally, heat capacity and thermal conductivity of the air were held constant. Several computer runs were conducted to check the model dependence upon changes in the heat capacity and thermal conductivity, which both changed with Changes in temperature. The energy equation for this region reduced to: 2 __ VT xi 32 and can be written in cylindrical coordinates as: 47 r I. r-l‘. ..3 l,_¥ « s p t-lh ul l. as.“ "tut ‘ 4.24 321‘ or azr pCp V7. in + + - —-—=O In the same method used previously, the computer algorithm could be written as: Cl Ti+1.j + Cz Ti-1,j + C31Ti,j+1 + C32 Ti,j-1- C4 Ti, j = 0 4.25 The coefficients, C1, C2, and C4 were the same as identified previously. The other two coefficients were identified as: C - l pCp VZ . 4.26 31 _ (A Z)2 2AA Z C32=(A Z)2 + 2AA Z Nevertheless, changes in densities and velocities of the air with position did not affect the results because of the continuity conservation law: .8 acovz) 42 DJ 3“ ..L s... tn; {-4.1 ‘mb .w .' "-". 9 f. r "a 1“\ 4' .1 bub—u“ The error involved within this algorithm was also second order. Thus, the error limit was the same as that modeled in the outer chamber. 4.5 Discharge Chamber (Section 4). The assumptions were made that there was angular symmetry, no angular or radial motion, ideal fluid with a linear velocity profile, constant pressure, steady state flow conditions, no viscous heating, and a 6.5 watt (net) heat transfer to the exiting fluid. The fluid simulated were helium and a helium-nitrogen mixture. Unlike the coolant chamber model, the viscosity, heat capacity, density, and thermal conductivity were not held constant. These transport coefficients were calculated using the statistical mechanics method discussed in Appendix D. This left two sets of unknown variables - the temperature and axial velocity. The energy (temperature) equation for this region was (in cylindrical coordinates): 4.29 or 3 rar azr C V ——-). + - rAH-+P- =0 P P Z0,Z [r3r[ 8r] 312] 21:1 1 in In the above equation, the Pin into the plasma was the average net power into the differential element. This was defined as the net power entering the microwave system subtracting out the power radiated to the cavity walls. The heat of reaction term was the ionization or recombination energy coming from the reactions within the differential Volume. The rate of reaction term, ri, was defined below, and derived from the continuity Cqmtion. 49 a 4.30 MWiaz( 21),) fl: With the substitution of the reaction rate, the energy equation would become: 4.31 or a rat azr AH, a cv—-x— -——v-p=o p P 232 [rel arl+azzl 2”,,an zPr)+ 1n Using the centered finite difference method, the computer algorithm for this could be expressed as: 4.32 Pi,jCPi,jVi,j i+1,j'Ti-1,j)_li. Ti.j+1+Ti,j-1-2Ti.j+ 2A 2 ,J (A r)2 Ti,j+1 -T1,j+1 +Ti+1,j +Ti—1,j 411,3] _ 211.jA r (A 2)2 V. . . . -V. . pl. . AH 1 pl , 1-1, 1—1, 2 l 1+1aJ 1+1] 1 .l .l +Plni j :0 lMW] 2A2 ’ The axial momentum (velocity) equation for this region was (in cylindrical coordinates): 50 .- twmexmce . 4 ' gratisxu 35- hasoeto amnhsBo thdumas r ".W‘ ‘h 1 handllml. 2 * 4.33 PVZaVZ- [ a{ravz)+a VZ]=O 82 TI rark 3r 322 Using the centered finite difference method, the computer algorithm for this could be expressed as: Vl+1,j'Vi-1,j Vi,j+l+vl,j-l'2V1,j 4.34 pi’j Virj -Tll,j + 2A 2 (A r)2 Vi.j+1 - Via-1 + Vi+1.j + V1-1, j - 2v... _ O 21.1,]. A r (A z)2 These two sets of equations could not be solved in the same way as that of the previous simulations. Both algorithms were non-linear and required an iterative solution using a method such as the Newton Method. The following variable sets were defined in this simulation: — T1 _ 4.35 3.4-.4. S\U¢~\ '30-an ~ULk‘V4w o n; the .1». 11 l» ftemp x) 4.36 Because these two vectors had two distinct regions, the Jacobian was broken into four regions, such as the temperature equation with respect to the temperature variables ('IT) and the temperature equation with respect to the velocity variables (TV). ’ g 7 4.37 T1" TV ' 4;) = . VT VV The Jacobian in the TI' region was calculated using the following: ”IT _ 2 AiJ 2 )‘iJ 4'38 1” _ (A z)2 +(4 r)2 1.7? =- lid. - A” 1,j+l (A {)2 21.1,]. A 1' 1T? -- Ai’j + li’j ”"1 (A r)2 ZIMAI rr _Pi.jCPi.jVi.j m J . .. 1+1.J 2A z (A z)2 JTT __ Pi. j Cm V1.1 - 31.1 i-l,j 2A Z (A Z)2 The Jacobian in the VV region was calculated using the following: vazpi,jlvi+l,j‘Vi-1,jl 21m 2111,]- . - + + 1’] 24 z (A r)2 (A 2)2 VV "Li “Li i’j+1—-(A 32-2434 1' VV =_ 1li,j _ 1'li,j l’J—l (A r)2 213.1 A r 53 4.39 4.40 4.41 4.42 4.43 4.44 4.45 It It: (‘1 vv.=Pi.jVi.i_ "Li 4'46 1+1.J 2A2 (A zjz va _ Pi.iVi.j "Li 4'47 i-l.j"' 2A 2 (A z)2 The Jacobian in the TV region was calculated using the following: TV Pi,j CPi,j 4-48 J 1,1 =Tfii+ijr i To1,j) TV _ AH] pli+1,j - 4.49 J. .-2 1+1.1 1 2ij A2 TV 1:2 AH] pli-1,j 4.50 1-1.1 , 2MW1Az Because there was no temperature in the velocity equation, the Jacobian in the VT region was set to zero. The error involved within this algorithm was second order for each iteration with each iteration converging quadratically. S4 WW 3 ('91.. .4 4 ,. V'uh. is W 0., f .31.. J16 a 2\l""'\ AM}. ‘ I fa." =M “*‘huL 3";“3 l‘. t.“ H u‘. .lJ33:",“‘H i3" , I Q.- 4.6 Validity of Assumptions. 4.6.1 Realbody Radiation. The assumptions made within this set of calculations appear to be accurate. Both the condition of the microwave cavity wall and the characteristic of the discharge tube material were accounted for in the calculations. As for the assumption that the air had no effect upon the radiation, the results of the calculations should not be too different had the emissivities of this air had been used in the calculations. The final assumption that should be accounted for was that of the hard body condition of the plasma. A more thorough set of calculations should not use this assumption, which is expected not to be too much different from that of these calculations. 4.6.2 Outer Chamber. The assumptions made within this chamber appear to be quite valid in that there was no particle motion within this section. However, tube wall boundary conditions were assumed using a linear and parabolic temperature profile with known inlet and outlet cooling air temperature. As shown through my simulation by changing this boundary condition temperature profile, the accuracy of these conditions was important in accurate calculations using the model equations. 4.6.3 Coolant Chamber. The assumptions made within this chamber also seamed to be quite valid. The gaseous flow should be close to having a linear velocity profile because this flow was turbulent (although barely) with a Reynold’s number of about 2100. Because it was flow through a linear tube, no angular or radial velocity 55 mum ' fine. if he 51 if? m. 31m 3PM - I \ni v ”Elm N \ - r'w‘hit,‘ ‘fi-u. N‘A A ._ “I ~ would have been expected. The temperature boundary conditions along the walls would have a significant effect upon the simulation calculations. The assumption that the heat capacity and thermal conductivity being constant would not be valid. These transport parameters would change with temperature. These changes would affect upon the results of the simulation. However, these changes would be expected to be smaller than that resulting from changes in the temperature boundary conditions. 4.6.4 Discharge Chamber. Unlike the previous two sections, many assumptions were made for calculations in this section. Unlike the coolant chamber, the gaseous flow should be close to having a parabolic velocity profile because this flow was laminar with a Reynold’s number of 2.8. However, as will be seen in the simulations, this boundary (velocity profile) would have an insignificant effect upon the results (i.e. temperature profile). Thus, for ease of calculation, the linear velocity profile was used. Assuming constant viscosity, heat capacity, density, and thermal conductivity was not done, as was done for the coolant chamber. However, these values were calculated through statistical mechanics assuming local thermodynamic equilibrium. From previous spectroscopic experiments, assumption of this equilibrium for atmospheric and low ionization helium plasmas appears to be quite valid [Dinkle, 1991]. Magnetohydrodynamic equations were not used because maximum ionization of the plasma was about 1%. Therefore, the majority of the species (about 99%) were neutral atoms and not directly affected by the electromagnetic fields within the plasma region. Like the previous two sets of simulations, the temperature boundary conditions would have an impact upon the simulations. 56 if": T10 if.“ “I; 8.45» Iii .'... .1: UFO 34"” N "L’ i O r"... ~. \1 ‘pl I E ‘ 4.6.5 NASA Program Simulations. The largest and most influential assumption made was not an assumption of condition, but one on dimensionality. These simulations were not done to provide an accurate portrayal of the system, but to provide an insight into the trends and magnitude of that portrayal. However, the assumptions made within this program module appear to be quite valid. 4.6.6 Summary. For simulations of the microwave cavity discharge system, the first and most important set of assumptions to be relaxed should be that of temperature boundary conditions. These boundary values could be determined experimentally by taking temperature measurements along the quartz tube wall. The velocity profile within both the coolant chamber and discharge section should be verified using trace measurements of visible or radioactive particles inserted into the fluid streams. As for the NASA program simulations, using the two-dimensional module would be the next generation of calculations required to accurately portray the operating performance of the thruster system. 4.7 Summary of Model Equations. The following table of equations contains those used to model the four sections within the microwave plasma system: 57 ll Really Out 00 g) / Table 4.2 Research Model Equations Section Equation Realbod Radiation y Qc =so(§2ApTg -Acré‘) Outer Chamber 321‘ a T 32T + + =0 a 1'2 I' a I' a 22 Coolant Chamber 82T 8T + 32'1“ _ PCp V2 31%) a 1’2 r a r a 22 2. a Z DischargeChamber C V 91-), 8 (1,31. +321. - p p z a z r a r\ 8 r a 22 ZAHii (V p-)+P =0 iMWiaz 21 in V 6V2- i rBVz +82Vz " Z az “rat a: 322 3-0 58 ‘ r “J. 1 . ‘- 5.3. :3. CHAPTER 5 MODEL SIMULATIONS 5.1 General. . The computer simulations for each section, along with the NASA computer program, were done on different computer systems. The real body radiation calculations were done using a calculator. Computational techniques were not needed. The calculations for the three other sections within the microwave cavity system were done using the enclosed FORTRAN programs on a VAX computer system (see Appendix A). These calculations could not be done on a personal computer because of the precision required in the calculations. The NASA computer calculations were done using their program on their VAX computer system network. 5.2 Realbody Radiation (Sectionfi. Using the more realistic equation, number 4.6, the plasma temperature was calculated and shown in Figure 5.1 [Touloukian, 1970]. Using the data in the Table 4.1, one could calculate the energy transported to the cavity wall for various materials and conditions. This assumed that the plasma surface temperature remained the same for each. Using the plasma surface temperature to be TS = 1500 K, a tabulation of energy transPorted was calculated and was shown in Table 5.1. As seen in this table, the surface material and conditions would have a large effect on the plasma. For example. I used a 59 l‘il ta will i ’5? ”1 why ”i“ FT.“ 1"! (Iv; 1K} '-_._-..-...m-<‘-1 r4:- 250 watt microwave source. If the threshold for sustaining a plasma was 200 watts, it would suggest that one could not maintain a plasma with a dull brass cavity wall. Thus, the material and condition of the discharge chamber would be important parameters for designing, maintaining, and operating electrothermal rocket thrusters. VRLASMA SURFACE TEMPERATURE (TM 012 mode, Helium Gas) 1600 1450- 1460- 1440— Tempem‘lure (K) 1420- 1 +00 . . . , . 300 460 560 not) 760 500 900 Presser re (torr) ——r 'I Figure 5.1 Radiation Plasma Surface Temperature liATERL- Gold Copper \ Silver Alumni Brass lunu Iron '\ QUIL- Table 5.1 Emissivity Values for Selected Materials MATERIAL CONDITION EMISSIVITY Absorbed Heat (watts) l====__l— Gold foil 0.009 2.25 plating 0.017 4.25 Copper polished plating 0.015 3.75 Silver plating 0.020 3.75 commercial roll 0.030 5 Aluminum foil 0.0294 7.35 Brass highly polished 0.030 7.5 polished 0.090 22.5 dull 0.22 55 oxidized 0.60 150 Yttrium film 0.35 87.5 Iron polished 0.078 19.5 Quartz fused crystal 0.760 190 Brick rough red 0.93 2325 Water smooth ice 0.97 242-5 61 A‘.» I‘- 5.3 Outer Chamber (Section 2). The domain of this region was the outer portion of the microwave cavity. As depicted in Figure 5.2, the radial length went from 25 — 89 mm, while the axial length went from 0 - 144 mm. Excluding the radiation heat transfer, the heat transfer through conduction was assumed to be 1 watt (see Figure 5.2). For this assumption, the following simple calculation was made: Q = h A AT = 1 watt Acavity = 0.10345 m2 Atube = 0.02262 m2 hair = 11.356W/m2K ATcavity z 1 0C ATtube z 4 0C 62 ZA. .......... 144 w R 0. \7 :§ .335 - l W 9 m 8 a.” .m m n 0.... w W... 5 2 £5 9:35 0 63 hint the abo ’0le be mad: 35313613011 6 )1 Printer Up an 166 a p21 immature \ 235.; "MW. [he Using the above temperature change calculations, the boundary conditions for the region could be made. Assuming a linear temperature gradient along the tube wall, the surface temperature was: T(25,z) = 301+95[1‘:7] 5-1 T(r,0) = 301 5.2 T(89,z) = 301 53 T(r,144)= 301 5,4 Two other sets of boundary conditions were provided to show the sensitivity of this parameter upon the numerical simulations. These additional boundary conditions assumed a parabolic temperature gradient instead of a linear one. The concave surface temperature was: 2 5.5 Whereas, the convex surface temperature was: hcse surfs gal mesh 5 .l schema 001131 2 5.6 T(25,z)=301+95[1-{144'z) J 144 These surface temperature gradients along the tube wall were shown in Figure 5.3. The grid mesh size and boundary conditions were varied to see their effects upon the results. A schematic of the grid mesh was shown in Figure 5.4. The following simulations were done: Table 5.2 Outer Chamber Simulations List GRID SIZE FIGURE T(57,72) 3x3 5.5 315.971 5x5 5.6 316.168 7x7 5.7 316.246 9x9 5.8 316.283 11x11 5.9 316.302 27x3 5.10 315.855 3x27 5.11 316.428 9x991 — concave 5.12 ‘ 309.520 9X9P2 — convex 5.13 323.046 65 Temperature (deg C) 4001—- 390; 360; 370'! 360: 350i 340i 330~ 320 310.1 300-4é Temperature (deg C) BOUNDARY TEMPERATURE (Outer Chamber Tube Wall) 4oo aeoj aso~ 370-: 380; 350' 3+0- 330 320; 310' zoo l U j l I 26' ' '4'6 '33 ' s'oj 160 réo' ' '140 Tube Length (mm) Figure 5.3 Boundary Temperature (Outer Chamber) Wham @QSLP WQWHOOU 0.0 mm 144 1 I. Cooling, Tube 3* :1, j; o 0 '23 Figure 5.4 Grid Mesh 67 Flg‘JI Temperature Gradient (outer chamber. 3x3 grld) 140- 1\ 120‘- i’ or 100- 3 LA BBQ 80- q 60- eng’th 310 ”E“ .5. Mo! L 40— 20- q O 25 1315125155 1615171585 Radial Length (mm Figure 5.5 Temperature Gradient (Outer Chamber, 3x3 grid) 68 £33 .1288 m... Temperature Gradient (outer chamber. 5x5 grid) 140- m M 100- ”E M -l {—15 50- E 1 E. eo- 40- q 20- 330 ngth .310 Axial Le D 25 '55 4151591315 1737155 Radial Length (mm) Figure 5.6 Temperature Gradient (Outer Chamber, 5x5 grid) 69 cameo; _e._x< res Flgu Temperature Gradient (outer chamber. 7x7 aria) 140 1204 i 100- l Len th immlg Axio 2555115155623 7955 Rodi i Length mm) Figure 5 .7 Temperature Gradient (Outer Chamber, 7x7 grid) 70 1 semen: 6.34. Figure 5. Temperature Gradient (outer chamber. 9x9 grid) O 25 '3'5 '4'5 ' 5'5 ' 6'5 251—5's Radiai Length (mm) Figure 5.8 Temperature Gradient (Outer Chamber, 9x9 grid) Temperature Gradient (outer chamber. 11x11 grid) 140 El 120- 100- :5 g)" 501 3 E - ,e 6— 60- .? . 40— 201 O'T—‘i'ltl'l'fj 25 :55 45 55 65 75 65 Radial Length (mm) Figure 5.9 Temperature Gradient (Outer Chamber, 11x11 grid) O Ifigure 5.1 Temperature Gradient (outer chamber, 27x3 grid) 140- 100-> m 50... E -l w 60" 1 4o- 340 330 320 ength 5‘10 Axial L 20- o , . ,_ . . , 25 3'57'5 5'5 6'5 7‘5 as Radiéalm L ngth Figure 5.10 Temperature Gradient (Outer Chamber, 27x3 grid) Figure 5. Temperature Gradient (outer chamber, 3x27 grid) 146- 120- a 100- E e 3 ’3 n - R D 3.35 5°_ 9’ :9 ~55. t... 613--1 Q 40- 20- O J 293315 3'5'6'5'5'5'6'5 Radial Length (mm) Figure 5.11 Temperature Gradient (Outer Chamber, 3x27 grid) Figure 5.1 Temperature Gradient (outer chamber. 9x9p1 grid) 140 \\\ racy-a 6 (WE 100 50 60 Axial Length (mm) 40 20 O [FT—'T'IT—IRT'i 25 235 4-5 55 85 75 85 Radi 1 Length mm) Figure 5.12 Temperature Gradient (Outer Chamber, 9x9p1 grid) “we: Temperature Gradient (outer chamber. 9x9p2 grid) 120 1 L111 100 " .3. gth E; . 50 3K] 3H) ’5? 5.. 60 Axial Len 4O 20 O n * IT—FT I""i 7—1 F r 25 35 4-5 55 65 75 85 Radi 1 Length mm Figure 5.13 Temperature Gradient (Outer Chamber, 9x9p2 grid) 76 156,131th these Si high temperature Tet" berm much finer, shotn in Figure 55 tithe graph. This PC 5.311159. the graph on another. suggesti larger than 919. As 0 indirection affecte radial direction. This it axial direction. Fi boundin' conditions :31 mesh sizes were i obscned at the point ( jfptndent upon the [U traction. As expected, these simulations predicted an ellipsoidal temperature gradient around the high temperature region of the tube boundary temperature profile. As the grid size became much finer, the temperature profile converged to a more predictable solution. As shown in Figure 5.5, a small temperature peak appeared near the upper right hand comer of the graph. This peak disappeared for finer grid meshes. As compared between Figures 5.8 and 5 .9, the graphical illustration of the temperature profile appeared to be identical to one another, suggesting that the grid size did not alter the calculations for a mesh size larger than 9x9. As compared between Figures 5.10 and 5.11, the grid step size in the axial direction affected the simulation calculation much more than the step size in the radial direction. This would be expected based upon the higher temperature gradient in the axial direction. Finally, Figures 5.12 and 5.13 illustrated the effects that temperature boundary conditions would have upon the simulations. The errors resulting from these grid mesh sizes were illustrated in Table 5.2 listing the predicted temperature data observed at the point (r = 57mm, z = 72mm). Thus, the algorithm error would be highly dependent upon the tube wall boundary condition and upon the step size in the axial direction. 5.4 Coolant Chamber (Section 3). The domain of this region was the cooling chamber between the quartz tubes. As depicted in Figure 4.14, the radial length of the region went from 16.5 - 23.6 mm, while the axial length went from 0 - 144 mm. Neglecting the radiation heat transfer, the heat conduction was assumed to be 125 watts. The velocity of the gas into the region was '7‘) 1.0551n/soc at a tem lg it the exiting tem constant with an ave temperature change 1 to about 375 C. Usi muldbe made. Asf linear. parabolic con rods that the maxim Lhiee101111115 the way than 1.055 m/sec at a temperature of 300 K. Based upon the heat capacity of air at 1.0467 kJ / kg K, the exiting temperature was 398 K. The thermal conductivity was assumed to be constant with an average value of 3.00 W l m K. Similar to the calculations of temperature change for the outer chamber, the temperature change for the inner wall came to about 375 C. Using that temperature change, the boundary conditions for the region could be made. As for the outer region, three sets of boundary conditions were used - linear, parabolic concave, and parabolic convex tube wall gradients. (An assumption was made that the maximum wall temperature was near that of the plasma - approximately three fourths the way down the tube.) Linear: T(r,0) = 300 5.7 T(r,l44) = 398 5.8 T(23.6,z)=300+98 3— 5.9 144 T(16-5,Z)=300+375[1—(7;—8] v oszslos 5-10 =398+277[144'z] v108szs144 78 144 =i/\ Inner Tube Wall Tube 521 A' 3 111' 5 O O .......................... >R 16f1 2313 Figure 5.14 Temperature Gradient (Coolant Chamber Sketch) Contrifl TP1(2316 Cortex: TD: (23.6 1p; (16.5. in surface temperaill 111111311003 were done Concave: z 2 5.11 Tp1(23.6,z)= 300 + 98(m] T (165 ) 300 3 z 2 5'12 .,z= +75—— vos $108 pl [108) Z 2 =398+277[I4:6‘Z) V108525144 Convex: 2 5.13 144-z T 23.6,z =300+981- p2 ) x [[144 j] 2 5.14 Tp2(16.5,z)=3OO+375[1-[l?:éz] ] v 0323108 z-108 2 =398+2771-[ 36 J V10852$144 The surface temperature gradient along the inner wall was shown in Figure 5.15. Again, simulations were done varying the grid mesh and boundary conditions. 80 Temperature (deg C) 'zo' '. 40' ' 'a'o' ' '8'0' "130' ' Bid ' '140 Tube Length (mm) Figure 5.15 Boundary Temperature (Coolant Chamber Table 5.3 Coolant Chamber Simulation List GRID SIZE FIGURE 3x3 5.16 9x9 5.17 9x9pl - concave 5.18 9x9p2 - convex 5.19 81 Figure 5. Temperature Gradient (coolant chamber. 3x3 grid) 140- o-§\\ 1 00 .C 1 465A 50 . 60- .15? 0 j I U Ivl U ' ll "J 16.5 16.5 23.5 22:5 Radial Length (mm) Figure 5.16 Temperature Gradient (Coolant Chamber, 3x3 grid) 82 Hgms Temperature Gradient (coolant chamber. 9x9 grld) 140- 120- 3% mo../ 50- 60- 4o— V a or 500 20 Pedal Length (mm) Q 43% zon/gmo 01,..rfim ...pa 16.5 15.5 26.5 22.5 Radial Length (mm) Figure 5.17 Temperature Gradient (Coolant Chamber, 9x9 grid) I:l'E-Ture 5.1 Temperature Gradient (coolant chamber, 9x9p1 aria) 140- J 120- t} 0 0 100'- KW .C 44 ‘i gm 50— B E _ “5‘70 —.._e ‘53 60-1 ‘5 >9 * e 404/ 20—1 0 I'" Tjfil'ij I T 1‘7 "T I 18.5 18.5 20.5 22.5 Radial Length (WW) Figure 5.18 Temperature Gradient (Coolant Chamber, 9x9pl grid) 84 Flgul‘e 51 Temperature Gradient (coolant chamber; 9x9p2 arid) ill-C)d 12mm “:tha c: (:3 ~39 § '08 J 80- | Length ,_, E .5 eo-J Axia 40W C) 2 O _ dbfi/ -l 315:0 O fjlj—FT'Tj—fi'rfi—‘I 16.5 18.5 20.5 225 Radial Length (mm) Figure 5.19 Temperature Gradient (Coolant Chamber, 9x9p2 grid) 85 like the simulation 1 ellipsoidal shape Wt! drape. This profile i tube with the temper 5.17. the temperature lordte 3x3 mesh grit cloth in Figures 5.1 condition had on the laundry condition p audition would r10t The domain 0 tom Figure 5.14, the l: omOto 144 mm. 1‘ liar- ' a. lhe urscositv : *l . he i v th order tempera marry, and electron to»: " ' noed m the first 31;“ i etube was assu interline m assum vi, «at ' Pehmental calc Like the simulation for the outer chamber, the temperature profile appeared in an ellipsoidal shape with a slight distortion towards the downstream end of the geometric shape. This profile was expected for a non-reacting coolant gas of a linear cylindrical tube with the temperature boundary conditions provided. Comparing Figures 5.16 and 5.17, the temperature profile for the 9x9 mesh had a lower temperature gradient than that for the 3x3 mesh grid. As expected, results of the change in the temperature profile shown in Figures 5.18 and 5.19 demonstrated the large effect the temperature boundary condition had on the simulation. As shown in Figure 4.19, the convex parabolic boundary condition produced unrealistic results, suggesting that this type of boundary condition would not exist. 5.5 Discharge Chamber (Section 4L The domain of this region was the discharge within the quartz tube. As implied from Figure 5.14, the radial length went from 0 to 16.5 mm, while the axial length went from O to 144 mm. Neglecting radiation heat transfer, the net heat absorption was 6.5 watts. The viscosity and the thermal conductivity for each point were calculated using the 5th order temperature polynomial calculated in Appendix H. The density, heat capacity, and electron density were calculated using the Statistical Mechanics subroutine described in the first part of Appendix H. The temperature boundary condition along the quartz tube was assumed to be linear (similar to the cooling chamber) and down the centerline was assumed to be sinusoidal with the maximum temperature selected equal to the experimental calculation of electron temperature in my master’s thesis research. The 86 inlet temperature w [based upon the net following were the t The “106in boun distant mass flow inlet temperature was assumed to be 300 K and the outlet temperature to be 1100 K (based upon the not heat transferred). For the simulation with a pressure of 400 torr, the following were the temperature boundaries: T(r,0)= 300 5.15 5.16 T(0,z)== 300+800[%) v OS 2 s72 =7550+6450$in[7t-(—z-é—69—0—-)] V 72323144 5.17 T(l6.5,z)= 300+ 1 root—Ed v 0: z 5108 =1100+300[144'Z] v 108323144 T(r,144)=1100 5.18 The velocity boundary condition was based upon the ideal gas law equation with a constant mass flow: 0 T 5.19 Flow = n R 87 [sing the dimensio pressure of 400 I0”, Three sets of figures 5.20 to 5.22 expected. this tempe of 'he plasma Asp fro: the plasma Th ru * ' ' ll axed 1n FigureS p; Using the dimensions of the tube, the volumetric flow rate of 572 SCCM, and a constant pressure of 400 torr, the velocity at the boundary conditions could be approximated to be: V=0.004233T (units are: r%mn) 520 Three sets of simulations were done, at pressures of 400, 600, and 800 torr. Figures 5.20 to 5.22 showed the temperature gradient in Kelvin for pure helium. As expected, this temperature profile appeared in an ellipsoidal shape, the same optical shape of the plasma. As pressure increased, the temperature gradient decreased more away from the plasma. The temperature plot of a 25% mole mixture of nitrogen in helium was portrayed in Figure 5.23 for 600 torr pressure. Not much difference was seen between pure and mixture components for temperature. Figures 5.24 to 5.26 showed the velocity gradient in m/min. Like that for the temperature, the velocity profile appeared ellipsoidal about the plasma. The fluid velocity in the plasma increased to about five times that outside it. Thus, it would be unexpected to see the neutral fluid flow around the plasma. Instead, it should all go through the plasma. As the pressure increased, the velocity decreased proportionally with the gas law relationship for pressure and volume. The 25% mole mixture of nitrogen in helium was portrayed in Figure 5.27. Like that of the temperature plot, not much changed in the velocity between mixture and pure components. 88 £964 .25... Temperature Gradient (discharge. 400 torr) 14o \ 9 120 9a K 0 8 100 .1: // a... 3 E 50 —.§. .5 150 .h. 0 Pa O O ..T d i e 1'2 1'5 Radial Length {mm} Figure 5.20 Temperature Gradient (Discharge Chamber, 400 torr) 89 cement. Bees, 5 TC Figu Temperature Gradient (discharge. 500 tart”) .: 'Eia re 3 E 5““ 150 4a 20 O 1"'1'r'r'flrT‘fi—r a 4 a 12 re Radial Length (horn) Figure 5.21 Temperature Gradient (Discharge Chamber, 600 torr) — - : FTgUfe Temperature Gradient (discharge. 800 tarr) 140- 2000 4%480 \ 120-%~ 0 ‘éga 1 5’: ..E E d .I/ l ...J ~52 g. 50- 4a- 20- 04'r_'l'r1—l"‘—l"1 0 4- 8 12 16 Radial Length mm) Figure 5 .22 Temperature Gradient (Discharge Chamber, 800 torr) 91 3 1.. .3. .91». F .21.:qu— .DmXxu. Temperature Gradient (discharge, 600 tarr. mixture) o ...,...,e..,... o 4 a 1216 Radial L n th tmmi 9 Figure 5 .23 Temperature Gradient (Discharge Chamber, 600 torr, mix) 92 facet. .a..x< Flgu Axial Length (mm) Velocity Gradient (discharge. 400 tarr) ‘m .. e . x8 Bow—~10“ 50' l 40- 20.1 ’ 1 O .m . . . . . . . . a 4 g 1'2 1‘6 Radial Length mm) Figure 5.24 Velocity Gradient (Discharge Chamber, 400 torr) 93 emcee 6.8.x £ Figu Velocity Gradient (discharge. 600 tarr) 14a- 5 'ngis‘ 1213-335? d> $3 I \9 ., 5 Axial Length (mm) O "’l"’l"'l"“7 O 4 8 '12 16 Radial Length (moo) Figure 5 .25 Velocity Gradient (Discharge Chamber, 600 torr) Axial Figure 5 Velocity Gradient (discharge. BOO torr) — fi 140 6 120- lbrial Length {mm} o "'T"'l"'l"'l a 4 e 1216 Radial Length (than) Figure 5.26 Temperature Gradient (Discharge Chamber, 800 torr) 95 Lvmcwl— .Ddes Figure . Velocity Gradient (discharge, 600 torr. mixture) —.¥ 100? BO- Axiul Length (mm) o m” 01TE1'216 Radial Length (MM) Figure 5.27 Temperature Gradient (Discharge Chamber, 600 torr, mix) 96 PlasmadlSCharge' sccnon. W3 25% 532. Unlike the I different between hat it would be e contamination. Th ntngen into the h min of 73. the dat Nonetheless. the If Figures 5.28 to 5.31 showed the electron density gradient in #ICC. Because of the large electron density gradient in the plasma, graphs of this were done using the plasma region, as opposed to the entire flow region. This profile showed an ellipsoidal shape and could directly be linked to experimental Optical (photographic) measurements taken of the plasma discharge. This had been compared to the experimental values in the following section. The 25% mole fraction mixture of nitrogen in helium was portrayed in Figure 5.32. Unlike the temperature and velocity, the electron density gradient was remarkably different between the mixture and pure component. This difference came from the fact that it would be easier to strip electrons from nitrogen than it would be from helium. 5.6 NASA Proggar_n Simulations. Using the One-Dimensional Equilibrium computer program, one could determine the effects on engine performance from pressure and energy changes along with propellant contamination. The propellant contamination effects were inferred from inserting nitrogen into the helium gas. Although the nozzle geometry allowed for an expansion ratio of 75, the data provided in the following simulations only went to a ratio of 10. Nonetheless, the trends were demonstrated. 97 |3|411IIQ|14| :3th. 3.3.». ,. FlSure 5 Electron Density Gradient (discharge. 400 torr) I EH4 140 120 100- 4: / I'd en. t l 3,5 50- ...g l .5? 60- 2.5 i, 40— l 20- 0 . O Hill 4 I—IIVI— a '1'2H16 Radial Length (MM) Figure 5.28 Electron Density Gradient (Discharge Chamber, 400 torr) 98 Electron Density Gradient (discharge. 400 torr) 130 125- 120-f9é‘ o _a 0' 0 1 I \3' ' N. 100 - 4929* . )4 95-1/ ”(v 90- Axiol' Length {mm} .\ 85 .— ...... l oliéiééfl. Radlal Length mm) Figure 5.29 Electron Density Gradient (Discharge Chamber, 400 torr) 99 Electron Density Gradient (discharge. 500 torr) 130 n h tmmig1 Axial Le Radial Length (mm Figure 5.30 Electron Density Gradient (Discharge Chamber, 600 torr) 100 Electron Density Gradient (discharge. 800 torr) l 130 120- % 115- A u “E :5 110- D‘s-M 3;“ 105- too-i 95- 90 ... m a 1[ {'3 5. 5 6 Radial Length (MM) Figure 5.31 Electron Density Gradient (Discharge Chamber, 800 torr) 101 Electron Density Gradient (discharge. 600 torr. mixture) 130 l 1 133+ 1254 1“ 120-3\\ ‘l 115-i “’5’“ Z 9I+ 9I+IEI SI+3 S 1104 fiscal Length (MM) is" I \ 9I+3 __l_ 1001 95— x9. 4 5» x 901/ 35.. .. . .1 . oiiflflflgro Radial Length (MM) Figure 5.32 Electron Density Gradient (Discharge Chamber, 600 torr, mix) 102 Specific Impulse Pressure Plot (Throat Ratio - 10. 4 Watt) 1700 1°B°i A a; . ‘3“; 1600- ”‘E no :2 i a” 1550- 150c . . , . , . (mo cfis 'Lo L5 :20 Pressure (Atmosphere) Figure 5.33 Specific Impulse Pressure Plot 5.6.1 Pressure Changes. The pressure changed from 0.1 - 2.0 ATM in a 4 - kWatt thruster. The specific impulse was plotted against these pressures at the position in Which the nozzle expansion ratio was 10. As seen in Figure 5.33, the specific impulse increased with pressure with a large increase for pressures less than 0.5 ATM. 5.6.2 Energy Changes. Five different energy regions, ranging from 250 - 4000 Watts, were inferred from the mole fraction and temperature data. The figures had data Plotted against were a ratio to throat. The throat was defined as l on the graph within the discharge chamber. The helium mole fraction gradient was plotted for all five energy regions. As seen in Figure 5.34, power less than 500 watts produced almost negligible ionization. Additionally, the temperature gradient plot was provided in Figure 5.35. These suggested that the instability of maintaining a plasma existed at powers less than 103 500 watts. For all five energy regions, the pressure gradient along the nozzle axis remained the same. Figure 5.36 showed that plot. As expected, the large pressure drop occurred at the throat. Likewise, the mach number gradient remained the same for each energy region. This gradient plot was provided in Figure 5.37. Finally, an important measure of engine performance was the specific impulse. Figure 5.38 had this plot for each energy region. For specific impulses greater than 1000, the simulations suggested that one needed to use at least 1 kWatt power. Helium Mole Fraction Gradient (5 Difforont Energy Region.) A A A A _ A 1'0 V Z J 260 Watt eoowmt 1 Watt c ZWOt-t .9 4kWat-h ‘3 D a: $3 .. o 2 07-i .5 .. CHAMBER EXIT Tr : 0.5- 0'5e'iri'i'k'é e 10 Area Ratio to Throat Figure 5.34 Helium Mole Fraction Gradient (5 Energy Levels) 104 Temperature Nozzle Gradient (S Differ-ant Energy Region.) 2.5E+ 4 Watt rm 0 Fe? 2'05.” can Watt 25 25a Watt :2 v 1.5E 5: i‘ '— 5000.0 Area Ratio to Throat Figure 5.35 Temperature Nozzle Geometry Gradient PreSsu re Nozzle Gradient (All 5 Different Energy Regions) 19°“ 0‘8~l A a 8g ' CHAMBER ” EXIT E 0.4- “ 4‘5.- _ 0‘2” 000 a . 4!- u i i fl # r 5 ' 1o Area Ratio to Throat Figure 5 .36 Pressure Nozzle Gradient 105 nic velocity) Mach Number (ratio to local so élbf sec / lbm) peclfic mpulse Mach Number Nozzle Gradient (All 5 Difforont Energy Region.) EXIT 3- CHAMBER _ (supersonic) (subsonic) i'i'e'é'é'ro Area Ratio to Throat 8 4 Figure 5.37 Mach Number Nozzle Gradient Specific Impulse Nozzle Gradient (5 Different Energy Regions) 2000 4- Watt 2 Watt 1 kWo-tt 600 Watt . 250 Watt 10 Area Ratio to Throat Figure 5.38 Specific Impulse Nozzle Gradient 106 5.6.3 Nitrogen Mixtures. To foresee the effects of propellant contamination, several simulations were performed using helium-nitrogen mixtures. Using the 4 kWatt power, the mass percentage of nitrogen was varied from 0-90%. The helium mole fraction plot was given in Figure 5.39. It was seen that the helium molecule was negligible upstream of the throat for a 90% mass nitrogen mixture. The other important species was the electron. As the nitrogen mass percentage increased, so did the mole fraction of electrons. For a 90% mass nitrogen mixture, half of the species was electrons. As seen in Figure 5.40, this was twice as many electrons than for pure helium. Again the specific impulse was important. Seen in Figure 5.41 was a plot of this versus the nitrogen mass percentage. It reached a minimum value around 50% and increased dramatically beyond that value. Although the specific impulse may have increased for high nitrogen mass percentages, the mach number decreased as illustrated in Figure 5.42. Helium Mole Fraction Gradient (5 Different Nitrog-n Mixtures) 1.0 H oxNz " H 13 N2 '1 H 10% N2 = 0'3“ H (50:: N2 _ g j 0—0 90!: N2 , O . Lt- 0.6- % - _ 2 04 ;—fi fi _ “S. : CHAMBER EXlT E ‘ I 0'2— k 0.0 o 3 . 3* i ' i ' é A ‘10 Area Ratlo to Throat Figure 5.39 Helium Mole Fraction Gradient (5 Mixture Levels) 107 Electron Mole Fraction Gradient (5 Different Nitrogon Mixture.) 095-¥ f C t z 0 4 i .9 ' . E5 . m L ‘ A _ ‘5 1g ; A .2 z r ‘ ' ‘3 g 0.2- EXlT “8' . H cox N2 L H 108 N2 ' H 1% N2 0.0 “7 °’,‘ N2. . . . . a 4 i i 1:» é é 1 0 Area Ratio to Throat Figure 5.40 Electron Mole Fraction Gradient Specific Impulse Mixture Plot (Throat Ratio = 10. 4— kWatt) ec lbm 'fic lr/npulsg 8 l 522; o ' 2'0 ' 4T0 ' a'o ' 8'0 ' 100 Nitrogen Mass Percentage Figure 5.41 Specific Impulse Mixture Plot 108 Mach Number Mixture Plot (Throat Ratio = 10. 4- kWatt) 5.05 5000 - 4.95 - 4090 '- Moch Number 4.85 - 4'80- 4075 I [ l I u l y I . O 20 4—0 60 80 1 00 Nltrogen Mass Percentage Figure 5.42 Mach Number Mixture Plot 5.7 Comparison with Experimental Results. Except for the width of the plasma, the data depicted in the graphs corresponded quite well with previous experiments, those of Whitehair and of Haraburda. The difference in the width of the plasma could be corrected by taking a more accurate value for the electron density in the simulations for the measurement. The values used were for an electron density of lxlOl4 electrons per cubic centimeter. Table 5 .4 listed the Volumetric measurement comparison between the simulation and experimental data. The data from this experimental research were for the strong discharge region. The data from Whitehair’s dissertation research were for no-flow in 37mm inner diameter tubes [Haraburda, 1990 and Whitehair, 1986]. 109 Table 5.4 Simulation vs. Experimental Values. Pressure Plasma Power Width Length Volume Power Density (ton) (W) (W/CC) (mm (mm (cc) 400 165 45.83 13.8 36 3.60 Simulations 600 167 60.95 12.4 34 2.74 (theory) 800 169 78.60 11.0 34 2.15 400 165 - 18 36 5.35 Haraburda 600 167 - 17 34 5.10 (experiments) 800 169 - 16 34 4.85 g 474 441 79.28 - - 5.64 Whitehair 584 456 97.99 - - 4.72 (experiments) 760 469 128.7 — - 3.71 110 CHAPTER 6 SCALE-UP ANALYSIS 6.1 W Much work has been done at the laboratory—scale for the design of 3 Microwave Electrothermal Thruster (MET). However, the MET system is still not a mature technology. This is based primarily upon not much having been done to take this technology and scale it up to an operational thruster, performing the duties and tasks required of that thruster. As a result, this chapter has been added to highlight the tasks needed for that scale-up effort. These tasks, seven of them, were identified as scale—up issues. These issues were developed using information obtained from actual manufacturing experience, which could be applied to the scale-up effort for the MET system. An eighth item, six sigma, has been included to provide a systematic process towards addressing these seven issues. 6.2 Scale-up Issues. 6.2.1 Operability. Proving that the technology works in the laboratory is not enough. The MET system needs to be robust enough to be able to perform its tasks while in operation. The Operating limits need to be established and verified so that the system can operate with acceptable variance within the process. For example, we cannot expect to operate a MET system with a target of 400 torr with a variance of 10 torr that fails at Ill 395 torr. Also we cannot expect to operate the system at that same target pressure if the pressure instrumentation has an instrument error Of :l: 10 torr or greater. Basically, the MET system needs tO be analyzed thoroughly at expected operating conditions to ensure that the system will be robust enough tO Operate effectively when required. 6.2.2 Maintainability. What happens when the system fails, or a piece Of the system fails during operation while the thruster is in orbit? Does the entire satellite system fail as well? Mitigating actions need to be in place tO handle the common, or expected, mechanical / electrical failures. This includes using diagnostic equipment to determine if a system component has failed or 'will fail. During scale—up activities, plans should be developed to include methods or ways to maintain the MET system during operati on. 6.2.3 Controllability. Using the example Of operating the MET system at a target pressure Of 400 torr with an operating range Of 395 — 405 torr, one needs to determine the control strategy to ensure that this system is Operating within that required Operating range. This could be feed-forward, feed-back, or cascade types Of control systems using industrial-grade control systems, such as Programmable Logic Control (PLC) or Distributed Control System (DCS). The scale-up effort should include plans to develop and test the different control strategies. 6.2.4 Cost. Cost is always an issue for projects. One does not have an unlimited supply Of money or resources to develop the “perfect” system. An analysis 112 needs to be done to determine the cost-benefit for future research, including scale-up activities. The benefits for doing the research should clearly support the costs Of that research. For the MET system, this needs to be tied directly to the overall cost savings to a satellite system (including Operations) for using this new electric thruster system. Without this analysis, future funding for research of this MET system will not be expected. 6.2.5 Schedule. When does this scale—up effort need to be completed? After answering this question, plans should be developed tO ensure that this schedule is Obtained. Failure to Obtain the required schedule may result in not Obtaining the expected benefits Of the research. For example, if there were a potential Opportunity to attach a prototype MET system to a satellite platform being launched on a specific date in the future, the schedule needs to allow for the scale-up development to meet that date. Similar to cost, failure to meet this date may result in the MET system not being supported in the future. 6.2.6 Performance. Similar tO the operability issue, the MET system needs to Obtain a required level of performance. For example, the MET system may have a required specific impulse. If the system can meet that required performance level for no more than 20 hours, the system will not work for a system requirement Of a minimum Of 30 hours. The performance specifications should be Obtained from thruster and satellite vendors today, and should then be included into the scale-up effort. 113 6.2.7 Public Acceptance. Public acceptance is very important for the development Of the MET system. For example, if we Were to use nuclear power for the power source for the system, the system will not be supported if the public has a strong opposition to this type Of power in space. In essense, prior to designing the system for scale—up, one needs tO ensure that the various components Of the MET system comply with existing laws and regulation, and that the public has no strong Opposition to any Of those components. Failure to conduct this analysis early in the design effort may result in a wasted scale-up effort. . 6.2.8 Six Sigma. Six Sigma (60) is a rigorous implementation Of existing quality principles and techniques, which has a focus of eliminating errors or defects. Sigma (0') is the Greek letter used by statisticians to measure process variability, or deviations in the process. Historically, statisticians have found that processes shift up to 1.50 from the target. Using this process shift, the following sigma levels can produce the following associated defects per million (DPM): 114 Table 6.1 Sigma Significance Sigma Level ' DPM 1.50 500,000 2.00 308,300 2.50 158,650 3.00 67,000 3.50 22,700 4.00 6,220 4.50 1,350 5.00 233 5.50 J 32 6.00 3.4 Traditional manufacturing companies Operate at a 30 quality level. TO achieve 60 quality level Of performance, one does not Optimize the process to be twice as gOOd as 30 quality level. Instead, one has to become 20,000 times better. Using common situations, 30 quality level has the following results [Pyzdek, 1999]. Virtually nO modern computer would function. 10,800,000 healthcare claims would be mishandled each year. 18,900 US Savings bonds would be lost every month. 54,000 checks would be lost each night by a single large bank. 4,050 invoices would be sent out incorrectly each month by a modest- sized telecommunications company. 540,000 erroneous call details would be recorded each day from a regional telecommunications company. 270,000,000 erroneous credit card transactions would be recorded each year in the US. As can be seen by these numbers, this level of quality would not be acceptable for the MET system. For the design and fabrication Of a MET system, it would be wise to follow this proven process for scale-up, or one similar to it. 115 Using the normal distribution curve, an example Of two typical problems (large variation and Off-centered processes) that affect many processes are plotted in Figure 6.1. Two solutions are listed on that figure: 1) reduce the variation; and, 2) center the target Of the process. These curves are plotted within the lower specification limit (LSL) and upper specification limit (USL). The normal distribution, or the famous bell-shaped curve, has the following equation. ‘ l" ' l2 l '5 T f(x)= e for -oo ********************************************************* PARAMETER (NN = 10) REAL R(NN*NN,NN*NN), Z(NN*NN,NN*NN), B(NN*NN), T(NN*NN). + C1, C2, C3, C4, DR, DZ, A(NN*NN, NN*NN), RN INTEGER NR, NZ, I, J, N, X, IPATH, MN CHARACTER*8 FILE 1 PRINT *,’INPUT NAME OF FILEz’ READ 8, FILE 8 FORMAT (A8) OPEN (9, FILE = FILE, STATUS = ’NEW') IPATH = 1 PRINT *, ’HOW MANY R NODES?’ READ *, NR PRINT *, ’HOW MANY Z NODES?’ READ *, NZ N = NR * NZ DR = 64. / (NR+1.) DZ = 144. / (NZ+1.) C3 = 1./(DZ*DZ) C4 = -2./(DR*DR) - 2./(DZ*DZ) ****************************‘k'k'k‘kt'k'k'k'k‘k‘k‘k'k'k**************** * 'k * Set A matrix and B vector to zero. * * * **************************************** B(I) DO 6 J _ 1,N A(I,J) = o 6 CONTINUE 5 CONTINUE **** ****************************************************** 'k * * Set up A matrix and B vector. ‘k * * ***** ***************************************************** D010I=1,NR DO 20 J = 1, NZ R(I,J) = 25. + DR*I Z(I,J) = J * DZ C1 l./(DR*DR) + 1./(DR*R(I,J)*2.) C2 1./(DR*DR) - 1./(DR*R(I,J)*2.) = - * NZ + J IF(I1EQ11) B(X) = -C2 * (301.+95.*Z(I,J)/144.) IF(I.EQ.1) A(X,X+NZ) = C1 IF(I.EQ.NR) A(X,X-NZ) = C2 IF(I.EQ.NR) B(X) = -C1 * 301. IF(I.NE.1 .AND. I.NE.NR) THEN 162 A(X,X+Nz) = C1 A(X,X-NZ) = C2 ELSE ENDIF IF (J .EQ.1) B(X) = B(X) - C3*301. IF (J .EQ.NZ) B(X) = B(X) - C3*301. A(XIX) = C4 IF (J .EQ.1) A(X,X+1) = C3 IF (J .EQ.NZ) A(X,X-1) = C3 IF (J .NE.l .AND. J .NE.NZ) THEN A(X,X+1) = C3 A(X,X-1) = C3 ELSE ENDIF 20 CONTINUE lO CONTINUE ********************************************************** 'k t * Solve the A T = B problem. * * at ********************************************************** CALL GAUSS (A, B, T, N, MN, IPATH, RN) ******************************‘k'k************************** 'k 'k * Print out data. * ‘k t **************************‘k******************************* PRINT 500 D0 100 I = 1, NR+1 IF (I.EQ.NR+1) THEN PRINT 600, 25., 0., 301. WRITE (UNIT = 9, FMT = 600) 25., 0., 301. PRINT 600, 89.. 0., 301. WRITE (UNIT = 9, FMT = 600) 89., 0., 301. D0 50 J = 1, NZ RN = 301. + 95.*Z(l,J)/144. PRINT 600, 25., Z(1,J), RN WRITE (UNIT = 9, FMT = 600) 25., Z(l,J), RN PRINT 600, 89., Z(l,J), 301. WRITE (UNIT = 9, FMT = 600) 89., Z(1,J), 301. 50 CONTINUE DO 55 J = 1, NR PRINT 600, R(J,1), 0., 301. WRITE (UNIT = 9, FMT = 600) R(J,1), 0., 301. PRINT 600, R(J,1), 144., 301. WRITE (UNIT = 9, EMT = 600) R(J,1), 144., 301. 55 CONTINUE ELSE DO 110 J = 1, N2 x = (I-l) * NZ + J PRINT 600, R(I,J), Z(I,J), T(X) WRITE (UNIT = 9, PMT = 600) R(I,J), Z(I,J). T(X) 110 CONTINUE ENDIF 100 CONTINUE PRINT *,'IERROR =', IPATH . PRINT *,’DO YOU WANT ANOTHER SIMULATION (1 YES)? READ *, x IF (x .EQ. 1) GOTO 1 500 FORMAT(T5,’R’,T20,'Z’.T40,’T’) 163 600 FORMAT(T1,F8.4,T16,F8.4,T36,F8.3) END 164 A.4 COOLANT CHAMBER 165 * * t * * 6 5 ‘k * * ‘k * PROGRAM TWO ********************‘k************************************* * This program is designed to determine the temper— * ature profile of the air coolant in the discharge * cavity of the microwave electrothermal thruster * diagnostic chamber. The temperature profile is * known for the boundary of the region. This program * will use the centered finite difference method * to solve the heat equation for the axial and radial * dependant temperature assuming steady state and * angular independent conditions. * * 'k *t******************************************************* PARAMETER (NN = 10) B(NN*NN), T(NN*NN), RN REAL R(NN*NN,NN*NN), Z(NN*NN,NN*NN), C1, C2, C31, C32, C4, DR, DZ, A(NN*NN, NN*NN), INTEGER NR, NZ, I, J, N, X, IPATH, MN CHARACTER*8 FILE PRINT *,'INPUT NAME OF FILE2’ READ 8, FILE FORMAT (A8) OPEN (9, FILE = FILE, STATUS = ’NEW’) IPATH = l PRINT *, 'HOW MANY R NODES?’ READ *, NR PRINT *, 'HOW MANY Z NODES?’ READ *, NZ N = NR * NZ DR 7.1 / (NR+1.) DZ 144. / (NZ+1.) C31 l./(DZ*DZ) - 0.1785/DZ C32: l./(DZ*DZ) + 0.1786/DZ C4 = -2./(DR*DR) - 2./(DZ*DZ) ********************************************************* * * 4,. Set A matrix and B vector to zero. * ********************************************************* DO 5 I = 1,N B(I) = 0. DO 6 J = 1,N A(I,J) = 0. CONTINUE CONTINUE *‘k‘k****************************************************** * * Set up A matrix and B vector. * ********************************************************* DO 10 I = 1, NR DO 20 J = 1, NZ R(I,J) = 16.5+ DR*I Z(I,J) = J * DZ C1 = l./(DR*DR) + 1./(DR*R(I,J)*2.) C2 = 1./(DR*DR) - l./(DR*R(I,J)*2.) X = (I—l) * NZ + J IF(I.EQ.NR) B(X) = -C2 * IF(I.EQ.1) A(X,X+NZ) = C1 IF(I.EQ.NR) A(X,X-NZ) = C2 IF(I.EQ.1 .AND. Z(1,J).LE.108.) B(X) = -Cl* (300.+98.*Z(I,J)/l44.) 166 (300.+375.*Z(1,J)/108.) .. _C1* IF(I.EQ.1 .AND. Z(1,J).GT.108.) B(X) — (398.+277.*(l44.-Z(1,J))/36.) + IF(I.NE.1 .AND. I.NE.NR) THEN A(X,X+NZ) = Cl A(X,X-NZ) = C2 ELSE ENDIF IF (J .EQ.1) B(X) = B(X) - C32*300. IF (J .EQ.NZ) B(X) = B(X) - C3l*398. A(X,X) = C4 IF (J .EQ.1) A(X,X+l) = C31 IF (J .EQ.NZ) A(X,X—1) = C32 IF (J .NE.l .AND. J .NE.NZ) THEN A(X,X+1) = C31 A(X,X-1) = C32 ELSE ENDIF 20 CONTINUE 10 CONTINUE ********************************************************** * 'k * Solve the A T = B problem. * * *- ********************************************************** CALL GAUSS (A, B, T, N, MN, IPATH, RN) ********************************************************** 'k * * Print out data. * 'A' * ********************************************************** PRINT 500 D0 100 I = l, NR+1 IF (I.EQ.NR+1) THEN PRINT 600, 16.5.0., 300. 9, FMT = 600) 16.5,0., 300. WRITE (UNIT = PRINT 600, 23.6,0., 300. WRITE (UNIT = 9, FMT = 600) 23.6,0., 300. D0 50 J = 1, NZ IF(Z(1,J).GT.108.) RN = 398.+277*(l44.-Z(1,J)) /36. + IF(Z(1,J).LE.108.) RN = 300.+375.*Z(1,J)/108. PRINT 600, 16.5,Z(1,J), RN WRITE (UNIT = 9, FMT = 600) 16.5,Z(1,J), RN RN = 300. + 98.*Z(l,J)/144. PRINT 600, 23.6,Z(1,J), RN 9, FMT = 600) 23.6,Z(1,J), RN WRITE (UNIT = 50 CONTINUE DO 55 J = 1, NR PRINT 600, R(J,1), 0., 300. WRITE (UNIT = 9, FMT = 600) R(J,1), 0., 300, PRINT 600, R(J,1), 144., 398. WRITE (UNIT = 9, FMT = 600) R(J,1), 144 , 398, 55 CONTINUE ELSE DO 110 J = 1, NZ x = (I-l) * NZ + J PRINT 600, R(I,J). Z(I,J). T(X) WRITE (UNIT = 9, FMT = 600) R(I,J). Z(I,J). T(X) 110 CONTINUE 167 ENDIF 100 CONTINUE PRINT *,’IERROR =’, IPATH PRINT *,’DO YOU WANT ANOTHER SIMULATION (1 YES)?’ READ *, X IF (X .EQ. 1) GOTO 1 500 FORMAT(T5,’R’,T20,’Z’,T40,’T’) 600 FORMAT(T1,F8.4,T16,F8.4,T36,F8.3) END 168 A.5 DISCHARGE CHAMBER 169 PROGRAM THREE *******************************‘k'k'k************************* * 'k * This program is designed to determine the temper’ * * ature and velocity profile of the discharge in the * * cavity of the microwave electrothermal thruster * * diagnostic chamber. The temperature profile is * * known for the boundary of the region. This program * * will use the centered finite difference method * * to solve the momentum and energy equation for the * * axial and radial dependant temperature assuming * * steady state and angular independent conditions. * * 'k * Units: T — K * * P - Atm * * CP — BTU / mol K * * TC - BTU / m min K * * VIS - BTU min / m‘3 * * MW - lb / mol * * RHO, RT-lb / m“3 * * NE — # / cm‘3 * * R — mm * * Z - mm * * V — m / min * * * **i*****************************************~k************** PARAMETER (NN = 10, NT - 64) REAL*16 T(NN,NN),V(NN,NN),B(2*NT),RHO(NN,NN),RT(NN,NN,3), CP(NN,NN),NE(NN,NN),MW(NN,NN),VIS(NN,NN),P,XTOL,FR, A(2*NT,2*NT),G(2*NT),X(2*NT),RNORM,XERR,TC(NN,NN) REAL R(NN),Z(NN),PI,DR,DZ,RN INTEGER NR, NZ, I, J, K, L, N, IPATH, MN,MAXI,IERR + + IPATH = 1 FR = 0.5 * PRINT *, ’HOW MANY R NODES?’ * READ *, NR * PRINT *, ’HOW MANY Z NODES?’ * READ *, NZ P = 600./760. PI = 3.141593 NZ = 10 NR = 10 MN = 2*NT MAXI = 50 XTOL = 0.001 * N = NR * NZ DR = 16.5 / (NR-l.) DZ = 144. / (NZ-1.) ********************************************************** 'k * * Set up initial Temperature and Velocity profile. * ******************‘k‘k************************************** DO 10 I = 9,10 * * R(I) = (I-l) * DR Z(I) = (I-l) * Dz T(I,10) = 1100. + (l44.-Z(I))*300./36. l 0 CONTINUE 170 DO 20 I = 1,8 R(I) = (I-l) * DR Z(I) = (I-l) * DZ ' T(I,10) = 300. + 1100. * Z(I) / 108. 20 ‘ CONTINUE DO 30 I = 1,5 T(I,1) = 300. + Z(I) * 800. / 72. 30 CONTINUE DO 40 I = 5,10 T(I,l) = 7050. + 5950. * SIN(PI*(Z(I)-90.)/36.) 40 CONTINUE DO 50 J‘= 2,9 DO 50 I = 1,10 T(I,J) = (T(I,10)-T(I,1))*(J-1)/9. + T(I,l) 50 CONTINUE DO 60 I = 1,10 DO 60 J = 1,10 V(I,J) = 1.693 * T(I,J) / (760. * P) 60 CONTINUE ********************************************************** * 'k * Set A matrix and B vector to zero. * 'k * ********************************************************** DO 70 I = 1,MN B(I) = 0. DO 70 J = 1.MN A(I,J) = 0. 70 CONTINUE ********************************************************** * 'k * Solve using Newton Iteration. * * * ********************************************************** DO 100 K = 1, MAXI PRINT *, K DO 110 I = I,NN DO 110 J = I,NN ***********************************1?************************ * *- * Calculate the Density and Heat Capacity for each * point (along with viscosity and thermal conductivity 1' 'k *************************************************‘k********** CALL SM(T(I,J),P,RHO(I,J),CP(I,J),NE(I,J),MW(I,J).FR, RT(I,J,1),RT(I,J,2),RT(I,J,3)) + TC(I,J) = 9.875E—3 + 7.728E-6 * T(I,J) TC(I,J) = TC(I,J) + 1.391E-9 * T(I,J) ** 2 TC(I,J) = TC(I,J) - 1.9E-l3 * T(I,J) ** 3 TC(I,J) = TC(I,J) + 8.417E-18 * T(I,J) ** 4 VIS(I,J) = 4.47E-10 + 1.447E-11 * T(I,J) VIS(I,J) = VIS(I,J) + 1.895E15 * T(I,J) ** 2 VIS(I,J) = VIS(I,J) - 6.83E—20 * T(I,J) ** 3 VIS(I,J) = VIS(I,J) - 2.88E-24 * T(I,J) ** 4 110 CONTINUE ******‘k***************************************************** 171 * * * Set up the X vector. * * * ************************************************************ DO 120 I = 2,NN-l DO 120 J = 2,NN~1 LN I + (J-Z) * 8 * l X(LN) = T(I,J) G(LN) = X(LN) 120 CONTINUE DO 125 I = 2,NN-l LN = I + (J-2) * 8 + 64 - 1 X(LN) = V(I,J) G(LN) = X(LN) 125 CONTINUE ********************************************************** 'k * * Calculate the Jacobian, or A Matrix. * 'k 'k ********************************************************** 'k * Set up the VV area. 1' DO 130 I = 2,NN-1 DO 130 J = 2,NN—1 LN I + (J-2) * 8 + 64 — l A(LN,LN) = RHO(I,J) * (V(I+1,J)-V(I—1,J))/(2.*DZ) IF(I.LT.NN-l) A(LN+1,LN) = -VIS(I,J)/(DR*DR) - VIS(I,J)*1000./(2.*R(J)*DR) 1 IF(I.GT.2) A(LN-1,LN) = -VIS(I,J)/(DR*DR) - + VIS(I,J)*1000./(2.*R(J)*DR) IF(J.LT.NN—l) A(LN+8,LN) = RHO(I,J)*V(I,J)/(2.*DZ) - + VIS(I,J)/(DZ*DZ) IF(J.GT.2) A(LN-8,LN) = -RHO(I,J)*V(I,J)/(2.*DZ) — + VIS(I,J)/(DZ*DZ) 130 CONTINUE * * Set up the TV area. * DO 135 I = 2,NN—1 DO 135 J = 2,NN-l LN I + (J-2) * 8 + 64 - 1 A(LN-64,LN)= RHO(I,J)*CP(I,J)*MW(I,J)*(T(I+l,J)-T(I-1.J))/ (2.*DZ) - 254742. * RT(I+1,J,1) / (2. * DZ) + + IF(J.LT.NN-l)A(LN-56,LN) - 43079. * RT(I+1,J,2) / (2. * DZ) + + + 87685. * RT(I+1,J,3) / (2. * DZ) IF(J.GT.2) A(LN—72,LN) = —254742. * RT(I—1,J,l) / (2. * DZ) + + -43079. * RT(I-1,J,2) / (2. * DZ) + -87685. * RT(I-1,J,3) / (2. * DZ) 135 CONTINUE 172 * Set up TT area. * DO 140 I = 2,NN-l DO 140 J = 2,NN-1 LN = I + (J-2) * 8 - 1 A(LN,LN) = 2.*TC(I,J)/(DZ*DZ) +2.*TC(I,J)/(DR*DR) IF(I.LT.NN-l) A(LN+1,LN) = - TC(I,J)/(DR*DR) + - TC(I,J)*1000./(2.*R(J)*DR) IF(I.GT.2) A(LN-1,LN) = - TC(I,J)/(DR*DR) + + TC(I,J)*1000./(2.*R(J)*DR) IF(J.LT.NN—l) A(LN+8,LN) = RHO(I,J)*CP(I,J)*MW(I,J)*V(I,J)/ + (2.*DZ) - TC(I,J)/(DZ*DZ) IF(J.GT.2) A(LN—8,LN) = - RHO(I,J)*CP(I,J)*MW(I,J)*V(I,J)/ + (2.*DZ) - TC(I,J)/(DZ*DZ) 140 CONTINUE ***************************************************************** * 'k * Set up the B vector. * * * ***************************************************************** DO 150 I = 2,NN-1 DO 150 J = 2,NN-1 LN = I + (J-2) * 8 — 1 B(LN) = RHO(I,J)*CP(I,J)*MW(I,J)*V(I,J)*(T(I+1,J) + -T(I-l,J))/(2.*DZ) B(LN) = B(LN) - TC(I,J)*(T(I,J+1)+T(I,J-1)- + 2.*T(I,J))/(DR*DR) B(LN) = B(LN) - TC(I,J)*(T(I,J+1)-T(I,J—1))/ + (2.*DR) B(LN) = B(LN) - TC(I,J)*(T(I+l,J)+T(I-1,J)- + 2.*T(I,J))/(DZ*DZ) B(LN) = B(LN) - 254742. * (V(I+1,J)*RT(I+1,J,1)- + V(I-1,J)*RT(I-1,J,l)) / (2. * DZ) B(LN) = B(LN) - 43079. * (V(I+1,J)*RT(I+1,J,2)- + V(I-l,J)*RT(I-1,J,2)) / (2. * DZ) B(LN) = B(LN) - 87685. * (V(I+1,J)*RT(I+1,J,3)- + V(I-l,J)*RT(I-1,J,3)) / (2. * DZ) IF(NE(I,J) .GE. 1.E10) B(LN) = B(LN) + 1336000. 150 CONTINUE DO 155 I = 2,NN-l DO 155 J 2,NN-1 LN = I + (J-2) * 8 + 64 — 1 B(LN) = RHO(I,J)*V(I,J)*(V(I+1,J)-V(I-1,J))/ + (2*DZ) B(LN) = B(LN) - VIS(I,J)*(V(I,J+1)+V(I,J—l)- + 2.*V(I,J))/(DR*DR) B(LN) = B(LN) - VIS(I,J)*(V(I,J+1)-V(I,J-1))* + 1000./(2.*R(J)*DR) B(LN) = B(LN) - VIS(I,J)*(V(I+1,J)+V(I-1,J)- + 2.*V(I,J))/(DZ*DZ) 155 CONTINUE 173 ***************************************************************** * 'k * Solve for the linearized Ax = B. * * * ********************************‘k'k'k****************************** CALL GAUSS(A,B,G,MN,MN,IERR,RNORM) DO 160 L = 1,MN X(L) = X(L) + G(L) IF(X(L) .LT. 0.) X(L) = —X(L) 160 CONTINUE * XERR = 0. * DO 170 L = l, * XERR = XERR + (X(L)-G(L))**2 *170 CONTINUE * XERR = XERR ** (0.5) * IF(XERR .LE. XTOL) GOTO 200 100 CONTINUE 200 CONTINUE ****************************‘k***i************************** * 'k * Write out the Temperature Profile. * * * ******************‘k**************************************** OPEN (1, FILE = ’D6B.T' , STATUS = ’NEW’) OPEN (2, FILE = 'DGB.V’ , STATUS = ’NEW’) OPEN (3, FILE = ’DGB.NE’, STATUS = ’NEW') DO 300 I = 1,NN DO 300 J = 1,NN WRITE (UNIT=1, FMT=900) R(J),Z(I),T(I,J) PRINT *, R(J),Z(I),T(I,J) 300 CONTINUE *****************************1k***************************** * i . . . * * Print out the velOCity profile. * * *********************************************************** DO 310 I = 1,NN DO 310 J = 1,NN WRITE (UNIT=2, FMT=900) R(J), Z(I), V(I,J) PRINT *, R(J), Z(I), V(I,J) 310 CONTINUE *********************************************************** * * * Print out the electron density profile. * t * *‘t‘k ******************************************************** DO 320 I 1,NN DO 320 J 1,NN WRITE (UNIT=3, FMT=900) R(J), Z(I), NE(I,J) PRINT *, R(J), Z(I), NE(I,J) 320 CONTINUE 900 FORMAT (2F7.2,T20,E9.3) END 174 A.6 STATISTICAL MECHANICS 175 PROGRAM ONE *********************************************************** * * * This program is designed to determine the mole * * fraction and the thermodynamic properties of the * * helium and nitrogen mixture plasma fluid using * * statistical mechanics (partition functions) using * «k the Newton Method to solve the non-linear set of * * equations. This program calls in previously * * defined files of energy levels for all species * * (ions and neutrals) being considered. * * * * VARIABLES: * * E(l) = Total Energy of Fluid * * E(2) = Enthalpy of Fluid * * E(3) = Entropy of Fluid * * E(4) = Chemical Potential * * E(S) = Heat Capacity * * EA() = Energy of individual species * * EI() = Energy of helium * * EII() = Energy of helium (+1) * * F(1-8) = The equations. * * FERR = Function error * * FTOL = Function tolerance * * FR = Mole Fraction of N2 to He * * G(l—8) = Second iteration of variables. * * H = Planks constant * * J(x,y) = The Jacobian of the four eqns. * * K(l-5) = Equalibrium constants * * KB = Boltzmand constant * * M(l-8) = The mass of the partical * * MAXI = Maximum number of iterations. * * NE = Number Electron Density * * P = The pressure * * RHO = Density Ratio * * Q(1-8) = Electronic partition function * * T = The temperature of the electron * * WI() = Degeneracy of helium * * WII() = Degeneracy of helium (+1) * * X(l) = The mole fraction of electrons. * * X(2) = m.f. of helium * * X(3) = m.f. of helium (+1) ion * * X(4) = m.f. of helium (+2) ion * * X(S) = m.f. of nitrogen * * X(6) = m.f. of nitrogen (+1) ion * ' X(7) = m.f. of nitrogen (+2) ion * * X(8) = m.f. of nitrogen (+3) ion * * XERR = Variable error * * XP(1-8)= Partial of X() * * XTOL = Variable tolerance * * Z = Compressibility Factor * ”k * *********************************************************** DOUBLE PRECISION X(8),G(8),F(8),J(8,8),T,M(8),P(2),Q(8),RHO, + XTOL,FTOL,EI(2,36),EII(4,36),XERR,FERR,E(5),EA(8),XP(8),Il, + PI,BETA,RNORM,K(5),KB,H,N,NE,Z,NI(36),NII(36),NIII(21),IP(5) INTEGER I,L,WI(36),WII(21),WWI(36),WWII(36),MAXI,IERR,MN, + WWIII(36), WWIIII(21) PRINT *,’TYPE IN MOLE FRACTION OF NITROGEN:' 176 READ *, FR * PRINT *,'TYPE IN TEMPERATUREz' * READ *, T P(l) = 1.013E6. PRINT *,’TYPE IN PRESSURE (atm):’ READ *, Il P(2) = P(l) * 11 T = 10000 *********************************************************** * * * Set constants. * 'k * *********************************************************** N = 6.022E23 H = 6.6262E-27 PI = 3.141592654 M(l) = 9.1095E-28 M(2) = 6.6473E-24 M(3) = M(2) - M(l) M(4) = M(3) - M(l) M(S) = 2.3261E-23 M(6) = M(S) - M(l) M(7) = M(6) - M(l) M(8) = M(7) - M(l) KB = 1.3806E-16 MAXI = 100 XTOL = 1.D-50 FTOL = 1.D-50 *********************************************************** * * * Read in energy levels for helium and helium (+1). * 'k * *******************************t'k************************** OPEN (7. FILE = ’HEI.EXC', STATUS = 'OLD’) OPEN (8, FILE = ’HEII.EXC’, STATUS = ’OLD’) DO 10 I = 1,36 READ (7,*) WI(I), EI(1,I) 10 CONTINUE CLOSE (7) DO 20 I = 1,21 READ (8,*) WII(I), EII(1,I) 20 CONTINUE CLOSE (8) *********** ***********************************r************ * * 177 * Read in energy levels for nitrogen and nitrogen (+1)* * * *********************************************************** OPEN (6, FILE 'NI.EXC’, STATUS = ’OLD') OPEN (7, FILE 'NII.EXC’, STATUS = ’OLD’) OPEN (8, FILE 'NIII.EXC', STATUS = ’OLD’) OPEN (9, FILE 'NIIII.EXC', STATUS = ’OLD') DO 15 I = 1,36 READ (6,*) WWI(I), EI(2,I) 15 CONTINUE CLOSE (6) DO 25 I = 1,36 READ (7,*) WWII(I), EII(2,I) 25 CONTINUE CLOSE (7) DO 28 I = 1,36 READ (8,*) WWIII(I), EII(3,I) 28 CONTINUE CLOSE (8) DO 29 I = 1,21 READ (9,*) WWIIII(I), EII(4,I) 29 CONTINUE CLOSE (9) *****************************************************‘k***** 'k 'k * Set initial values for X() and XP(). * 'k * *********************************************************** X(l) = 1.E-5 X(2) = (l.-FR) * 1. X(3) = (l.-FR) * 1.E-5 X(4) = (l.-FR) * 1.E-20 X(S) = FR * 1. X(6) = FR * 1.E-5 X(7) = FR * 1.E—20 X(8) = FR * l.D-25 XP(1) = 1.E—3 XP(2) = —l.E-3 XP(3) = 1.E-3 XP(4) = .1E-3 XP(5) = -1.E-3 XP(6) = 1.E-3 XP(7) = .1E-3 XP(8) = .01E-3 ************************************************ «k at * OPEN FILES FOR PRINTING AND BEGIN * ITERATIONS / CALCULATIONS. * 'k ‘k * ************************************************ PRINT 901 PRINT 902 1 IP(l) 198305. IP(2) 438900. IP(3) 117345. IP(4) 238847. IP(S) 382626. T = T 200. *********** ************************************************ + HIIIIN H t * ' * * Solve for electronic partition function. 178 * * **********************************************************'k EA(l) EA(2) EA(3) EA(4) EA(S) EA(6) EA(7) EA(8) Q(1) Q(2) Q(3) Q(4) Q(S) Q(6) Q(7) Q(8) II II II II II II II II OOOOOOOO II II I! II I) II II II OOOOHOON * HELIUM DO 30 I = 1,36 BETA = -1.98705E—16 * EI(1,I) / (KB*T) Q(2) = Q(2) + WI(I) * EXP(BETA) EA(2) = EA(2) - WI(I) * BETA * EXP(BETA) * T / 273.15 30 CONTINUE EA(2) = EA(2) / Q(2) DO 40 I = 1,21 BETA —l.98705E-l6 * EII(1,I) / (KB*T) Q(3) Q(3) + WII(I) * EXP(BETA) EA(3) = EA(3) - WII(I) * BETA * EXP(BETA) * T / 273.15 40 CONTINUE EA(3) = EA(3) / Q(3) * NITROGEN DO 35 I = 1,36 BETA = -1.98705E-16 * EI(2,I) / (KB*T) Q(S) = Q(S) + WWI(I) * EXP(BETA) EA(S) = EA(S) - WWI(I) * BETA * EXP(BETA) * T / 273.15 35 CONTINUE EA(S) = EA(S) / Q(S) DO 45 I = 1,36 BETA = -1.98705E-16 * EII(2,I) / (KB*T) Q(6) = Q(6) + WWII(I) * EXP(BETA) EA(6) = EA(6) - WWII(I) * BETA * EXP(BETA) * T / 273.15 45 CONTINUE EA(6) = EA(6) / Q(6) DO 48 I = 1,36 BETA = -1.98705E-16 * EII(3,I) / (KB*T) Q(7) = Q(7) + WWIII(I) * EXP(BETA) EA(7) = EA(7) - WWIII(I) * BETA * EXP(BETA) 48 CONTINUE EA(7) = EA(7) / Q(7) DO 49 I = 1,21 * BETA -1.98705E-16 * EII(4,I) / (KB T) * T / 273.15 f * A Q(8) - Q(8) + WWIIII(I) EXP(BET ) , T / 273.15 EA(8) = EA(8) - WWIIII(I) * BETA * EXP(BETA) 179 49 CONTINUE EA(8) = EA(8) / Q(3) *********************************************************** 'k i * Solve for individual energy of species. * 'k * ******‘k**************************************************** DO 50 I 1.7 EA(I) EA(I) + 3. * T / (2. * 273.15) 50 CONTINUE *i’********************************************************* * * * Solve for mole fraction ground state and K's. * 'k * *‘k********************************************************* * HELIUM IP(1) -1.98705E-16 * IP(1) / (KB * T) IP(1) = EXP (IP(1)) IP(2) = -1.98705E—16 * IP(2) / (KB * T) IP(2) = EXP (IP(2)) BETA = 2. * PI * M(1) * KB * T / (H*H) K(1) = BETA ** (1.5)*KB * T * IP(1) * Q(1) * Q(3) / Q(2) K(2) = BETA ** (1.5)*KB * T * IP(2) * Q(1) * Q(4) / Q(3) * NITROGEN IP(3) = -1.98705E—16 * IP(3) / (KB * T) IP(3) = EXP (IP(3)) IP(4) = -1.98705E—16 * IP(4) / (KB * T) IP(4) = EXP (IP(4)) IP(S) = -l.98705E-16 * IP(S) / (KB*T) IP(S) = EXP(IP(5)) BETA = 2. * PI * M(1) * KB * T / (H*H) K(3) = BETA ** (1.5)*KB * T * IP(3) * Q(1) * Q(6) / Q(S) K(4) = BETA ** (1.5)*KB * T * IP(4) * Q(1) * Q(7) / Q(6) K(S) = BETA ** (1.5)*KB * T * IP(S) * Q(1) * Q(8) / Q(7) * PRINT *, K * PRINT *, Q **** ******************************************************* * 'k * Solve for mole fraction using Newton Iteration. 'k * * ******* **************************************************** DO 100 I = 1, MAXI '- 2.*(X(4)+ X(7)) + 3.*X(8) F(l) X(l)+X(3)+X(6)+ )-X(6)-X(7)-X(8)+l. P(2) ;-X(l)-X(2)-X(3)-X(4)-X(5 F(3) =-X(l)*X(3)/X(2)+K(1)/P(2) F(4) =-X(l)*X(4)/X(3)+K(2)/P(2) 180 110 "I on V II II II F(7) 0 L 1,7 CONTINUE FERR = (FERR) ** (0.5) =—X(1)*X(6)/X(5)+K(3)/P(2) =-X(1)*X(7)/X(6)+K(4)/P(2) =(FR-l.)*(X(5)+X(6)+X(7)+X(8)) + FR*(X(2)+X(3)+X(4)) (l)*X(8)/X(7)+K(S)/P(2) FERR + F(L)*F(L) IF (FERR .LE. FTOL) GOTO 200 1 J(l,1) J(1,2) J(1,3) J(1,4) J(l,5) J(1,6) J(1,7) J(1,8) J(2,1) J(2,2) J(2,3) J(2,4) J(2,5) J(2,6) J(2,7) J(2,8) J(3,1) J(3,2) J(3,3) J(3,4) J(3,5) J(3,6) J(3,7) J(3,8) J(4,l) J(4,2) J(4,3) J(4,4) J(4,5) J(4,6) J(4,7) J(4,8) J(5,l) J(5,2) J(5,3) J(5,4) J(5,5) J(5,6) J(5,7) J(5,8) J(6,1) J(6,2) J(6,3) J(6,4) J(6,5) J(6,6) J(6,7) J(6,8) J(7,l) J(7,2) 0. -1. -2. X(3) / X(2) —X(l) * X(3) X(l) / X(2) (4) / X(3) OXOOOOO -'(1) * X(4) X(l) xxm (6) / X(5) o<3c3x<3c>o<3 Q(1) * X(6) X(l) / X(5) (7) / X(6) OOOOXOO - (1) * X(7) X(l) / X(6) 0. 0. —FR (X(2)*X(2)) (X(3)*X(3)) (X(5)*X(5)) (X(6)*X(6)) 181 J(7,3) J(7.4) J(7,S) J(7.6) J(7,7) J(7,8) J(8,1) J(8,2) J(8,3) J(8.4) J(8,5) J(8,6) J(8,7) J(8,8) DO 115 G(L) 115 CONTINUE MN = 8 CALL GAUSS(J, F, X, 8, MN, IERR, RNORM) DO 116 L = 1,8 X(L) = X(L) + G(L) IF (X(L) .LT. 0.) X(L) = -X(L) 116 CONTINUE XERR = 0. DO 120 L = 1,8 XERR = XERR + (X(L) - G(L)) ** 2 120 CONTINUE XERR = (XERR) ** (0.5) IF (XERR .LE. XTOL) GOTO 200 100 CONTINUE 200 CONTINUE -FR -FR 1. - FR 1. - FR 1. - FR 1. - FR X(8) / X(7) OOOOO -X(1) * X(8) / (X(7)*X(7)) X(l) / X(7) - 1,8 X(L) fill" M "II" "II" H "II" M * PRINT *,x * PRINT *,F * PRINT *,Q ********************************‘k‘k‘k'k'k'k'k‘k'k'k'k‘k'k'k'k‘k********** * ‘k * Calculate the electron density and compressibility * * factor and density ratio. * 'k * ********************************************************** NE = X(l) * P(2) / (KB * T) Z (1.-FR) * M(2) + FR * M(S) z z / (M(1)*X(l) + M(2)*X(2) + M(3)*X(3) + M(4)*X<4> + + M(5)*X(5) + M(6)*X(6) + M(7)*X(7> + M(8)*X(8>) RHO = P(2) * 273.15 / (P(l) * T * Z) * ******************************r************************** * * * * Calculate energy and enthalpy. * *- ***** ***************************************************** E(l) E(2) E(3) E(4) E(S) ll HOOOOO = 1,8 3(1) + Z * X(I) * EA(I) 182 210 CONTINUE E(2) = E(l) + T / 273.15 *********************************************************** Solve for partial of mole fraction with respect to temperature at constant pressure using Newton Iteration. ***** ##II'II’N’ *********************************************************** DO 300 I = 1, MAXI F(1) =-XP(1) + XP(3) + XP(6) + 3. * XP(8) F(1) = F(1) + 2. * (XP(4) + XP(7)) P(2) =-XP(1) - XP(2) — XP(3) - XP(4) F(2) = P(2) - XP(S) - XP(6) - XP(7) - XP(8) F(3) =-XP(1)/X(1) + XP(2)/X(2) - XP(3)/X(3) F(3) = F(3) + EA(3)*273.lS/T/T + l./T F(3) = F(3) + EA(1)*273.15/T/T + l./T F(3) = F(3) - EA(2)*273.15/T/T - l./T F(4) =-XP(1)/X(l) + XP(3)/X(3) — XP(4)/X(4) F(4) = F(4) + EA(4)*273.15/T/T + 1./T F(4) = F(4) + EA(1)*273.15/T/T + 1./T F(4) = F(4) - EA(3)*273.lS/T/T - 1./T F(S) =-XP(l)/X(1) + XP(5)/X(5) — XP(6)/X(6) F(5) = F(S) + EA(6)*273.15/T/T + 1./T F(S) = F(5) + EA(1)*273.15/T/T + 1./T F(5) = F(S) - EA(5)*273.15/T/T - 1./T F(6) =-XP(l)/X(l) + XP(6)/X(6) - XP(7)/X(7) F(6) = F(6) + EA(7)*273.15/T/T + 1./T F(6) = F(6) + EA(1)*273.lS/T/T + 1./T F(6) = F(6) - EA(6)*273.15/T/T - 1./T F(7) =(FR-1.)*(XP(5)+XP(6)+XP(7)+XP(8)) F(7) = F(7) + FR*(XP(2)+XP(3)+XP(4)) F(8) = -XP(1)/X(l) + XP(7)/X(7) - XP(8)/X(8) F(8) = F(8) + EA(8)*273.15/T/T + 1./T F(8) = F(8) + EA(1)*273.lS/T/T + 1./T F(8) = F(8) - EA(7)*273.15/T/T + 1./T FERR = 0. DO 310 L = 1,7 FERR = FERR + F(L)*F(L) 310 CONTINUE FERR = (FERR) ** (0.5) IF (FERR .LE. FTOL) GOTO 400 J(l,1) l. J(1,2) 0. J(1,3) —1. J(1,4) -2. J(1,5) J(1,6) J(1,7) J(1,8) J(2,1) J(2,2) J(2,3) J(2,4) J(2,5) J(2,6) J(2,7) J(2,8) J(3,1) Z/X(l) 183 315 316 J(3,2) J(3,3) J(3,4) J(3,5) J(3,6) J(3,7) J(3,8) J(4,1) J(4,2) J(4,3) J(4,4) J(4,5) J(4,6) J(4,7) J(4,8) J(5,1) J(5,2) J(5,3) J(5,4) J(5,5) J(5,6) J(5,7) J(5,8) J(6,1) J(6,2) J(6,3) J(6,4) J(6,5) J(6,6) J(6,7) J(6,8) J(7,1) J(7,2) J(7,3) J(7,4) J(7,5) J(7,6) J(7,7) J(7,8) J(8,l) J(8,2) J(8,3) J(8,4) J(8,5) J(8,6) J(8,7) J(8,8) DO 315 G(L) CONTINUE -1./X(2) 1./X(3) :/x(1) i./x<3) ./X(4) i/xm 1./X(S) ./X(6) I/xu) OOOOHOOP—‘I OOOI—‘OOOOI—‘l CHOOOOO OCDC)O He+ + e- * * 'A' * * * k2 * * He+ + e- --------- > He * * * * * * * * NOTE: At equilibrium, kl = k2 * * * * * * Units: T = K * * P = Atm * * Kl = Ionization Rate * * K2 = Recombination Rate * * E0 = Ionization Energy * * SIGMAO = Ionization Cross Section Area * * NE = # / cm“3 * * EI() = Energy of helium * * H = Planks constant * * K(l—2) = Equalibrium constants * * KB = Boltzmand constant * * M(1-3) = The mass of the partical * * Q(1-3) = Electronic partition function * * WI() = Degeneracy of helium * * X(l) = The mole fraction of electrons. * * X(2) = m.f. of helium * * X(3) = m.f. of helium (+1) ion * * * *********************************************************** DOUBLE PRECISION T,RHO,NE,EI(36),EII(21),K(2),I2,F(3), + N,PI,BETA,KB,M(3),I1,RATE1,K1,K2,E0,SIGMAO,XTOL,FTOL, + EA(3),P(2),Q(3),TIME,PRESS,J(3,3),X(3),G(3),XERR,RNORM INTEGER I, L, WI(36), WII(21), MAXI, IERR, MN OPEN (4, FILE ’FORWARD1.RXN’, STATUS = ’NEW’) OPEN (5, FILE ’REVERSE1.RXN', STATUS = ’NEW’) OPEN (6, FILE ’TIME_OF1.RXN’, STATUS = ’NEW') ’OLD’) ’OLD') ’HEI.EXC', STATUS ’HEII.EXC’,STATUS OPEN (7, FILE OPEN (8, FILE *********************************************************** t t * * Set constants. * *********************************************************** N = 6.022E23 PI = 3.141592654 M(1) = 9.1095E-28 M(2) = 6.6473E-24 M(3) = M(2) - M(1) 191 KB 1.3806E-l6 E0 5. * 1.6022E—l9 SIGMAO = 7.E-l4 T = 8000 PRESS = .4 P(l) 101325 P(2) P(1) * PRESS *****‘k‘k‘k'k'k************************************************* * t * Initiate values. * 'k * it*‘k'kt****************************************************'k X(l) = 0.00001 X(2) = 1.00000 X(3) = 0.00001 MAXI = 100 XERR = 1.E-5 XTOL = 1.E-5 FTOL = 1.E—S *********************************************************** 'k . 'k * Read in energy levels for helium and helium+1. * ‘k * *********************************************************** DO 10 I = 1,36 READ (7,*) WI(I), EI(I) 10 CONTINUE CLOSE (7) DO 20 I = 1,21 READ (8,*) WII(I). EII(I) 20 CONTINUE CLOSE (8) Il = 198305. IZ = 438900 ******** ********************trt**************************** * 'k * Solve for electronic partition function. * * * ******** *************************************************** Q(1) = 2. 0(2) = 0. Q(3) = 0. BETA = —1.98705E—16 * EI(I) / (KB*T) Q(2) = 0(2) + WI(I) * EXP(BETA) 30 CONTINUE DO 40 I = 1,21 BETA -l.98705E—l6 * EII(I) / (KB*T) Q(3) Q(3) + WII(I) * EXP(BETA) 40 CONTINUE 192 L...~ *********************************************************** 'k 'k * Solve for mole fraction ground state and K's. * 'k * *********************************************************** I1 = -l.98705E-16 * Il / (KB * T) I1 = EXP (I1) BETA = 1.8E+10 * T K(l) = BETA ** (1.5) * KB * T * I1 * Q(1) * Q(3) K(l) = K(l) / Q(2) PRINT *, K(l), BETA, KB, T, I1, Q(1), Q(3), Q(2) ***‘k'k*********************************************i'k'k'k‘k'kt'k‘k * * * Calculate the reaction rate for the ionization. * * * *********************************************************** KB = KB * 1.E-7 RATEl = ( 8. * KB * T / ( M(2) * PI ) ) ** 0.5 RATEl = RATEl * ( 1. + E0 / ( KB * T ) ) * SIGMAO RATEl = RATEl * EXP( - E0 / ( KB * T ) ) RATEl = RATEl * 1.E6 ************************'k********************************** * * * Calculate the particle density at ideal conditions. * t * *********************************************************** RHO = 7.24E16 * P(2) / T PRINT *, RHO *********************************************************** * * * Solve for mole fraction using Newton Iteration. 'k 'k * * ********************************************************** F(1) =-X(1) + X(3) P(2) =-X(l) - X(2) - X(3) + 1. F(3) =-X(1) * X(3) / X(2) + K(l) / P(2) FERR = 0. DO 110 L = 1,3 FERR = FERR + F(L)*F(L) 110 CONTINUE FERR = (FERR) ** (0.5) IF (FERR .LE. FTOL) GOTO 190 J(l,1) = 1. J(1,2) = O. J(1,3) = —1 J(2,1) = 1. J(2,2) = 1. 193 J(2,3) = 1. J(3,1) = X(3) / X(2) J(3,2) = -X(1) * X(3) / (X(2)*X(2)) J(3,3) = X(l) / X(2) DO 115 L = 1,3 G(L) = X(L) 115 CONTINUE MN = 3 CALL GAUSS(J, F, X, 3, MN, IERR, RNORM) DO 116 L = 1,3 X(L) = X(L) + G(L) IF (X(L) .LT. 0.) X(L) = -X(L) 116 CONTINUE XERR = 0. DO 120 L = 1,3 XERR = XERR + (X(L) - G(L)) ** 2 120 CONTINUE XERR = (XERR) ** (0.5) IF (XERR .LE. XTOL) GOTO 190 100 CONTINUE 190 CONTINUE *********************************************‘k'k'k‘kr‘k'k'k‘t'ki't'k'k * 'k * Calculate the reaction rate for the recombination * * using the equilibrium conditions. * * * *********************************************************** RATE2 = X(2) * RATEl / X(l) *********************************************************** 'k i' * Vary the electron density and calculate both the * ionization rate and recombination rate. * * * * *********************************************************** NE = 0. DO 300 I = 1,100 NE = I * RHO / 1.E+6 K1 = ( RHO - 2. * NE ) * RATEl K2 = NE * RATE2 WRITE (UNIT = 4, FMT = *) NE, K1 WRITE (UNIT = 5, FMT = *) NE. K2 300 CONTINUE ************************ * *********************************** 'k * Calculate the electron density formation vs. * time. * * ****** *****************r*********************************** TIME = 0. NE = 0. C2 =-2. * RATEl - RATE2 C2 = C2 / 1.E-6 DO 400 I = 1,100 194 TIME = I * 1.E-5 NE = 1. - EXP (C2 * TIME) WRITE (UNIT = 6, FMT = *) TIME, NE 400 CONTINUE PRINT *, RATEl, K(l) END 195 ANNEX B ATOMIC ENERGY LEVELS 196 B.l HELIUM 197 He Atomic Energy Levels He He+ (0 Energy Levels (cm'l) (0 Energy Levels (cm’l) 1 0 2 0 3 159850.318 2 329179.102 1 166271.70 2 329179.572 5 169081.111 4 329184.945 3 169081.189 J 2 390140.622 1 169082.185 J 2 390140.761 .1 3 171129.148 4 3901424 3 183231.08 4 3901424 1 184859.06 6 3901429 5 185558.92 32 4114770 3 185559.085 50 4213530 1 185559.277 72 4267170 7 186095.90 72 4299516 F 5 186095.90 72 4320509 3 186095.90 72 4334902 L 5 186099.22 72 4345197 3 190292.46 j 72 435281.4- 64 190900 ] 72 4358610 64 193600 ] 72 4363116 64 195200 j 72 436669.4 64 195950 1 72 4369575 64 196656 IONIZATION 438908.670 64 196941 J 57 197175 I 56 197375 1 36 197525 36 197640 | 35 197743 30 197815 L 27 197815 L 27 197924 L 27 197965 L 27 198000 I: 27 198030 ] 24 198056 J L 9 198077 IONIZATION 198305 198 B.2 NITROGEN 199 N Atomic Energy Levels N+ Energy Levels (cm-1) 0 15315.7 32687.1 47167.7 92250 109220 1441891 149000 149188.74 1551299 164611.60 166600 1667657 168893.04 170600 174212.93 178274.17 1 86600 187092.20 187460 188900 1893360 190121.15 196600 2021699 202800 203200 2035328 20553507 206000 JET—”T 2063275 209750 209926.92 210270 210750 Oxo’J‘mtgv—e31116303630uxlxoamgr—mxomwawwwxomxoGMHuee N (0 Energy Levels (cm‘l) 4 0 10 19227 6 28840 12 83320 6 86176.8 12 88160 2 93582.3 20 94800 12 95500 4 96751.7 10 96825 6 97785 10 99660 12 103680 6 104185 6 104630 28 104700 14 104850 12 104900 20 105000 10 105130 2 10647 8.6 20 106800 12 107000 4 107447 .2 12 109900 4 1 10050 28 110250 20 1 10300 4 1 10230 14 1 10350 12 1 10375 10 1 10460 10 1 10530 4 112310 12 112600 IONIZATION 1 17345 211030.90 200 IONIZATION 2388467 N+2 0) Energy Levels (cm‘l) 6 0 12 57250 10 101027 2 13100335 6 145900 4 1868023 6 230406 6 _, 245680 10 267240 12 287600 6 297200 2 3010882 I 6 309160 20 309700 6 311700 4 3142240 12 317765 FIDO 14 320287 10 321000 2 3270568 28 330300 20 332820 I 2 3337131. 10 334550 ] 12 336270 J 14 339800 10 341947 | 6 342700 14 3427520 18 343116 10 354517 12 3549557 18 355214 12 368600 10 373360 [:_ 6 374775 IONIZATION 3826255 201 N+3 Energy Levels (cm'l) 0 67200 130695 175500 188885 235370 377206 388858 404521 405900 419970 ‘ 429158 465300 473032 480880 484450 487542 494300 48315 MMOWGWWWMGOWer—LAOWOP—e 499708 1*T‘n—T—C'7-T— 21 499851 ICHNELAIICHV 623851 N+4 N+5 0) Energy Levels (cm‘l) (0 Energy Levels (cm'l) 2 0 1 o 6 80600 3 3385890 2 456134 9 3438300 6 477800 [g 3 3473790 10 484415 3 4016390 2 606337 3 4206810 6 615150 4452800 10 617905 2 673882 [7 6 678297 10 679725 2 709947 6 712464 ! 10 713289 N+5 14 713327 (1) Energy Levels (cm-1) 30 713335 2 0 2 731432 2 4034535 6 732993 2 4034605 10 733516 4 4035412 _. 14 733547 64 4782200 44 733552 _] 64 5043700 I IONIZATION 7 89532.9 1 1 IONIZATION 5379860 I 202 APPENDIX C A TWO-DIMENSIONAL KINETICS PROGRAM SIMULATION 203 same name ooze ozma .omnm: mz<¢ea ozmq .m«.o.om.o.mn.on:x.flumma>fiu .o.mu9m2 >EHU zomzdu .oom wPHDm .Emmmfim 2m<3mmmm .ZOHmmm> and .AXQE. 2enmzmo oo+moooo.ouonemaom mo+mooo~.ouamau ezmommm oo+moooo.ouE\Oo 0000.0 00.0 0 000.0 000H0.0 00000.~ z amam oooo.o oo.o o ooo.o oooam.o cocoo.H m: amen UU\U x Own AOZ\Jmqmq<=92m 00+m00000000.0 00+m00000000.0 00+000000000.0 m m0-m00vvmm~h.0 00.M00000000.0 m01m00vvmmah.0 2 00+m0mmmmnv~.0 00+000000000.0 00+mommmmhv~.0 m: Invom AH.H.eom .~.H.mom ox\m=oe<-oxo mo+mooovmmov.o oo+moooooooo.o oo+moooooooo.o ox\.x one..qoz-ox. omom: AH.EE: .mlmmz >ma<=62m mmaetz ezHeomeem Dupe m>HeomEmm oooooo.o 0 E00 +mz mm mam ~z we max (2 mw\e a +2 mm\e S 2 wm\o a 4m: om\oHS m: om\ofla m wm\o 3 :memxm mHze 2H omemommzoo oszm mmmummmo Poem mo .Ho maa<> moo» oz 0 memuqmzH «00.0 no 0mm.0 04> no 00000 Umm\fim .xdfimo 0000.0 E<\m< 0000.0 0000.0 mMDZDZ 00¢: M0000.0 00000.0 02 4902400 00-000.0 00(M00.0 mumaumxmma.x 00-0v0.v 00(Mv~.v NPm\01mma.Dz 0.0000 0.0000 umm\:.Jm> 200 0000.0 0000.0 .0. «£240 0N00.0 0000.0 ammvwfio (2240 0000.0 0000.0 0hmvm<0 mu 0000.00 m000.00 0x0000\440 .00 0000.0 0000.0 0.600\>400 00000.0( 00000.01 9.maa\>000 000.0 000.n 93 do: .2 001000.0 00-000.0 AMNFHaxuv 200 0000.00 mmv0.00 Axvhovxadu .m N.000m0 0.00000 qudo .1 0.00000 206 WZOHEHDZOU QmZOHmm< 00¢ mom 0000000. 2420 mmmd mam: WZOHEU~H Hv.va -.omm vo.gm~ m0.~v~ m~.mHH mm.mmm zmq\umm-mm4 .H mm.o~0~ om.ooh~ ma.mpw~ vo.mvw~ mo.~mmfl mm.omv~ po.~¢m~ 0~.mmmfi oo.om~m am.mvnv nfi.nanm on.wmmfl 4\omm-mmq.o¢>H mow.H 00m.” 0mm.» 00¢.” pan.fl ”mo.” 0mm.o mam.o m0~.o nm0.° mo~.o Nflm.o mo 000.0 cam.“ mnm.~ mam.” vac.” mmm.H mon.~ mmv.a mnH.~ 000.. mmo.m mam.“ u<> mu momvn wwmvn momqm menu” mmmvm ommvm mwMVM women mmmvm mmmvm meme" mwmvm omm\9m .mcemu coo.mp ooo.m~ ooo.oH oooo.m coco.“ coon.H oofio.fl ooo~.H oooo.m mooo.q oooo.m coco.“ e<\m¢ ovmo.o~ mpmo.0 vane.m vaa>.n aoo..~ mnmm.fl smfld.~ mohm.o ema~.o pmvfl.o on~0.o oooo.~ oooo.o mmmzpz :04: H.momH n.0mmm n.~mdm H..0ov m.vvpm m.pm~p o.mvmm o.momm 5.000m o.mmpm «.mbha o.~omm ~.oowm umm\z.am> zom 0000.0 0000.0 boom.fl 0000.0 come.“ mmmm.~ mwom.H 0000.” -ww.~ mflmw.fl «Hmw.fi ommo.H mflmw.a .m. 42:40 ommw.a pmmo.a 0mm0.0 mmmw.d Hemm.H momm.d memo.fl mnmm.fi novm.fl vfivo.a mflvm.« mvmm.H mfivw.0 .x..u.\qsHmzmo oo+moooo.cnoHeHsom mo+mooo0.onqwsm ezmommm oo+moooo.oum\00 oooo.o oo.o o ooo.o oooHo.o ooooo.m z ampm 0000.0 00.0 0 000.0 00000.0 00000.0 fix amsm uo\u x own noz\g90m00m0> o 00-000000000.0 00-0mmmm0000.0 00-000000000.m 90x0 0°-mmmfifimfinw.m Ho-mm-mommm.0 mo-m0m0~mmm0.o 940009 “c-00000mm0m.m ~o-m0o0m~mmm.~ mo-moofippmm~.o mmmzAD. M0000.Hu 00000.“- 0N000.H: HHMH0.0- mahm0.H:_0hmm0.a- cavv0.ds 000v0.~u hmhv0.H: 00590.“- 00000.Hv Hmhvo.n- B.mao\>AQO hmo.q «00.v 0No.v 0m0.m 000.0 0m0.m 000.0 mv0.m 0m0.m 0mm.n 0fih.n 000.0 83 002 .z m0-mom.N 00um00.0 v0-mmm.~ voumho.v volmah.0 mo-mvn.H Monm~0.m 00-m0N.N m0-m0m.m motm0m.~ 00:m09.~ moxNOv.~ Nvmw.mu mvm0.ma Nvm0.nu ~vn0.na Nvmw.n« mvmm.na mvmw.ma Nvmw.m« Ammbdecv 2mm mvmw.MH memo.m~ Nvme.MH mvm0.m~ .x..uvxqfium2mo 00+M0000.0u0~9<¢ mmDPxHZ UHOEm 00+m0000.0noHEH00m n0+m000H.0uqum fizmummm 00+m0000.0nm\00 0000.0 00.0 0 000.0 000H0.0 00000.N 2 Juan 0000.0 00.0 D 000.0 00000.0 00000.“ m: Juan uuxo x own JOZ\J9Hmzm0 02m? m9ndm4<=92m 00+M00000000.0 00+m00000000.0 00+m00000000.0 w no:fl00v¢0mah.0 00*m00000000.0 m0-m00vv0nab.0 Z 00+m0N0mnhvm.0 00+m00000000.0 00+MONOMMFQN.0 m: Anvom .H.vaom AN.vaom OxmeOBfluwxo 00+m000NvHON.0 00+m00000000.0 00+M00000000.0 Uxxax OMD.AJO:»Ux. 003m: .vam: .Nvmm: >AJH90mmmm swam m>Hfiummmm 000000.0 n mOO +Nz 0m MDm N2 00 max :2 00x0 A +2 00x0 4 2 mm\w J +m: 00\0Hq m: 00\0ad m 00x0 4 2m9m>m WHIP 2H DmmeHmZOU OZHmm mmHUmmmo 900m mo .dm .deom .mo mob zm>HU mDJ<> mDOm OZ 0 COCCCCCCOCOCCCCCCCOCCCfiCOCCIdCCCCGCCCCCCCCOC‘CCCCOCOCCICtCOCCCCCCCCCCCOO.CCCCCCOOCCCCCCCCCICCCCCCCCCCCCCCCCCCOC...CCCCCC N mZON mOm mmqbomxum mmbmmmmm D24 Oufidm u 000.0 000.0 000.0 v00.a 0-.H 000.0 000.0 000.0 0a~.0 v0~.0 000.0 000.0 LU 000.0 000.H 000.0 000.0 000.0 000.H 000.0 000.H 00H.N 000.0 000.0 000.0 04> mo 0v~00 00000 0va00 00000 0qa00 0va00 0v~00 0¢H00 00000 00000 0va00 00000 umm\00 .04900 000.00 000.0N 000.0H 0000.0 0000.0 0000.0 0000.H 0000.0 0000.0 0000.v 0000.v 0000.0 dem< 0000.0 a00a.0 «000.0 00H0.~ 0000.0 0000.0 NO0H.H 0000.0 00H0.0 v00H.0 00H0.0 0000.0 0000.0 mmmzbz :04: 00000.0 00000.0 0H000.0 00000.0 00000.0 00000.0 00000.0 00000.0 00000.0 00000.0 00000.0 00000.0 00000.0 02 0002400 ~0IM0v.0 00:000.0 dotm0o.u H01000.fl H01000.H Holm00.H 00:000.a H0u000.a “0-0v0.u 00:000.a ~0-m00.d ~01m00.d 00‘000.H mwmeumxmma.¥ 00um0H.a 00:000.H 001m00.N 00:000.N 00-000.0 00:000.0 00-000.0 00n0~0.0 00:000.0 00nfl00.0 00n000.0 00vW¢0.0 00-000.0 «memnmmq.02 0.0000 0.0000 0.0N0v 0.0H00 0.0000 0.00v0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0v00 0.0000 omm\z.qm> 200 0000.0 0000.H 0000.0 00H0.H 0000.H 000H.a 0HOH.H 000H.H 0000.0 NO0H.a ~00~.H v00a.H NO0H.H Amy 42240 0000.0 0000.0 c000.H 0000.0 0000.0 0000.H 0000.0 ~000.H 0000.0 0000.0 0000.0 0v00.~ 0000.0 .mmvmdo 4:240 0000.0 NHO~.A 0000.0 000~.H 0000.0 0000.H «000.0 0000.0 0000.« 0000.“ 0000.H 0000.0 «000.0 .mm_mHDOm 0000.0 00.0 0 000.0 000H0.0 0000.0 00.0 0 000.0 00000.0 UUxU x 000 AOZxJ9HmZmo 0209 wk<fim >00n k0 U<> mo Umm\bm .mdfimv 0<\m¢ 0000.0 KMQZDZ =U 200 0000.0 .00 42240 H000.~ Ax..UV\A 0 210 . .anH ”.momn o.mafim v,movv m.m~mm m.mfifiw m.hm~m o.mmmo c.5vmw “.mv . . wmmwmm “mow." wmmm.H mqmm H haw» d mmaw a omnm H ovmm H mvam.~ amen.” muommw MmMMww “mmwmw ommxz.qm> 20m comm.H mmew.~ mmow H «mom.H mumm.~ ommw.~ “www.fl ”umm.~ ommm.n mmwm.« mmwo.~ wwwo.fl ammm.~ .m. (zzce mann.a 50mm." on~.H afimm H m~n«.~ nvnm.a gnaw.“ vn.~.~ hnv~.~ omu~.~ ~m.~.~ oom~.~ yawn.” .mm.mqn. uno.q m~o.. mam.n mmm.n ohm.m “ha.” wmo.v ”hm.m H- gMqu\>qo. o: .: ooo.Hu ooooo.au oc bno.v hmo.v Hmo.v fimo.v bmc.v n01m~H.M MOINMN.M MOIMmN.n no:moo.m MOuMNN.M Ammfiuaxo. ZNQ no,mom.~ no-mom.H no-mm~.~ ommv.m~ mmmv.- amm..- .x..o.\gfiumzmo oo+moooo.ouoHHnbom mo+wooon.onqmbm Ezmommm oo+mooo . n o oooo.o 09.0 0 ooo.o oooHo.o ooooo.m z o c:W\0o oooo.o 00.0 0 000.0 ooomm.o 00°00." at “man UUxo x Own qoqu9Hmzmo mzmb madem >mqmq<=92m oo+moooooooo.o OO+moooooooo.o 00+moooooooo.o u mo‘mmovvaHb.o oo+moooooooo.o moymmovvmmH>.o z oo+Mo~ammsvm.c oo+moooooooo.o oo+mo~mmmnv~.o m: .Hvom .a.H.mom .N.H.mom oxmeOEma¢=Hzm mmDBXHZ BZHFUmmmm JmDm m>HPUmmmm oooooo.o n mOo +NZ am mam N2 «h max :2 moxb a +2 mmxb a z mm\w A +@: om\0Hq m: omxoaa m wmxm a Sufim>m MHz? 2H QMKMQHmZOU UZHmm wwwuwamo 90mm mo .HU MDJ¢> MQOm OZ 0 CCCC¢CC§CCOOCIC‘CJCCQOCC.OOCCCCOCCCCCCOC.‘OOCliCCCIiCO.CC¢ICC...Oil.C.CCOIOOCCCCCOOCCOOOIlCCCCCCCCOOiCC’CICCCOCC‘C’CCO.C m WZON mom mmabomrom mmbmmmmm QZ< OHEEHMOUMN> AZMNOmh UZHBKHZWV mom moz<fim «mm. Dmmmonzoo mum: ZUH23 WEUDQOmA 000000.0 vmm000.o 0vn000.0 hmhoH0.o Nuwmh0.0 N00000.o 000000.o avhN00.o mmoooo.o 0m50~0.0 ommme.0 vam«o.o bb0.n 0~.mh me.>0~v v00.0 hvo.m wmeom mmma.v mom«.0 vwwmo.o 000000.0 5vm000.0 nm0000.0 00H~00.o vmmmmm.o H00000.0 000000.0 ommmoo.0 50H000.0 0HH~00.0 nmmmmm.o ombmoo.o 0H0.v mn.vmo mm.mm0~ mmh.0 chm." mmwwm 0000.“ 0000.0 mmwwo.0 ame .3955 mmmz<=u 20H949w o H maoz +~2 - qazonbmooco oooooo.o mz vmwmoo.o *2 mvmooo.o z vmmoHo.o +m: m00050.0 m: «00000.0 m monaom no o<> mo ummka .xaawu demd oooo.o ammzpz =uu ppm.d can." she.” FHM.H fimo.a Hmm.o mam.o mh~.o nnfl.o mo~.o NHm.o mm mam.” «hm.fl mmm.~ vmq.~ awn.” mom.H wm¢.~ mmfi.~ wwo.v nmo.m mam.” 04> mu ovamm e¢fimm wcflmm magma wVHmm wvamm wvflmm wvflmm wvflmm wvamm wVHmm umm\em .xaemo ooo.m~ ooo.o~ ocoo.m oooo.~ ooo~.~ ooHo.~ ooo~.H oooo.~ oooo.~ mooo.m coca.“ eaxm< “vac.“ mmmo.m mamb.m moov.~ mhmm.d hmHH.H warm.o ~mm~.o mmvfl.o mmafi.o coco.“ oooo.o mmmzbz mu zom powm.fi poem.a memo.“ mmmm.fi ”www.fl Nwmw.H owmo.~ ammo.” mmwo.a ~mww.a mmwm.~ .m. 4:240 5000.fi 00v~.H wmvm.u wmvm.a mmvN.H 00vm.~ vam.fi mmvm.~ mmvm.~ vmvN.H vmvm.a Hmvm.d vmvm.~ .x..o_\qhumzmo 00+m0000.0uOHFHDOM m0+m000a.0uqmam szummm 00+m0000.0um\00 0000.0 00.0 0 000.0 000H0.0 00000.n 2 Juan 0000.0 00.0 0 000.0 00000.0 00000.H m: dem UU\O x vmo JOZ\A¢U Amsoz mmmv (Jazmom J 20m Hmww.H owmm.H amww.H wmwo.H vme.H Hmwe.H mmmo.H Hmmm.H «mmm.H mmmm.H mmwm.H HHHH.H Hmwm.H .am..m. <=zgo. pHo.H “Ho.H HHo.H hmo.H “mo H wHo H mmo v Hmo.v Hno.v HHo.v HHQ.» Hmo.~ Hmo.v ezqmmnqflh -m .n mo-mmm.~ Ho-mHo.~ Ho-maH.H no-mmH.H mo-mmm.~ Ho-mav.m Ho-mmo.m mo-mvm.m Ho-mHm.m Ho-mmo.m no-mmm. - . ”MHH+MH HHmm.HH Hmmm.HH Hamm.HH HHmm.HH Hmam.HH Hmmm.HH Hmam.HH HHmH.HH Hmam.HH HHmm.HH Hmam.HHm ”mmwmem ”wwwmwmwmuzma H.Hm- m.m- H.HHH o.~HHH m.mn~n v.memm H.oumw H.5mmm ~.mOHm H.Hnmm ~.mmmm o.HnHH o.ooooH o\afiHmZmD 00+m0000.0uOHEHDOM m0+m000H.0ndm3m FZNUZNQ 00+m0000.ouk\00 oooo.o oo.o u ooo.o o00H0.o ooooo.~ z amen oooo.o oo.o o ooo.o ooomm.o ooooo.H m: qmam UUxU x DMD JOEdeU .mEOZ mmm. fidbzxom Q(UHZMIU NEHmZNQ mzwfi m&(fim hmdézfizm 20H90<¢h 93 «Hmm n.vH u umo ZOHm2HFUmmmm +m2 00 max +0: 00x0a4 00+mo0000000.0 no-0009vmmah.0 00+mommnm5vm.0 .N.H.mom 00+mo0000000.0 .m0mmx qm0m m>HFUmmmm m2 0h mbm m: 00x0dq 2w90>m WHIB 2H ommeHmZOU OZHMQ MMHUmmmo 90mm m0 .(m .deOm .mO mom 2m>HD m044> 000» 02 quxbbm 2H >mqAJ4ZEZN 000000.0 u k00 :2 mmxb A m mwxw A ...“.OCOlfiiliiiCiC‘flthO.I...‘CCC£CII.CCCCC.CQICCOIiifiiICOCIOICiCGO0C0.0‘CfifiifilCCOIC‘COOCCCOCC....CCCCCCC.ICCCCICCCCCICC 0 214 ooooHo.o oooooo.o oooooo.o oooomm.o oooooo.o vaHoo.o oooooo.c oooooo.o mmmmmm.o oooooo.o pmo.o vo.nma o~.~vm nom.H oHo.H Huan 000.mh 0000.0 0000.0 UUxO >9HWZMD . 2 oo.o o ooo.o oooHo.o wmmmm.m m: oo.o o coo.o oooaa.o_ mq<2ezm oneo o Ho-mH~HmHme.e mo-momHHHHHm.H Ho-mmoHoHon.H eme Ho-memHHHmm.w Ho-movaHHno.H oo-m>m~HaHHH.m emo.v wmo.v «no. . ow.HHH HH.HHH Mw.Hom Hm.HHH Hp.wHo wo.mwm “Mo H Hma H HHo.H HHo.H HHo.H HHo.H .xH2.az Ho: ~mm.H mmm. mmm wH.mHm wH.~m . .mm~ Hm owH HH.HH «H.Ho HH.~HH - . Hoo.H H one H mHm.H H HH amp mm Hmm Hm.vva H . . mmHHQmm mm; H HonH ”WWW” Hem.H me.H mmm M omm.o va.o mHH.o HMHmwmm memmam MMWHWH onmm-mmq.uc>H . . . . o ooo mm ooo.oH ~wmww.m mman. mman MMMHM MMMmM «HH.~ mwo.v mmo.m aa~.H oas my ooo n ooo~.H ooHo.H ooom. ~mHmH. HmHmH mmHmH monH omm\am .mHPUmmmm amDh m>HBUmmmm 000000.0 n moo +~z 00 max N2 00 max :2 wmxb a +m: 00x0"; 0: 00x0flq m 00x0 4 Imem>m max? 2H DMKQQHWZOU wszm mmHUmmmo 900m :0 .4m .E4m0m .LO 00m zm>HU m004> 000m 02 0 m WZON mom mmqbomxvm mmbmmmmm 024 OHF4K 4mm4 moo WB4JDUJ4U .cc.uca«ccccccccccocc «....ocaocccccc.00ctacctc.ocacv¢¢ccc¢cccacacouocccc.oc#ccccocacti-cocoa...cccccccccacui¢ccocoacccc. MZOHEHGZOU GmZUHmm4 004 mom m000000. 24:9 cocooo.o Homooo.o pm.mHa mm.vmm oH.oom Hm.a~a mH.mHm pH.wmm ppm.H wmm.H HHH.H cam.H H»H.H mHm.H HHHmH pmbmH HHHHH ooo.m~ ooo.oH oooo.m wmmo.p vaHo.m Hmmn.n o OHHH H.0HHH m.o-~ Howw H Hmoo.H wam.H «Hmm H ”HH~.H mHH~.H Hno v Hmo.H HHo.H mo-mmm.n Ho-m~o.~ Ho-mH~.v “mam HH Hmmm.HH Hmam.HH wen «.mHm m.HmvH onv . Ham mHHH ”mommmH wHoo o HHHo.o v a r ame o mpm HHH Hm 00.00h ah.wvm me.H vmv.a finbmn 0000.N Naov.m H.~m~m hwwo.u 0~MN.H Hmo.v Ho-moH.H ~mmm.~a «.mvmm mFZ4DHXO J4BOP 2H PZ4DHXO m0 024 meDm J4EOE 2H dem LC 20H904xh 9:0Hm3 mmma mum: WZOHEU4ML 0008 mmOIZ 800 Dmmmonmzoo m NZ 000000.0 +2 ommmoo.o +2 mhwm00.0 2 mHmmmm.0 mo.mmo pm.mHm Ho.oom mm.owH Hm.bp om.Hw po.Hmp vo.mmn mo.wmm mH.Hv~H NH.mmHm mw.mvmm nmo.H Nam.o mHm.o mp~.o HHH.o ooH.o 0mm.H mom.H wmv.H mMH.~ mwo.v mmo.m nman uman anmH HHHmH anmH HHHmH ooo~.H ooHo.H ooom.H oooo.~ Hooo.v oooo.m ommm.H anH.H nmnm.o momm.o mmvH.o mmHH.o m wean H.0mvv m.mh0m m.mn~m m.mmmm v.mmmm hmmeH mwom.H mmow.H mmmo.~ «mom.H mmmo.H mHmw H mHmm.H Hmmm.H mmmm.H vmmm.H vmmm.H Hmo v Hmo.q Hmo.v Hmo.v Hmo.v Hmo.v ”MmMMMHN no:mmv.m mo-mmo.m mo-mvo.m no-mmm.m no-mwm.m a.momm ”mam HH Hmmm.HH Hmam.HH Hmmm.HH Hmmm.HH ommv mpmm m wmmm o.vohm o.Hmmm H.0mmm ohH~.o comm. pawn mmom mhmm comm momm.v Man.o opmp.o mHHm.o mmmm.o woma.o eme van N amon.H nmbo.H cpHo.H moHo.H <2 mass eme eme stm eme ewa 2H2 oHosm oo.moooo.ouoHe<2 muzmq<>Haom +NZ 2 m: hr.mhv mv.0mh -0.0 mmN.H hmhma 0000.“ 0000." n.0m0v v000.H Hmmm.d Hm0.v m0:mma.m m0:M00.n anam.HH «.Nnvn mmmm mhmv.0 bmm0.N E40m26 MPOZ :2 om: mm: :UHIZ mBUDQOmm A4ZOHBHQQ40 000000.0 w: mZOHEU4dm mm4z H00000.0 N2 v0m000.0 m MZOHEU4ML mac: quxomm:bmq .H q\0mm:me.04>H m0 04> mu 0mmxbm .m4PmU P4xm4 0000.0 «@0202 :04: 0.0mmm 0mm\z.qm> 20m mm00.~ .m. (2240 vmm~.n «30.0.xq40 .mo Hmo.v 93 002 .2 .mmBHJ\00 2mg Hmam.H~ Ax..00\440 .m 0.0000H 0x440 .z mmnm x own .9 000.H 294 .m 0000.H m\Um zmm24=0 0 00+0000~.0u000m BZMUmMQ 00+00000.0nm\00 1216 hm0.v hm.vhw 00.000 N00.H «Ho.u vvaH 000.mh mmav.0H mvmwo.0 N0:mwo.a h0:mom.fi h.vmw m0:flmm.v vwhh.0~ v.NNN: ban H000.0 hm.wwv0 Eme 000.0 000.0 000.0 0H0.v 0H0.0 000.0 00.000 Hv.000 00.H~0 n~.vmm 00.H00 H0 000 00.H00 H0.H00 00.000 00.HHO ~0.Hpm 00.000 000.H omm.H 000.H th.H 000.H H00.o 00m.H 00m.H 0H0.H vmv.H mmm.H HOH.H vvmmH HHmHH cvmmH vaHH vvmmH evaH 000.00 000.0H 0000.0 0000.~ 000~.H 00H0 H 0000.0 0000.0 m000.n H000.~ 0000.H 00HH.H 00H00.0 00H00.0 00000.0 00000.0 0H000.0 0H000.0 ~0-m~0.H molm0m.~ ~0-mmm.n ~0-mnm.m ~0-mmm.0 ~0-mHm.0 00-0H0.m 00-mmm.m 00-000.» 00:mmH.H 00-mmv.H 00-2m0.H 0.0H0 H.0mmH 0.000H m.000~ 0.~00~ 0.0HNH 0000.H 0000.H 0000.H 0000.H 0000.H ~000.H 0000.H 0000.H 0m00.H 0000.H 0000.H 0000.H Hom~.H momm.H 00H~.H 000~.H HHm~.H ~HH~.H vom~.H 00H~.H 00HN.H momm.H HHm~.H 0HH~.H 0000 H 0000 H 0000.H 0000.H 0000.H 0000 H 00000.H- 00000.H- 00000.H- 00000.H: 00000.H- 00000.H- 000.0 000.0 0H0.0 0H0.0 0H0.0 000.0 00-2H0.H 00-000.n 00-2m0.0 m0-20~.~ H0-m~m.0 H0-200.0 0000.0H 0000.0H 0000.0H 0000.0H 0000.0H 0000.0H m 00. 0.H0~ 0.0mm 0.000H 0.0mmm 0.0000 mmmo.o ~00 . mmn mmvH mnmm 000m 0m.H~HH 0H00 0 H~H0.0 0000.0 00H~.0 00H0.0 eHx 00 000 ~mm.H0 000.0H 0H0m.0 0000.0 00w +me eme stm pme eme oooomwooo ow>aHmzmo 00+20000.0noH0<2 mmoetz uHoem 0000.0 mm.0 0 000.0 000H0.0 uo\0 2 M00 0 000.0 00000.0 >eHmzmo 0220 202\2<0 22002 mmm. mfitfim >ma4zfizm ZOHEU4¢L 93 onmz< mxm qumso onsHmOmzoo onmquHzom 02H20mm< 20242200220 pm 00.200000000.0 H0-m000000H0.0 00+20~00000~.0 .Hvom 00.00000HH00.0 00002 hmo.v hm0.v 00.0aN 00.naa mm.v00 hm.m00 mam.0 mhm.0 mmv.m mm“.m vvmma vvaH 000~.H 0000.N momm.0 000~.0 mmv00.0 0Nv00.0 N0-M0q.0 N0:me.0 00:Mmm.~ 00:fl00.a m.m00m 0.005n 0000.m mmmw.~ vm00.H vm00.H Namm.u NHmN.H mhmm.a mvvm.H v000.H 0H00.a H0000.H: v0000.n: hm0.v bm0.v n0:3m0.m N0:m00.n 0000.0H vmhb.0a 0.H0vv 0.0vmv Hmmm mmmv «000.0 mm~m.0 mm0m.~ mmho.n Pme mem m n mZON 00+moooooooo.o 00.000000000.0 oo+moooooooo.o .H.H.mom 00+moooooooo.0 2H0002 >m0.v wm.mm hh.HHhH nMH.0 000.v vvmna mmmm.n 0mv~.0 Hmvwm.0 N0 mam.0 00:mmm.~ n.0mmm Nmm0.H mmww.a NHMN.~ Hmvm.m MH00.H m0000.n: bm0.v N0:WNH.H v0hb.0~ 0.0000 hwnv vm00.0 00~0.H anm 00+00000.0u0~94¢ muzmq4>Hbom 0H0.0 0H0.H “no.0 .2H2003 002 00.00 00.H0H 200\uwm-000 .H 0H.H~HH H0.0Hm 0\umm-002.0<>H 00H.0 HH0.0 00 Hmo.m 00~.H 04> 00 HHmHH HHmHH ummHFm .24000 H000.m 0000.H e<\m< 0HHH.0 0000.H 0000.0 200202 2022 HH000.0 0H000.0 HHH00.0 02 0002420 H0-m~0.0 00-000.0 ~0-mqm.0 2000-0H200.2 00-000.H 00-000.H 00-0H0.H ~00\0-000.02 0.~00H 0.0HHH 0.0000 ummH2.0m> 200 «~00.H 0000.H HHH0.H .0. 42240 0000.H 0000.H 0000.H 220.020 42220 HHHH.H HHHH.H HHH~.H 200.020 00 0000.H 0HH~.H ~00~.H 220200.040 .00 HH00.H 0000.H HH00.H 0.900\>00. 00000.H- 00000.H- 00000.H- e2000\>000 0H0.0 “no.0 HHO.H 03 002 .2 ~0-m~H.H H0-mHH.H 00-0HH.H .mmeHgHo. 200 0000.0H 0000.0H 0000.0H .20200H020 .0 0.0000 0.000H 0.0000 oxqau .2 FHHH ~0~H mmHH 2 000 .0 0000.0 H00H.0 000.H 204 .0 00H0.H 0000.~ 0000.H 0\00 Fme B40x=6 mmmS4ZU 0 H0.2000H.0u0000 0220200 00.00000.0u2\00 00000.0 2 0002 00000.H m: 0000 4002200 0<0H2220 mJ4=sz m 2 w: U¥\m2094:0¥0 ox\.x 0m002402:030 >mq4zfizm 217' v0-m00.0 m0:mwm.~ no mam.v m0:mmw.0 m0:mmo.0 m0-mwo.H N010ma.a N0-mma.a m0:mvn.h N0:mn~.a HxMEHq\U. sz 00-200.0 00-2Hm.H 00-200.H . . . 0 .0H 0099.0H 0099.0H 0099.0H 0099,0H 0099.0H 22.20.2220 .0 0029.0H 0099.0H 0009.0H Howwon “wwwva ”wwmmmH MMMMHMH wwHwov 9.000. 0.0000 0.9900 H.~00H 0.0002 0H2<0 .2 0.0HN... n.~.ml m.hON N. QwNm mmmv & ONO .h. 0 9H0H 0HHH 0NHH 0HHH . 9HH 000 H00 H09 HmoH. 00mm. 09 H. . . HH00.0 0000.0 0900.0 000.H 292 2 H000.0 0000.0 0H00.0 H~H0.0 0000 0 99HH 0 09H» 0 0909 0 0000 w 0 H0 H 00H0.H Hmmo.« 0000.H 2202 . . . . . . . 9090. 9 . 0H.0900 00 HHHH 00 090 002 H0 090 HH 0002 H HHmH N HH0H H 0 H22 9H22 920229 2222220 0.me 838 833 .0me BHXM 0.me Puxm mem 9.me a. . a +0 0. "2:00 00+20000.0u>9H0220 00.20000.0uoH9<2 22092H2 0H09m 00+20000.0H0H9<2 20220<>H002 H0+2000H 0-0202 92202MM00W0~ 0N0 00202 0000.0 00.0 0 000.0 000H0.0 00000.H 22 0202 0000.0 00.0 0 000.0 00000.0 20 00Ho 2 020 202\04u .2902 2200 2002202 09Hw2wo AIME. mfi4fim >mq4zg 200.905..” 0.3 4Hmm b.va I 0&0 20Hm24mxm OZHmDD ZOHPHMOQZOU Zmnomu OZHrSmm4 mUZSLMOmem meUOm J4omfimmomzfi m u mZON m 02 2 22 mon920000<0 02H9222020 920202229 2H 02220H2200 22H0200 0 2020-920H220 00-2HH009H9H.0 "222 09 9022H 202 22202 ~90mm.0 0H 20220.._2229\9..22202n02 2202 229 202 220220. 92220222 >9Hmo0mH> 0 Ho-200002000.0 00-2H09NHH00.H 90-2HH00990~.~ 9me H0-2H0000HH0.0 «0-20HHHmmH9.9 00-20000H090.H 920229 H0,2omHHHH00.0 No-2mH000000.0 00-20H00HHH0.H 2222220 .m OmQ:Umm\&mqv Amc.Pm\Umm:mma0 22 2 02 20H929m monfi4mFZMUZOU EDHKquHDOm 20mm 0094000440 mmmfimmmomm Fmomm24xfi szOxm 0 H 2220022 02H92H2m. 202 002290 .20. m92<0H20 02909 2H 92<0H20 20 022 02202 02909 2H 0202 20 20H90<22 920H23 2902 0.2035200 0220824 .3 :2 +22 2 . 2 202 . mmwmw0.0 000000.0 000000 0 000000.0 00000M0m0mwm002229 0222 2222 mon90<22 2202 20023 902 02220H0200 2223 20H23 09000022 0220H9H0020 00000N.w 0000H0 0 0000H0.0 0000H0.0 0000H0.0 0000Hm.0 000000.0 000000.0 000000.0 000000.0 000000.0 H00000 0 0H0000.0 +~2 000000.0 000000 0 000000 0 000000.0 000000.0 000000.0 000000.0 000000.0 000000.0 H~H000.0 HOH000.0 900000.0 000000.0 «2 000000.0 wmwmmm.0 000000.0 000000.0 000000.0 000000.” mwmm00.0 000000.0 000000.0 000000.0 000000.0 000000.0 9H0000.0 +2 0 000000.0 000000.0 000000.0 000000.0 0000M0.0 00H000.0 022000.0 090000.0 9m0000.0 HH0000.0 HH9000.0 2 0 0 000000 0 000020.0 000000.0 000000.0 000000.0 000020.0 22 00000 . . 20H90<22 0022 330M.” wwwwmwfi wwwwooé 0000000 0000000 0000000 00000 m M00000 0 000000.0 0000mm.” WHHHoo 0 HO0H00.0 HHHH00.0 HHvHow.m mwm000.0 000000 0 000000.0 000000 0 000000.0 ~00000.0 .02 m””0000 000000.0 000000.0 0m0000.0 000000.0 000000.0 000000.0 000w00.0 00HH00.0 HHHH00.0 HHHH00.0 0HOH00.0 0HHH00.0 «2 000 0 000000.0 000000.0 mamm00.0 000000.0 000000.0 H00000.0 0200m0.0 000000.0 000000.0 000000.0 000000.0 000000.0 .2 00 0 022000.0 000200.0 002000.0 HHmmmM.0 00H000.0 00H000.0 H00000.0 «00000.0 HH~000.0 2 0 000000 0 H00000.0 000000.0 900000.0 90000m.0 22 mZOHEU4dL 0402 i218 mm.mrw mm.hho Noo.H Haw.” vmmna ooo.mh Hume.oa N.vnm Howw.a momm.a hmo.v mZOHEHDZOU szommm< 44¢ mom moooooo. 2dr? mmmq mmmz MZOHEUH LU 04> ho ummxfim .m49w0 9 20m .m. 42:40 AxVAU.\J A' + e' electron excitation A' ... e‘ ———> A + e' electron de-excitation A + B -———> A' + B neutral excitation A + hv ———> A' radiative excitation A' ——9 A + hv radiative de-excitation 220 e-s mm e—> Electron-Neutral l .—->.---> G—> Ion—Neutral mom» 6» Neutral-Neutral G—> m» L Radiation Figure D.1 Collisions Disassociation: A2 + B ———> 2A + B collisional disassociation 221 Ionization: A +e' —->A+ + 2e' electron ionization A+ + B —-> A" + B+ + e' ion ionizatation A + B —-> A+ + B + e' neutral ionization A‘+hv—)A++e' Capture: A+B’ —-—>AB' Recombination: A*+e' —->A+hv Transfer: A*+B ——->A+B+ A'+B —>A+B' Superelastic: A'+B —-—>A+B+lskin radiative ionization ion - neutral capture radiative recombination charge transfer excitation transfer superelastic collision 222 Elastic: A + B ——) A + B elastic collision The last type of collision was elastic, in which no energy was transferred. In reality, a small amount of energy would have been transferred. The equation for that energy would be (M 1 was the mass of particle l and M2 was the mass of particle 2) [Cherrington, l979k 2M1M2 D.l (M1+M2)2 ‘Eelastic = D.1.2 Cross-Section. As shown in Figure DI, the cross-section of the particles were important in determining whether they would collide and react. Each type of collision (inelastic or elastic) had its own collisional cross-section. This collisional cross- section was important for calculating collisional rates, which typically could be done using a distribution function. These calculations were described in subsection D.3.4. D.1.3 Coulomb Forces. As previously mentioned, some particles (such as ions) could exhibit coulomb forces upon one another. These forces were important when dealing with electrons and ions. For example, electrons repelled one another but attracted positive ions. The kinetic energy of a particle could be represented as [Reed, 1973]: 223 E_ mvc 13.2 Defining "b" as the impact parameter (or separation distance) of two particles approaching one another, the centrifugal potential could be defined as [Reed, 1973]: _mv(2, b2 D'3 Vc 2r2 The effective potential was the combination of the coulomb potential (V) and the centrifugal potential [Reed, 1973]. Therefore, the kinetic energy could be rewritten as [Reed, 1973]: 2 D5 When this energy was equal to the effective potential, one obtained a solution for a minimum radius (or closest approach). If this minimum value was within the confines of the collisional cross-section, one obtained a collision; otherwise, one obtained a 224 reflection. Examples of potential plots were illustrated in Figure D.2 for both repulsion and attraction potentials. D.2 Charged fiarficle Motion in BM. Field. D.2.l EM. Field (TM012). The following were the electric and magnetic field equations for the microwave resonance cavity used in my thesis. The TM 012 mode was used. The derivations of these equations were described within my thesis [Haraburda, 1990]. Classical Potential Plats (Neutral —- Neutral Collision) 5— '75 "3 a; 25 fire a NO 1... . D. : Hard—Sphoru & °‘3 fl 1 La. 6.12 ”1? rut—Square Well —2;! n.3,. #-"6"'é"'jl Radius (Angstroms) Figure D.2 Classical Potential Plot 225 _ +. E2 = Eme(kc1’) COS(m00e(_JBPz) D.6 - 2 E _JBpRoEl m Pmnr Pmnr Pmr GB Z) D7 r-T— —Jm — - Jm+1 cos(m90e p Pmn r R0 R0 R E9 = 2 o sin(m9i eOsz) PmnRo -jmngmElJm[Pmnr] . D'9 Hr = 2 sin(m9i eOsz) PmnRo ' 2 . D.10 HG = J08R0E1[2Jm[Pmnr]- PmnrJm+l[Pmnr]]cos(meoc(’sz) 2 R R R Pmn r o o 0 Hz =0 D.11 D.2.2 Equations of Motion. In a collisionless environment, the Lorentz force equation was used to characterize particle motion in an E.M. field Dance], 1966]. F:q(€+?x§) D.12 226 The force could be related to mass (m) and acceleration (a) of a particle [Halliday, 1978]. Fzma D.13 The fundamental definition of acceleration was the time change of velocity. Therefore, the Lorentz force equation could be written as [Jancel, '1966]: - _ _ _ D.14 m[9dTv]=q(E+va) Unfortunately, one needed an equation that would account for collisions. The following was the Langevin equation, which modified the Lorentz equation using the collisional frequency (vm): ICELangevin = FLorentz -vmm v D.15 Recalling the phasor transformation of the derivative operator (d/dt = jar), the Langevin equation could be written as: (vm +iw)m?=qfi+q?xfi ”‘6 227 Particle motion in a collisional plasma discharge with an electromagnetic field was characterized using this Langevin equation. D.2.3 Power Absorption. Power was transported to the plasma from the electromagnetic field. Power was defined as the time derivative of the energy. Of importance was the time average power. This average power was proportional to the cyclic time integral of the force equation. This integral should be done over a cyclic interval [Cherrington, 1979]. 2%) _ D.17 If Vrn = 0 (Lorentz force), the above integral would become zero. Thus, power was transferred only in the presence of collisions. Physically, power was transported through collisions because of the random component of velocity (C). The velocity field could be expressed as the sum of the average velocity and the random velocity [Cherrington, 1979]. - "' D.18 Note that the average of the random velocity was zero. The energy of a particle could be written as [Cherrington, 197913 228 _2 _ __ _ _ mV _mvg+mc(2)+m(V0C) D.19 2 2 2 1 E: Although collisions did not alter the average velocity (or kinetic energy), collisions did alter the random component. This change in the random component was non-zero in the integral. Thus, power was transferred. The equation describing this power transfer was written as [Cherrington, 1979]: -1, D20 D.3 Distribution Function. Energy, mass, and momentum transport phenomena were involved in this system. To predict those transport processes, one should solve the conservation equations for each. The following described the important theoretical areas involved in developing an accurate fluid flow model. D.3.l Statistical Mechanics. Because of the high temperature region and the low number & high cost of accurate diagnostic equipment to measure thermodynamic Properties of plasmas, theoretical methods must be used to predict them. Thermodynamic Properties could be obtained through the use of statistical mechanics. The following 229 subsections discussed the use of partition functions and the pertinent mathematical formulas for calculating several thermodynamic properties. D.3.l.l Partition Functions. Using statistical mechanics required the use of partition functions (qn), which were functions of temperature and volume. Canonical partition functions were defined as the sum of energy functions [McQuarrie, 1976]. Q(N,V,T)= )3 el'EiB) 13.21 ' 1' Whereas, E: was defined as the sum of energy particles. Ej = 2 8i,j D22 1 These could easily be rewritten into partitions. ...-aria i Therefore, the canonical partition function could be written as the product of partitions. 230 i . D.24 QN = H Clr QN = [q(V,T)]N distinguishable particles D.25 N indistinguishable particles D.26 V,T QN = ———q 12, ) For helium, the molecular partition function was only a product of its translational (Citrans) and electronic (Clelect) partition functions. Using quantum mechanics in rectangular coordinates (x, y, z), the energy states of a particle could be written as: h 2 2 2 a 13.27 8:8 2[nx+ny+n‘z‘] an,ny,nz=l,2,3,~- m a The translation partition could be written as: oo oo 00 (—Ba) distinguishable particles D.28 Qtrans = 2 2 X e nx=l ny=l nz=l indistinguishable particles D.29 n: 3 ‘ltrans =[ 020 36%)] 231 Statistically, over an infinite number of points, this summation could be rewritten as an integral over the n points. l [437,2 “2] ) D.30 co 2 [e 8ma dn 0 Qtrans = 3 D31 2Pka A (ltrans: T V The electronic partition function could be written as: . D.32 (lelect = 2 me,i 6GB 81) 1 By fixing the ground state energy (8i) to be zero, we obtained the following expression for the electronic partition function: - D.33 Clelect = 2 men 9H3 A8161) i 232 For nitrogen, the molecular partition function became more complicated in that the product included the vibrational (Civib) and the rotational (QrOt) partition functions. The vibrational and rotational partition functions were defined as: qVib =e(-B h V)? e(.Bh V n)dn D34 0 1 = for >>hv Bhv )6 and, qmt =0]? e('AB J(J+1))d{J(J +1» D.35 0 1 :— f >>A AB 0. )6 with "A" being defined as the rotational constant with a value of 2.001 cm'1 for N2. Values obtained through these calculations assumed ideal gas conditions and local thermodynamic equilibrium. D.3.1.2 Chemical Equilibrium Reactions. Thermodynamic values were calculated assuming equilibrium conditions. The following ionization reactions were used in these calculations: 233 He <—KL-> He” + e‘ He+ 649—) He++ + e' The equilibrium constants could be related to the mole fractions of the species as such: X X. — D.36 K 1 = P He+ C XHe X X _ D.37 K2 = P HC++ C X + He The values of these constants could be obtained using statistical mechanics and exPerimental values for ionization energies (Aei) [Bromberg, 1980]. D38 q .. K1: e('BA81)L_ QHe D.39 ( 5A5 )qH ++ q — K2 = e - 2 —C_e__ qu-r- D.3.1.3 Species Mole Fraction. As previously mentioned, this provided two non-linear algebraic equations and four unknowns - the mole fractions of 234 the species. Assuming particle conservation and electrical neutrality, one could derive two additional equations. 2 Xi = l ' D.40 i X - - X — = . D.3.1.4 Average Molecular Weight. An important parameter in any chemical reaction was that of molecular weight. With the known molecular weights of the species and the calculated mole fractions, one could calculate the average molecular weight of the plasma using the following formula: (M) = .2 X, M: 1142 l D.3.1.5 Compressibility Factor. Using the ideal gas relationships and compressible fluids, the compressibility of the plasma became important. The following relationship defined the compressibility factor (Z) for calculating additional thermodynamic properties. Mo 13.43 235 Here, Mo was defined as the molecular weight of the fluid at 273.15 Kelvin. D.3.1.6 Plasma Density. The density of the plasma should be known to be able to develop model equations. The functional relationship between density, pressure, temperature, and compressibility was: _ p0 PTO D.4-4 " ZPOT Here, p0, To, and Po were the values at standard state conditions of 273.15 Kelvin and 1 ATM. D.3.1.7 Energy / Enthalpy. It was essential to know the energy of the plasma so that one could do heat transfer modelling. The energy level of an individual species could be calculated using the following [McQuarrie, 1976]: a ln(QN )) 13.45 2 E'=NkT ‘ [ at . D.46 (Di Aei C(-BA81)] 236 The individual energy levels could be used to calculate the overall energy of the plasma as: 13:22 Xi E1 D.47 i Energy may also be expressed in terms of enthalpy. The enthalpy of an individual species could be calculated using the following [McQuarrie, 1976]: 1 D48 8T Hi=Ei+NkT D49 Likewise, the overall enthalpy of the plasma could be calculated as: H=E+ZRT D.50 D.3.1.8 Entropy. In addition to knowing the energy of the plasma, Which allows one to obey the First Law of Thermodynamics, one must know the entropy of the plasma. The Second Law of Thermodynamics could not be violated in the modelling of this plasma fluid flow. Entropy could be calculated for each species using the following [McQuanie, 197613 237 si =NkT[alg($N)]+kln(QN) D.51 Si —" + I N This could be solved and expressed in terms of partition functions [Dow Chemical, 1969]. 1.43879R -A- (’BASi) 1153 Si=-23-Rln(Mi)+%Rln(T)+ T 2‘”! 813 +ln(qe1ect)+C i Gelect Using the formula from the JANAF Tables, the constant (C) was -l.164956. This caused the entropy to be zero at a temperature of 0 Kelvin. This condition satisfied the Third Law of Thermodynamics. Similar to energy, the total entropy could be expressed as: 8:22 Xi Si D.54 l D.3.l.9 Chemical Potential. Assuming chemical equilibrium, the stoichiometric (vi) sum of the chemical potentials (114) must be equal to zero. i 238 Using statistical mechanics, the chemical potential could be expressed as [McQuarrie, 1976]: “i =-kTaln(QN) 13.56 M Using partition functions, this reduced to: . 13.57 pi Lung?) Mathematically, this reduced to: 3 D58 2 w] n -km(q,,,,,)+mn(p) lli='lefl D.3.1.10 Heat Capacity. To model the changes in energy levels with temperature changes, one must know the heat capacity of the fluid. This parameter was defined as the change in enthalpy with respect to temperature at constant pressure. ’ 8T p 239 The heat capacity of the individual species was calculated as [McQuarrie, 1976]: C” i clelect _ 5 d [ 2 mi As: 4'3 A81 )] D.6O Unlike that of energy and entropy, the overall heat capacity must account for the chemical reaction and compressibility changes. The overall formula for this came down to (Lick, 1962]: ' Z H - D.6l (Cpl=22 Xi Cp,1+zz[di] Hi l >2[dx‘] M1 ldT MoidT D.3.2 Boltzmann Equation. The Boltzmann distribution function was a function of position, velocity, and time [f(r,v,t)]. A mathematical derivative identity of this function was [Cherrington, 197911 df are: are} are? D.62 __=____+_=__+_=__ dt atat are: avat The following were identities for substitution into the above mathematical identity: 240 a t identity D.63 _ = 1 a t a} — velocity D.64 — = v t a __ V f radial gradient of function D65 3 r r E _ V f velocity gradient of function D66 3 3 V a 3 F acceleration D.67. at - m d f 8 f collision term D.68 d t - 6 t C Substitution of the above yielded the Boltzmann equation: - D.7O 6f - F 5f — — V f= — at+v0Vrf-i-m- v [5‘1 This equation was a conservation equation of the function over a differential element. 241 D.3.3 Conservation Equations. One obtained the conservation equations by integrating the Boltzmann equation over the velocity space. The Boltzmann equation was first multiplied by a phase function, ¢(v), before the integration. The following was that integral [Cherrington, 1976]: .. — _ _ D.7l m qt?) flaw-v.1: +5-vvr d3v=m¢6 5—f .13., — o" t m — 5: V v C D.3.3.l Continuity Equation. MASS COULD BE NEITHER CREATED NOR DESTROYED. For the continuity equation, one would set the phase M function to one. After integration, one would obtain the following equation: _ D.72 a—n—+Vonvo=<§—f—n> ~ 3 t 6 t C D.3.3.2 Momentum Equation. THE TIME RATE CHANGE OF MOMENTUM OF A BODY EQUALS THE NET FORCE EXERTED ON IT. For the momentum equation, one would set the phase function to mass times the velocrty. After integration, one would obtain the following equation: - ._ D.73 M+Vc(nm)»Fn=<§-:nm V>C 3t St 242 D.3.3.3 Energy Equation. ENERGY COULD BE NEITHER \— CREATED NOR DESTROYED' IT COULD ONLY CHANGE IN FORM. For the k energy equation, one would set the phase function to mass times the square of the velocity divided by two. After integration, one would obtain the following equation: _ _ — - — r — D.74 576m.+gp].%v.hn(a)+mnuuu.spv,]-nF.v=%<§_tnmw>c D.3.4 Collisional Processes. The Collisional processes were affected by the interactions between particles. The following five binary interactions were provided for a thorough discussion of their effects upon the transport phenomena within this model. Tertiary and higher level collisions were not provided because their frequencies were much less than binary and will be neglected from the model [I.ick, 1965]. D.3.4.1 Neutral-Neutral. The collision reaction involved the collision between He and another He molecule. Many potentials existed for these reactions. Only the classical potentials were analyzed. 0 The Hard Sphere potential was the least realistic. It assumed no attraction beyond its radius and an infinite repulsion within the radius. This potential was defined as [McQuarrie, 1976]: 243 = °° V r<6 13.75 (D {0 V r>o o The Square Well potential was a little more realistic than the Hard Sphere model in that it allowed for a transition period of attraction beyond the radius. This potential was defined as [McQuarrie, 1976]: cc V rlo 0 The Lennard-Jones 6.12 potential was more realistic than the previous ones in that it assumed a decreasing level of attraction as the radius increases. It also predicted a maximum attraction at a specific radius. This potential was defined as [McQuarrie, 1976]: mitt-1%)] Experimental parameters for helium were provided from previously obtained data: a = 1.52 x 10'22 J rm = 2.963 x10"10 m o o: 2.57 A 244 A plot of these three potentials was provided in Figure D.2. D.3.4.2 Neutral-Ion. The collision reactions involved the collision of He with the He+ and H6” ion. For slightly ionized plasmas, the He-He'H' collision could be neglected. Therefore, only the reaction between He and He+ were analyzed. A weak interaction exists. This was the polarizability of He by the He+ and was defined as [McQuarrie, 1976]: 2 D78 Whereas, 0t was defined as the polarizability of He. 0 a = 0.2051 A This potential was not important until high temperatures were obtained. The other interaction was charge exchange- He + He+ —————> He+ + He 245 The difference between this and the neutral-neutral collision was the collision cross- section. The potentials used would be the same. D.3.4.3 Neutral-Electron. The collision reactions involved the collision of He with an electron. This collision was very similar to that of the Neutral-Ion collision. Likewise, a weak interaction existed for the polarizability of helium and was not important until high temperatures were obtained. This collision was treated like that of the Neutral-Neutral one with a different collisional cross-section. D.3.4.4 Charged Particles. The reactions involved the collision of He+ with another He+ ion or electron, and of electron with another electron. Regardless of the species, the potential for the collision involved Coulomb forces. The Coulomb potential was defined as [McQuarrie, 1976]: D79 (”$11 4neor Whereas, AD was the Debye length and Q was the charge. Because the electron velocity was much faster than that of the ions, the ions were considered stationary. Therefore, the Debye length could be defined as [Nicholson, 1983]: 246 D80 113 =[ kTe 80 f6 D.3.4.5 Excited Species. The reaction involved the collisions of Heat and He+*. The potentials for these reactions should follow that of the neutral- neutral collision. However, the collisional cross-section would be larger [I.ick, 1965]. Nevertheless, the effect of excited species was small; because at temperatures below 20,000 Kelvin, there were a negligible amount of them. D.3.4.6 Collision Integral and Rates. As mentioned earlier, the Boltzmann equation was important for solving the transport phenomena within this model. An important component of that equation was the collision integral. This integral was needed to calculate the transport coefficients. This integral was defined using Sonine polynomial expansion coefficients. The values "I" and "s" were dependent upon the considered transport coefficients and the degree of approximation. This integral was defined as [McQuarrie, 197613 D.81 m Q(l,s)=(41tkT)% T T e[‘72) V(25+3)(1—coslx)b db dy 0 0 247 whereas "b", "x", and "7" were respectively impact parameter, deflection angle, and reduced initial relative velocity. The deflection angle and reduced initial relative velocity were defined as [McQuarrie, 1976]: m ]}5 D82 oo D.83 xzn-Zb] dr y rm 2 4(1) b2 2 ’ 1‘ 2'7 mv r Figure D.3 contained the scattering plot geometry of this deflection. D.3.5 Transport Coefiicient The calculation of these transport coefficients required the collision integral as previously discussed. For this model, five of the coefficients were important. They were electrical conductivity (0), thermal conductivity (2»), mobility (K), viscosity (1]), and diffusion coefficientw). 248 Impact Parameter _ ____fiz _______ Deflection Angle Figure D.3 Collision Path D.3.5.1 Electrical Conductivity. This referred to the ability of the material to conduct electricity. It was a scalar multiple of the electric field (E) which related to the current density [Panofsky, 1962]. J=oE+VxB D34 249 Using the Langevin equation, a model using elementary gas discharge interactions (Nicholson, 1983]: dit Y)=-e(E+YxB)-mvvm D85 the electrical conductivity could be calculated from the following [Chenington, 1978]: 2 V3afo D.86 °° av ‘dV 3m 0vm+jco Mn 6:- D.3.5 .2 Thermal Conductivity. This referred to the ability of the material to conduct thermal heat. It was a scalar multiple of the gradient of temperature when related to the heat flux (q) [Bennet, 1974]. q=_WT D.87 D.3.5.3 Mobility. This referred to the steady-state velocity attained by an object under the action of an external unit force. According to Anderson, one could relate the drift velocity of charged particles (such as ions) to the electric field by [Anderson, 1990]: 250 dexE D88 Thus, one could relate the mobility to the electrical conductivity by: D.89 D.3.5.4 Viscosity. This referred to the fluid’s resistance to flow. The force that caused this resistance was directly related to the differential change in velocity with respect to position change by this scalar value [Bennet, 1974]: F=-nv‘\7 ' 13.90 D.3.5.5 Diffusion Coefficient. Once one knows how electromagnetic and thermal energy was transmitted and the mobility and fluid resistance of particles in the fluid flow, one would still need to know in which direction particles want to go. Diffusion coefficients provide one with that information. These coefficients provided a direct relationship between molar fluxes (F) and concentration gradients in space. The molar flux of a particle could be expressed as [Cussler, 1984]: fi = Ci (Yr -) D.9l 251 with "ci" being the concentration of species "i" and "" being the average molar velocity. According to Fick’s First Law in a binary system, the molar flux could be expressed as [Cussler, 1984]: E1 =-cD1,2 Vxl 13.92 For example of electrons lost by diffusion to cavity walls and using the combination of the continuity equation with the momentum equation, one could obtain the following mathematical expression for the number of electrons (rte) [Chenington, 1979]: 1 3616?) 13.93 Vm o’t ne V-t’ne Kc E’De Vile ' with the "vm" identified as the collision frequency. Assuming constant temperature and steady-state conditions, one could approximate the electron diffusion coefficient to be: kTe 13.94 De = me Ve This diffusion coefficient was based upon free electron diffusion, which could only occur at low pressures when coulomb effects could be neglected. At higher pressures, the ions could affect the flow of the electrons. Likewise, the diffusion coefficient for ions could be approximated to be: 252 k1", 13.95 mi Vi Because the mass of the electron would be much smaller than that of the ions and neutrals, the electrons would diffuse away at a much faster rate. However, the electrical attraction of the ions would hinder their diffusion, but increase that of the ions. As a result, a space charge field (Es) would be established. This ambipolar effect was illustrated in Figure D4. The new fluxes for the species would be written as: f6 =-De Vile ‘ne Ke ES D.96 fi ='Di Vni +ni Ki Es D97 Assuming quasi-neutrality within the plasma, 253 Ambipolar Diffusion Plot (Electron 6:: Ion Diffusion) Particle Dens‘ Radial Distance (from axial center) Figure D.4 Ambipolar Diffusion the continuity equation could be expressed as: a n _ D.98 254 or rewritten as: a _ _ =Di Vzn-nrci VOES 49' ES OVn D.99 By rewriting the equatiOns and adding both equations together, one could obtain the following: a n 2 D. 100 —— = D V n a t a with the ambipolar diffusion coefficient being defined as: _ De Ki +Di KC DJOl Da The same analysis could be done to characterize multiple species diffusion. However, in the pure helium plasma system for temperatures less than 20,000 Kelvin, one would not have multiple species present. As discussed in Chapter 8, one could neglect the double Charged helium ion concentration at these low temperatures. 255 D.4 Chemical Kinetics. The steady-state calculations mentioned previously assumed equilibrium conditions in the chemical reaction. He e—l—(J—e He+ + e' However, a more realistic set of calculations may be done by relaxing this assumption and by considering the kinetics (or reaction rates) involved in the ionization and recombination reactions. -He —51—> He+ + e’ Ionization He+ + e' i—a He Recombination D.4.l Reaction Rate. A first approach calculation of the kinetics involved within the reactions would be to take a maxwellian distribution of the particles. 2 D.102 The collisional rate from equation D.81 could be calculated as [Chenington, 1979]: 256 21: kT “3 “We _mv2 D.103 (G(VM=4n[m J/zm [2kT]d dv Using a maxwellian energy distribution function defined as [Chenington, 1979], f(,)_ 2i: ]}/ [ )721 117;] 6.104 kT The ionization collisional rate (in m3/sec) could be approximated to [Cherrington, 1979]: -30 ] 13.105 («wimp {111—ij a, [$ij [lfilik—Tc‘ C The ionization energy (so) and the ionization cross sectional area (00) were extracted from Cherrington. -l9 50:8.011 x 10 J 00: 7.0x1014c m3 257 The ionization reaction rate (Rion) in numbers per second could be calculated from the ionization collisional rate rate (ion) by multiplying it by the density (p) of the species. Rion = PHe (OM/hon D.106 Likewise, the recombination reaction (Rrec) could be determined: Rrec = pHe+ (0(v)v>rec D'107 At equilibrium, the following relationship holds: RlOIl =RI'CC D.108 Therefore, one could calculate the recombination collisional rate using the following rearranged relationship at equilibrium: '1 D.109 __ p‘13qu 1 (G(V)’>ion (C(V)v>rec - pequil He+ 258 As mentioned previously, double ionization for low temperatures (less than 20,000 Kelvin) could be neglected. Thus, the ion (He+) density would be equal to the electron (e‘) density. One could calculate the recombination rate using the the following: uil D. l 10 p Ilqe (06,)" > ion (6(v)v> = . rec equrl P _ e D.4.2 Reaction Time to Equilibrium. An important element in these rate calculations would be to determine the time required to reach equilibrium conditions. This would be important to determine how close the model simulations would be to the actual conditions. A simple first approach calculation could be done by taking a small elemental volume within the plasma discharge and treat it as a batch reactor. The overall reaction rate would be equal to the forward reaction minus the reverse reaction. - D.111 de Roverall = Rion ' Rrec = —dt With the substitution of equations D.106, D.107 and D.110 into DJ 1 1, one could obtain the following overall rate: . 2 Roverall = PHe (G(V)V>jon ' pe’ ion 'Pe- (0(V)V>rec This results in the following solution. -vm[p 066M... ape- <00»)... we. con...) Pe- ”“6 2 (o(v)v)l-on + (0(V)V>rec 0 260 Starting with a zero electron density, the following solution results. -(2v>..,+v>...lt W V l-e p3- = 2 (0(V)V>ion + <0(V)V)rcc 261 APPENDIX E FLUID TRANSPORT PHENONENA The conservation laws for fluid flow characterized the transport phenomena for the fluid. For example, the continuity equation characterized the conservation of particles. Similar things could be said of the motion and energy equations. E.1 Flow Functions. For analysis of fluid flow, the velocity of the fluid particles was an important parameter. The conservation laws were generally derived using the velocity parameter. For some situations, it may be easier to transform the velocity function. For example, a two component velocity could be transformed into a one component function (w). This function was known as a "Stream Function." The following were transformation equations to stream functions for various coordinate systems [Bird, 1960]: 262 COORDINATE mm Rectangular Cylindrical (r, 0) Cylindrical Spherical VELOCITY COMPONENTS _ -1 8w i-Zsin(0)ae 1311! V 9 = rehab)? 263 13.] B.2 E3 E4 135 E6 E7 E.8 In a physical sense, the stream function characterized the average path of a particle in a flowing fluid. Another important flow function was the "Vorticity," (a). The vorticity was the curl of the velocity. w=VxV E9 In a physical sense, the vorticity was twice the angular velocity of a particle in a flowing fluid. E.2 Conservation Laws. The following laws were generally used to characterize the transport phenomena of fluid flow. The laws were the same as those described in the Plasma Transport Phenomena chapter previously. E.2.1 Continuity. This equation was derived by writing a mass balance over a differential element. The accumulation rate of mass was equal to the difference of the rate of mass into and out of the element. This equation could be mathematically written as [Bird, 1960]: §a_f=-V.¢);) EJO 264 E.2.2 Motion. This equation was derived by writing a momentum balance over a differential element. The accumulation rate of momentum was equal to the difference of the rate of momentum into and out of the element plus the sum of forces acting upon the element. This equation could be mathematically written as [Bird, 1960]: 5%tz)=-v.(osc)+i= 3“ The force, F, could be seen as the combination of pressure, viscous transfer, and gravitational forces. In mathematical terms, this was: F=-VP-Vo:+pg 512 The "t" was known as the stress tensor. If the divergence of that tensor was zero, one obtains the famous Euler equation (first derived in 1755). For Newtonian fluid, the stress tensor terms could be expressed as: E.13 Tisi :-21]-85vi—i+(-§-ll-C)(VOV) . E.l4 ‘t- --—_ ...—...avi +9.3]; “1' " ii] iii 265 The "C" term was the bulk viscosity term, which would be identically zero for low density monatomic gases. The equation of motion could be rewritten as: aithF-V-QW)-VP-V.:+p§ E-15 For constant density (p) and viscosity (n), the above equation could be reduced to the famous "Navier-Stokes" equation. E.2.3 Energy. This equation was derived by writing an energy balance over a differential element. The accumulation rate of energy was equal to the difference of the rate of fluid energy into and out of the element plus the sum of energy transported to the element. This equation could be written mathematically as [Bird, 1960]: 2 2 _ _ E.l6 % p[c+L2—]]=-Ve[p[e+x2—]v]+E The added energy, E, could be seen as the combination of pressure, viscous transfer, gravitational and heat transfer] generation energy. Mathematically, this was: E.l7 E=-v.(o?)-v-C.C)+pG-§)+pd 266 For calorically perfect gases, the internal energy (e) could be written as: e=CvT ' EJ8 Thus, the overall energy equation could be written as: ait[p[c,r+123]]=-v.[p[cvr..gflmqstvlac).peg)”; E.19 E.3 Compressible Fluid Flow Variables. Basic assumptions were normally made when analyzing physical systems. In the case of fluid flow, it was often desired to assume incompressible fluids. This was normally valid for liquid fluid. However, for gaseous and plasma fluids, incompressibility could not be assumed. To reduce the complexity and rigidity of the conservation equations, transformations to new variables could be done. For compressible fluid flow, sound speed (a) v mach number (M). and heat capacity ratios (Y) were important. The heat capacity ratio was defined as [Anderson, 1990]: 13.20 Cp 267 E.4 Method of Characteristics. A powerful numerical method in solving compressible fluid flow equations was the "Method of Characteristics." This method assumed two compatible equations intersect each point in space. To get these equations, one must rewrite the square of the sound speed [Anderson, 1990]. 2 BP E.21 a = — lap. For one-dimensional flow, the continuity and momentum equations could be re-written respectively as [Anderson, 1990]: 1 3P 3P 3v E22 ——-——+v—— —-=0 a2(3t 8x 8x 3v av 123:0 E23 at VEI 63x Adding these equations together yielded: E24 268 On the other hand, subtracting these equations yielded: [2+ v-a 3:]-—1—-[§—l:+ v-a 3_P -0 E25 at 3x pa at 3x ‘- Respectively, a solution of these two equations yielded the following two differential relationships between "x " and t"": dX V 'l- a 13.26 9d- (Eight/die v —- a) t """""""""""" :33“; C+ characteristic line 1 1’ . \‘ +(dX/dt= v + a) 1’ I s I I s I j x 2X Figure B] Characteristics 269 In fact, this relationship was the step-size physical definition of the time differential of position (x). These two separate relationships formed the basis of the two characteristic equations intersecting a point in space ( see Figure B. 1). It was more advantageous to use pressure and density instead of position and time. Thus, these two characteristic equations could be written as: dv+-d£=0 E.27 pa d O E.28 v--——= pa Knowing information on some points in space, one could numerically predict other points using these characteristic lines. Integrating these equations along the characteristic lines, one could obtain useful constants, commonly called "Riemann Invariants." dP E29 J+ =V+I— pa (1? E30 J_=v-I—— pa 270 F . . . . or deterrmmng flow patterns wrthrn supersonic nozzles, method of characteristics could be sed. ' ' u As portrayed in Figure B.2, the nozzle could be subdivided and calculated by interlaping the lines. At the throat with Mach 1 being assumed, one could calculate the fluid profile within the nozzle downstream thereof. 3 / 4 A/ § centerline Figure B.2 Characteristic Nozzle 271 APPENDIX F RADIATION TRANSPORT PHENOMENA A major loss mechanism in energy transport to the propellant was energy loss to the walls of the discharge chamber. Energy could be transported through conductive, convective, and radiative means. This chapter only analyzed the radiative energy losses to the wall. The following parameters were used for this analysis: the cavity material was unpolished brass, the plasma geometry was an oblate ellipsoid, and the wall temperature was 300 degrees Kelvin. Table Fl Radiation Model Parameters TYPE PARAMETER Cavity Geometry: Cylinder Cavity Diameter: 17.8 cm Cavity Lengths: 7.2 cm (TM011 mode) 14.4 cm (TM012 mode) Cavity Material: Unpolished Brass Plasma Geometry: Oblate Ellipsoid Propellant Flow: 572 SCCM Propellant Gases: Helium and Nitrogen Wall Temperature: 300 K 272 BODY W BODY NUMBER NUMBER Q Figure F.1 Radiation Heat Transfer Model F.l Blackbody. A blackbody was a body that absorbed all incident radiation. No reflection occurred. Additionally, a blackbody was a body that emits radiation based upon its temperature. This emitted radiation increases with temperature so that energy transfer between two bodies was from the lower temperature body to the higher temperature one. The energy transport was illustrated in Figure F.l [Siegel, 1981]. 273 El=f(T1) RI 132 = f(T2) F.2 Q=Ei-E2=f(T1)-f(T2) F3 In addition to temperature, radiation energy was also a function of wavelength 0»), direction (941)), and surface area (A). Therefore, this energy could be written as the following function: E = f(T,7L,9,(p,A) 13.4 Assuming we had a perfectly smooth surface, the directional component of this function dropped off. Through quantum calculations, we could obtain Planck’s spectral distribution of emissive power [Siegel, 1981]. 21tclA F.5 _ ['2 _ ASIAT -1 _ l 274 BLACKBODY RADIATION EN ERGY (Temperature .— 300 K) 3.5 3.0- TM ’3' - 012 '5 2.5- E . E8 2.0— g . E\ ‘°°" TMO11 % 1.0.. 35. - 0.5- .J 0.0 ._.++,fi,.,.,.,_,r j) 4- % é 1o 12 14 15 1a 2022 (WAVELENGT micrometer Figure F2 Blackbody Radiation Energy The CI and c2 constants were defined as: c1=hcg F.6 h F.7 c2: CO k A plot of this energy versus wavelength for a temperature of 300 K and area for both TM modes was provided in Figure F.2. The total energy emitted by the body was the area under this curve. Mathematically, it was the integral over all wavelengths of the wavelength dependent energy. 275 E(T)=TE(T,i)di F'8 o 00 F9 E(T)= J. ZKCI A d)“ 0 15 (1%] - 2- _ l Using the following transformation of variables: é__-c_2 F.10 AT The above integral became: 4 oo 3 Fl] 2 T A E(T)= 7:01 5 I dé C; 0 e(§) _ 1 This, then, reduced to: 4 4 F.12 T A E(T)= 201 7‘ 15c 276 This expression could be simplified by redefining a new constant. The following new constant was used: 5 RB We now had energy defined as a function of temperature with a new constant (0), known as the Stephan-Boltzmann constant [Siegel, 1981]. E(T)=O'T4 A F.14 Assuming blackbody radiation only within the experimental system, we could obtain an estimate for the average surface temperature of the plasma. This could be done for both the strong and weak surfaces. The two body system of concern was the cavity wall and the plasma (see Figure F .3). The energy (QC) absorbed by the cavity wall was: QC =Ep 'EC F.15 Using the above energy expression, F.16 4 Qc =°(APT3 'ACTCJ 277 Figure F3 Radiation Model Schematic 278 The area equations used were that for a cylinder and an oblate ellipsoid [Perry, 1984]. 2 W L 2 219 1-19 [3.18 Whereas, "19" was the eccentricity of the oblate ellipsoid and was defined as: 2 2 F.19 p 'Wp Lp 19: Quick calculations of this suggest that the average plasma surface temperature was about 1000 K. Unfortunately, this blackbody approximation to the system was not quite accurate. The plasma temperature was expected to be higher. F .2 Graybody. A graybody was a body that absorbed a fixed fraction of radiation. Reflection of radiation occurs. This required two new terms, call emissivity (e) and the other absorptivity (or). Emissivity was the fractional measurement of how well a body could 279 radiate energy when compared to a blackbody; whereas, the absorptivity measured the fractional amount of how well a body could absorb energy when compared to a blackbody. These values could have the following range [Siegel, 1981]: 03831 F20 This emissivity and absorptivity were functions of wavelength, temperature, and direction for a given body. a = fOL, T, 9,< I III an Next, one should calculate the J acobian of the function, ""F, which was the derivative matrix of the function set with respect to the "x" matrix. (1 F x ml) 6.9 -d; with each matrix element defined as: 285 3%) £113) %Q (3.10 1(8) (MI) (“2%) fldxfi") ‘ (1’51 4’52 . at) at) ..j «.9 _ dxr dxz an J Using the Newton Method, one could iteratively calculate the succeeding sets of solutions. This should quadratically converge [Johnson, 1982]. Because it was inefficient to calculate the inverse of the J acobian, one should rearrange the Newton Method equation to produce a less intense computational algorithm [J ohnson, 1982]. J[;“][;"+1-;“]=-r{§") 6.12 By defining a new variable set: —n+1=;n+1_-n G.13 one could obtain the following: I?! Jam = “III 0.14 286 which was in the same form as the linear system of equations. Thus, each iteration could be solved by using the techniques need for solving linear systems. I used this method for calculating the thermodynamic properties using statistical mechanics (see Appendix A.5). G.2 Data Curve Fitting. This referred to approximating a function to duplicate experimental (or actual) data. For nth order polynomial interpolation, one could approximate the function [Johnson, 1982]: n . f(x)=g(x)= 2W1 xl'1 G.15 i=1 Then, one could solve the following: G.16 >< ll 2| II '11! If the "x" matrix was not square (i.e. i x i), one could make it square by multiplying it by its transpose matrix. This was commonly referred to as a Least Squares technique (see Appendix A.2) [Johnson, 1982]. F G.17 Another popular type of equation would be orthogonal polynomials. These functions must satisfy the following integral relationship [Finlayson, 1980]: 287 b IW(X)gn (x)gm(x)dx =0 V n¢m and g(x)Ea—)b (3.18 a The orthogonal equation that I looked at was the Chebyschev Polynomials of the First Kind [Abromowitz, 1964]. Ti+1(x)=2xTi (x)-Ti-1(x) v 121 (3.19 with the first two defined as: To (2051 T1 (x)2 x G.20 or, these polynomials could be expressed using the following trigonometric function: Ti(x)=cos(icos'l(x)] V -leSl (121 G3 Classification of Partial Differential Equations. Partial Differential Equations (P.D.E.’s) could be classified in various forms. Each form had a different approach for numerically solving the equation. For this dissertation, a second order linear P.D.E. was used to identify three different types of equations. A generic P.D.E. for two variables had the form of [Hal], 1990]: 2 2 2 3f+b af +c3f+d§i+cfi+gf+h=0 G22 8x2 8x3y ayz ax 3y a 288 The types were categorized according to the following conditions on the above equation: b2 -4ac >0 hyperbolic b2 - 4 a c < 0 elliptic b2 -4ac=0 parabolic G.4 Galerkin Method. One method for solving P.D.E.’s was a procedure in which one used a set of trial functions [Finlayson, 1980]. was; a. MK) 6.23 For the Galerkin method, those functions had to be orthogonal to one another. Mathematically, this satisfied the following: I401 ¢jdX=O Vigil G24 For functions with more than one variable (as in the case of many P.D.E.’s), the set of trial functions were approximated to be separable. f(x,y)=_2 Z ai, j (0i (’0ij G25 1 1 Examples of orthogonal trial functions included the following. These functions had to satisfy the boundary conditions of the problem or domain. 289 P, = cos (i n X) Harmonic G26 P = cos icos'l (x) - 1 Chebyschev Polynomral G27 i ’-1 P1 = a X + b X1 + ”’+ Z Legendre Polynomial 6.28 (using Gram-Schmidt algorithm) To solve the P.D.E. over the desired domain, one could transform the P.D.E. into a linear set of algebraic equations with unknown trial function coefficients (ai). Illustrating this method, the following simple example was provided. Function: F = f(x) Domain: x = 0 —> 1 B.C.’s: f(0) = f(1) = 0 Trial Function: (Pi = 5i“ (i 7‘ X) O.D.E.: £35 + x SE = G(x) dx 2 ‘1" Number of Nodes: N TRANSFORMATION: (to set of algebraic equations) 2 N 2. N d F 2 3i d $1 =- 2 ai(i1t)2 sin (inx) (3.29 clx2 i=1 dx2 1:1 2 6.30 3i 5%?— : ai (in)cos (in x) %|% ”NZ II'M 1 1 1 1 One would multiply each side of the O.D.E. with one of the trial functions; and, one would do this for each function to generate a set of N algebraic equations. 2 I —2+Xd—x- ¢idx= 2 IG(X)¢idX 0'31 j=10 dx j=10 In matrix form, this reduced to: Isri S1,2 SLN'Iai' Ibr‘ G32 82,1 82,2 52,N a2 b2 _SN,1 SN2 SN.N]_aN1 The unknown values were the coefficients (ai) for the set of trial functions. The above set of integrals reduced to a set of algebraic equations, which could easily be solved by such methods as Gauss elimination. G.5 Finite Difference. This method for solving P.D.E.’s was much easier to apply. One would approximate each derivative using the Taylor Series Expansion of a function about a nodal point. The following were examples of approximations for derivatives (using h = A x and g : Ay) [Hall, 1990]: (”(8):- f(a + h)- f(a) forward G33 dx h 291 df(a)_=_ f(a)-f(a -11) h dx backward G.34 df f + h - f - h (if); (a 2) h (a ) centered (3.3 5 For second order derivatives and two variables, the following centered differences could be used respectively: d2f(a)= f(a+h)+f(a-h)-2f(a) G.36 dx2 112 22am): f(a+h,b+ g)+f(a-h,b- g)-f(a+h,b- g)-f(a-h,b+ g) (3.37 Bxa’y — 4hg Solving P.D.E.’s using this method was similar to that of Galerkin methods in that one must transform the P.D.E.’s into a linear set of algebraic equations with the functions at the nodal points being unknowns. Using the same example as before (see the Galerkin example), this method could be transformed as: 2 N F: dF' N G38 2 d J+x-——J—= G(Xj) ._ 2 J dx -1 j—l dx _] whereas, Fj §f(xj) v x0 =0,xN =1,andXi-1