INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproductioa In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9" black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. University Microfilms International A Bell & Howell Information Company 300 North Zeeb Road, Ann Arbor, Ml 48106-1346 USA 313/761-4700 800/521-0600 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R eproduced w ith perm ission o f the copyright owner. Further reproduction prohibited w itho ut perm ission. Order N um ber 9127441 Analysis o f solar and microwave therm al propulsion system s Venkateswaran, Sankaran, Ph.D. The Pennsylvania State University, 1991 UMI 300 N. Zeeb Rd. Ann Arbor, MI 48106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R eproduced w ith perm ission o f the copyright owner. Further reproduction prohibited w itho ut perm ission. The Pennsylvania State University The Graduate School Department of Mechanical Engineering ANALYSIS OP SOLAR AND MICROWAVE THERMAL PROPULSION SYSTEMS A Thesis in Mechanical Engineering by Sankaran Venkateswaran Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy May 1991 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. We approve the thesis of Sankaran Venkateswaran. Date of Signature Z-H-V Charles L. Merkle Distinguished Alumni Professor of Mechanical Engineering Thesis Adviser Chair of Committee c f- / S ’ f t Harold R. Jacobs Professor of Mechanical Engineering Head of the Department of Mechanical Engineering Stefan T. Thynell Assistant Professor of Mechanical Engineering 3-H'S/ J^Hn M. Cimbala Associate Professor of Mechanical Engineering \<L j Robert J. Heinsohn Professor of Mechanical Engineering 12 A ^I William G. Pritchard Professor of Mathematics Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ABSTRACT Theoretical models are developed to assess the conceptual feasibility of em ploying solar and microwave energy to heat a propellant gas and generate thrust. Both these advanced electrothermal propulsion concepts promise high specific im pulses and are, therefore, being considered as alternatives to chemical propulsion for space-based applications. The present research addresses two major feasibility issues: (1) the efficient coupling of electromagnetic energy to the propellant, and (2) the performance capabilties of these systems. Physical modeling of the direct absorption of solar radiation in gases involves the solution of the thermal radiation problem within the thruster. The absorption of the incident beam is modeled by a ray-tracing technique while the reradiation of the heated gas is treated by a thick-thin model. The thick-thin model is calibrated by the P -l model and an exact integral method. Spectral dependence of the radiative properties is handled by a band model. Results for different thruster sizes indicate th at high energy conversion efficiencies are possible if a combination of direct and indirect heating is used; part of the incident radiation is absorbed directly by the gas while the remainder is absorbed on a regeneratively cooled wall. Performance predictions for the thruster indicate specific impulses in the range of 500 to 650 s for moderate gas pressures (1 to 10 atm); at higher pressures, better performance is possible but structural considerations are likely to play a limiting role. The analysis of the microwave problem involves the coupled solution of the Navier-Stokes equations for the gasdynamics and the Maxwell equations for the mi crowaves. Time-marching techniques are employed for the solution of both equation sets. Local thermal equilibrium is assumed for the plasma. Computational results are in good agreement with experimental data for resonant cavity plasmas — peak Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. plasma temperatures are between 9000 and 10,000 K and high coupling efficiencies (up to 100%) are predicted. Parametric studies are performed to extend the results to high power and high pressure regimes. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TABLE OF CONTENTS Page LIST OF FIGURES.......................................................................................... viii LIST OF TABLES................................... xiii LIST OF SYMBOLS........................................................................................ xiv ACKNOWLEDGMENTS ............................................................................. xvii Chapter 1. INTRODUCTION................................................................................ 1.1 1.2 1 Advanced Propulsion Systems ............................. ................ Electrothermal Thrusters ...................................................... 1.2.1 Resistojets ................................................................... 1.2.2 Thermal Arcjets ......................................................... 1.2.3 Laser Thermal P ropulsion......................................... 1.2.4 Solar Thermal Propulsion ......................................... 1.2.5 Microwave P ropulsion................................................ 1.3 Review of Modeling in Laser P ropulsion............................. 1.4 Direct Solar Thermal P ropulsion.......................................... . 1.4.1 Review of Modeling of Direct Solar Propulsion................................................... 1.4.2 Synopsis of Present Solar Research ......................... 1.5 Microwave Thermal Propulsion ............................................ 1.5.1 Coupled Navier-Stokes and Maxwell Solutions. . . . 1.5.2 Synopsis of Present Microwave R esearch ................ 1.6 Overview of the Thesis ........................................................... 1 4 '4 5 5 7 8 10 14 2. SOLUTION OF THE MAXWELL EQUATIONS .......................... 27 2.1 2.2 Electromagnetic Field T h e o ry ................................................ 2.1.1 The Maxwell Equations ............................................. 2.1.2 Electromagnetic Force and E n e rg y ........................... 2.1.3 The Electromagnetic Wave Equation ...................... 2.1.4 Microwave Waveguides ............................................... 2.1.5 Microwave R esonators................................................. Time-Dependent Solutions of the Maxwell Equations — 2.2.1 Introduction to Maxwell Solution Methods ........... Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 17 20 20 22 25 27 27 30 32 33 37 38 38 vi 2.2.2 Vector Representation of the Maxwell Equations............................................................... 2.2.3 FDTD Solution Procedure ...................................... Characteristic Boundary C onditions..................................... Electrical Conductivity of the P la s m a ................................... 2.4.1 Electrical Conductivity Expressions........................ 2.4.2 Modified Maxwell Equations ................................... 2.4.3 Equilibrium Ionization .............................................. 2.4.4 Collision Frequency ................................................... 2.4.5 Electrical Conductivity Calculations ...................... Demonstrative Solutions of the MaxwellE q u a tio n s 42 43 45 51 51 53 54 54 55 57 SOLUTION OF THERMAL RADIATION TRANSFER ............ 71 2.3 2.4 2.5 3. 3.1 3.2 Direct Solar Absorption in G a se s........................................... Radiative Properties Modeling .............................................. 3.2.1 Dimer and Atomic Line A bsorption......................... 3.2.2 Photoionization .......................................................... 3.2.3 Inverse Bremsstrahlung ............................................ 3.2.4 Complete Absorption Coefficient ............................ 3.2.5 Band Model for Non-Gray Gases ............................ Equation of Radiative T ransfer............................................... 3.3.1 The Thermal Energy Equation ............................... 3.3.2 The Radiative Transfer Equation ............................ 3.3.3 Direct Solar Absorption ............................................ Solution Techniques for Radiative T ran sfer......................... 3.4.1 Exact Integral Method .............................................. 3.4.2 The P -l Differential Equation Model ..................... 3.4.3 The Thick-Thin Model .............................................. Demonstrative Solutions .......................................................... 71 75 75 77 78 79 87 87 87 88 90 94 94 95 97 100 THE FLUID DYNAMIC MODEL ................................................... 109 3.3 3.4 3.5 4. 4.1 4.2 4.3 Solution of the Navier-Stokes E q u a tio n s............................... 4.1.1 The Navier-Stokes Equations .................................... 4.1.2 Time-Iterative Solution P ro ced u re ........................... Low Mach Number Formulation ............................................ Demonstrative Flow Field S olutions....................................... 4.3.1 Sphere-Cylinder Configuration................................... 4.3.2 Bluff Body Configuration............................................ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 109 109 112 115 117 117 125 5. FEASIBILITY STUDY OF DIRECT SOLAR PROPULSION .. 131 5.1 131 131 135 139 146 146 149 154 5.2 5.3 6. Direct 5.1.1 5.1.2 Direct Direct 5.3.1 5.3.2 5.3.3 5.3.4 Solar A bsorption.......................................................... Sizing of the Solar Thruster .................................... Issues in Direct Solar Propulsion ............................ Solar Propulsion Using Particle A bsorbers............. Solar Propulsion with Alkali V ap o rs........................ Details of the Analysis .............................................. Flow Field C haracteristics............................ Efficiency of Direct A bsorption............................... Comparison to One-Dimensional Predictions............................................................. 5.3.5 Solar Thruster Performance E stim ates................... INVESTIGATION OF MICROWAVE ABSORPTION IN GASES................................................................................ 6.1 6.2 157 160 165 Resonant Cavity Microwave P la sm a s................................... Basic Computational Results ................................................ 6.2.1 Microwave Plasmas in a Straight D u c t ................... 6.2.2 Plasmas in Sphere-Cylinder Configurations...................................... 6.2.3 Bluff-Body Stabilized Microwave Plasm as.................................................................. Model Validation with Experiments ..................................... 6.3.1 Coupling Efficiency..................................................... 6.3.2 Peak Temperature ..................................................... Effect of Parameter Variations on Microwave Plasmas ........................................................ 6.4.1 Effect of Variations in Input P o w e r......................... 6.4.2 Effects of Gas Velocity .............................................. 6.4.3 Effect of Gas P ressu re................................................ 6.4.4 Effect of Gas Pressure on Stability Boundary............................................................... 165 168 168 CONCLUSIONS AND FUTURE WORK ...................................... 214 7.1 7.2 Direct Solar P ro p u lsion........................................................... Microwave Absorption in Gases ............................................. 214 219 REFEREN CES................................................................................................. 224 APPENDIX 236 6.3 6.4 7. ABSORPTION BANDS OF ALKALIS............................. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 174 177 181 181 185 188 188 196 200 203 viii LIST OF FIGURES Figures Page 1.1 Schematic of a direct absorption solar thermal propulsion system................................................................................ 18 1.2 Schematic of experimental set-up for resonant cavity plasmas (from Balaam and Micci [28])............................................................ 24 2.1 Field distribution for TM and TE modes in a circular waveguide (from Gandhi [57])............................................................ 36 2.2 Field distribution for TM waves in a cylindrical resonant cavity (from Gandhi [57]).................................................................... 39 2.3 Staggered grid showing locations of electric and magnetic field components within a grid cell.................................................... 46 2.4 Geometry of microwave cavity showing the inlet waveguide and the co-ordinate system................................................................ 48 2.5 Direct current conductivity of a helium plasma at different temperatures and pressures................................................................ 56 2.6 Real part of ac conductivity for helium plasmas at different temperatures and pressures. Field frequency = 2.45 GHz 58 2.7 Assumed plasma profile at the center of the cylindrical cavity. Plotted lines represent log10<r,jc......................................................... 60 2.8 Analytical solution for the TM 012 mode in a no-loss cylindrical cavity................................................................................. 61 2.9 Temporal development of the electric field at a point in the no-loss cavity............................................................................ 62 2.10 Spatial distribution of computed electric fields inside the no-loss cavity. TM012 mode. / = 2.45GHz.............................. 63 2.11 Temporal development of the electric field at a point in the lossy cavity. Peak value of <7,fc = 1 m ho/m .................................... 65 2.12 Spatial distribution of the electric field in the lossy cavity..................................................................................... 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2.13 Temporal development of the electric field at a point in the lossy cavity. Peak value of <Tdc = 100 m hos/m ................................ 67 2.14 Spatial distribution of the electric field in the lossy cavity................................................................................. 68 2.15 Variation of coupling efficiency with cavity length. No-loss resonant length of the cavity, L 0 = 0.148 m ..................... 70 3.1 Black body intensity at the surface temperature of the sun. TSUIl = 5800 K........................................................................................ 72 Inverse bremsstrahlung absorption coefficient for Na,K and Cs. T = 4000 K, P 0 = 3 atm, yca = 0.05................................................. 80 Absorption coefficient of Na vs. wavelength. P 0 = 3 atm , 3*Na = 0.05........................................................................ 81 Absorption coefficient of K vs. wavelength. P 0 = 3 atm , ya = 0.05......................................................................... 82 3.5 Absorption coefficient of Cs vs. wavelength. P 0 = 3 atm, yca = 0.05........................................................................ 83 3.6 Mean absorption coefficient of Cs vs. temperature for different pressures and alkali mole fractions.................................... 85 3.7 Absorption coefficient of a 5% mixture of Na, K and Cs (1.67 % of each alkali) in hydrogen. P 0 = 3 atm ............................. 86 3.8 Solar ray model superimposed on a fluid dynamic grid. 240 rays are shown with a 51 x 51 grid............................................. 91 3.9 Intersection of a single solar ray with the fluid dynamic grid cells................................................................. 93 3.10 Geometry used in the computations with the boundary conditions............................................................................. 101 3.11 Comparison of exact, P -l and thick-thin solutions. k = 0.4 m ~ \ P = 1 MW, C R = 104 : 1 ............................................ 103 3.12 Comparison of exact, P -l and thick-thin solutions. k = 2 to " 1, P = 1 MW, C R = 104 :1 ............................................... 104 3.13 Comparison of exact, P -l and thick-thin solutions. k = 6 t o - 1, P = 1 MW, C R = 104 :1 ............................................... 105 3.2 3.3 3.4 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. X 3.14 Distribution of radiative flux along the peripheral wall. Non-dimensional heat flux, q* = 107 4.1 Fluid dynamic grid for sphere-cylindrical tube with the microwave cavity. Right side of the cavity shows the no-loss field pattern of the TM 012 mode................................................................................... 118 4.2 Convergence of the low Mach number algorithm for cold flow in the sphere-cylider geometry. 81 x 41 grid. Re=80.......................... 120 4.3 Streamlines in plasma tube for cold flow at various Reynolds numbers. Flow is from left to right................... 121 4.4 Effect of specified heat addition on flow through sphere-cylinder configuration. R e=8, P=150 W ......................................................... 122 4.5 Effect of specified heat addition on flow through sphere-cylinder configuration. Re=80, P=300 W ....................................................... 123 4.6 Effect of specified heat addition on flow through sphere-cylinder configuration. Re=80, P=300 W ....................................................... 124 4.7 Fluid dynamic grid for axisymmetric bluff body geometry. Flow is from top to bottom .......................................................................... 126 4.8 Convergence of low Mach number computations for the bluff-body geometry at different Reynolds numbers.......................................... 128 4.9 Bluff body computations for different Reynolds numbers. Re = 5, 60, 350,1000........................................................................... 129 5.1 Power requirements of a solar thruster for different thrust levels and gas temperatures.................................................... 133 5.2 Collector and focal spot sizes as a function of thrust size. C R = 1 x 104 :1 .......................................................................... 134 5.3 Temperature contours at two different inlet gas velocities (a) 0.1 m /s, (b) 0.4 m /s. t = 5^ .................................................... 141 5.4 Radiative flux on the wall for different gas velocities and wall emissivities. q* = 142 5.5 Collection efficiency and peak gas temperature as a function Boltzmann number. <max = 144 5.6 Thrust and specific impulse of the direct thruster using particles for absorption........................................................................ 145 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5.7 Solar thruster geometry with ray model superimposed on fluid dynamic grid................................................................................ 147 5.8 Temperature and Mach number contours for P — 5 MW, rh = 0.032 kg/s..................................................................................... 150 5.9 Temperature and Mach number contours for P = 5 MW, rh = 0.128 kg/s..................................................................................... 151 5.10 Temperature contours for P = 1 MW, rh = 0.021 kg/s..................... 152 5.11 Temperature contours for P = 10 MW, rh = 0.288 kg/s................... 153 5.12 Collection efficiency vs. average gas temperature. Solid lines show 1-D results (from Rault [44]). Dashed lines show present 2-D results................................................................................... 159 5.13 Overall efficiency vs. thrust for three engine sizes P = 1, 5 , 10 MW..................................................................................... 161 5.14 Performance of the direct solar thruster - specific impulse vs. thrust level for three power levels, P = 1, 5, 10 MW................... 163 6.1 Convergence of the coupled Navier-Stokes and Maxwell equations. Fluids grid 101 x 51, EM grid - (a) 41 x 41, (b) 121 x 121............................................................................................ 169 6.2 Representative flowfield solutions for a floating plasma in a straight duct. P 0 = 1 atm, rh = 4.2 x 10~5 kg/s, Re = 13, Pinc = 2 kW, Pref = 0.16 kW ................................................................. 171 6.3 Distorted electric fields for the straight duct solution in the TMoi2 mode. The electric field is normalized by Eref = 25000 V /m ................................................................................... 172 6.4 Representative solution of the sphere-cylinder configuration. P = 1 atm, rh = 1 x 10~5 kg/s, Re=10, Pinc = 3 kW........................ 175 6.5 Distorted electric field solutions for the TMoiz mode. The electric field lines are normalized by E ref = 25000 V /m ................... 176 6.6 Convergence of coupled Navier-Stokes and Maxwell equations for the bluff body geometry using the TM on mode................................. 178 6.7 Representative flowfield and electric field solutions for the bluff body geometry. P0 = 3 atm, P = 3 kW, Re = 1000.................. 179 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6.8 Variation of coupling efficiency with discharge pressure. Comparison of computations with experiments (from Balaam [101]]............................................................................. 183 6.9 Variation of peak plasma temperature with discharge pressure. Comparison of computations with experiments (from Balaam [101]]............................................................................. 186 6.10 Effect of input power on plasmas in a straight duct. P 0 = 1 atm, m = 4.2 x 10~5 kg/s............................................................................ 189 6.11 Effect of input power on plasmas in the sphere-cylindrical tube. P 0 = 1 atm, rh = 2.6 x 10-a kg/s...................................................... 190 6.12 Effect of input power on bluff body plasmas................................... 191 6.13 Effect of input power on peak plasma temperature for straight duct, sphere-cylinder and bluff body geometries ............. 192 6.14 Effect of gas velocity variation on straight duct plasmas. P0 = 1 atm, P = 3 kW........................................................................ 197 6.15 Effect of gas velocity on peak plasma temperature for the straight duct geometry. P 0 = 1 atm, P = 3 kW ............................. 199 6.16 Effect of gas velocity on sphere-cylinder plasmas. P0 = 1 atm, P = 3 kW........................................................................ 201 6.17 Effect of discharge pressure on straight duct plasmas. P = 3 kW.............................................................................................. 202 6.18 Effect of discharge pressure on sphere-cylinder plasmas................ 204 6.19 Effect of discharge pressure on bluff body plasmas......................... 205 6.20 Effect of discharge pressure on peak plasma temperature for the straight duct, sphere-cylinder and bluff body geometries........................................................................... 206 6.21 Stability Boundary for plasmas sustained in the sphere-cylinder. 208 6.22 Reduction in the threshold power with decrease in the thermal conductivity of helium. k0 is the original thermal conductivity........................................................................................... 211 6.23 Reduction in the threshold power with increase in the electrical conductivity of helium. <r0 is the original electrical conductivity.......................................................................... 212 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. xiii LIST OF TABLES Table Page 5.1 Sizing of the solar thruster for different thrust levels and concentration ratios...................................................................... 136 5.2 Energy balances for the three thruster sizes (P = l, 5 and 10 MW) for different mass flow rates................................................................ 155 5.3 Performance of the direct solar thruster for different sizes, mass flow rates and chamber pressure. (R t is the radius of the nozzle throat)..................................................................................... Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 158 xiv L IS T O F SY M B O LS A ,B ,D Jacobians of E,F and H flux vectors A Avogadro number B Magnetic field including effects of magnetization c/ Speed of light cP Specific heat at constant pressure CR Concentration ratio D Displacement electric field E ,F Dependent inviscid flux vectors E V,F V Dependent viscous flux vectors E, F Electromagnetic flux vectors E Electric field intensity fr Resonant frequency of microwave cavity 9 Acceleration due to gravity G Incident radiation h Planck’s constant H, H Source term vector H Magnetic field intensity h Total black body radiation intensity h Ionization potential Isp Specific impulse Ix Radiation intensity vector J Jacobian of transformation J Electric current density k Boltzmann number, absorption coefficient Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. XV Kc Thermal conductivity of gas Kr Rosseland diffusion coefficient 771 Mass flow rate of gas me Mass of an electron ne Number density of free electrons P Incident power, preconditioning matrix Po Gas Pressure V Poynting power flux vector 9 Electronic charge q Radiative flux vector Q Independent flux vectors for fluid flow field Q Independent flux vectors for electromagnetic field Qv Dependent vector for viscous flow Q Momentum transfer cross-section for atom-electron collisons R Gas constant t Time T Gas temperature Te Electron temperature Th Thrust u,v Velocity in x and y direction u ,v Contravariant velocity in ( and if direction V Swarm velocity of electrons in plasma Ve Nozzle exit velocity of the jet X Axial direction y Radial direction Zo Characteristic impedance of waveguide Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. xvi GREEK SYMBOLS 0 Propagation constant At Time step for time-marching algorithm € Perturbation quantity, electric permittivity Generalized nonorthogonal coordinates A Wavelength P Coefficient of viscosity, magnetic permeability V Frequency of photons Vc Collision frequency 03 Frequency of microwave field P Density Pe Electric charge density O’j e O de Electrical conductivity for a direct current field Tangential direction Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. xvii ACKNOWLEDGMENTS My greatest obligation is to Prof. Merkle, my adviser. He has been a constant source of motivation and guidance; I have learned a lot from him and have enjoyed working with him. I am grateful to him for the trust he placed in me and I hope th at I lived up to it. I would also like to thank Prof. Thynell for his help with the solar problem, Prof. Micci for introducing me to the microwave problem and Profs. Heinsohn, Pritchard and Cimbala for serving on my doctoral committee. I acknowledge the financial support provided by the Air Force Office of Scientific Research (AFOSR) who sponsored this research work. Finally, I owe a great deal to my friends and colleagues - especially Yun Ho Choi, Phil Balaam, K. Srikanth, Denise Cupps and Usha Prabhu. Their support and encouragement made the work seem much easier. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 CHAPTER 1 INTRODUCTION 1.1 Advanced Propulsion Systems During the past decade, considerable attention has been focused on developing advanced propulsion systems. This is primarily a consequence of renewed interest in space-based applications in this country. So far, chemical propulsion has been effec tively applied to satisfying most of the propulsive requirements in space. However, chemical systems are not capable of providing the high exhaust velocities required for some missions. Particularly, for applications such as space station keeping and orbital transfer vehicles, the use of chemical propulsion is not practical. Hence, the need for the so-called advanced propulsion systems is evident. Specific impulse (Isp) is the most commonly employed parameter for comparing propulsion system performance and is defined as the ratio of the rocket thrust to the sea-level weight flow rate of propellant. Chemical propulsion systems, in general, are capable of delivering high thrusts but their performance is limited to Isps in the modest 250-450s range. This limitation is fundamental in th at it is limited by the amount of energy released by the chemical combustion processes occuring in the rocket. Liquid oxygen-liquid hydrogen (LO X —Hz), which is one of the most energetic fuels available today, is capable of providing a maximum Isp of only about 450s. The goal of advanced propulsion systems is to provide thrust at high specific impulses. Various alternative concepts such as ion engines, nuclear rockets, arcjets, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. magnetoplasmadynamic (MPD) thrusters, laser, solar and microwave thrusters [1,2] are being considered. Each of these has a different thrust-specific impulse capabil ity. W ith the exception of nuclear propulsion, these advanced systems generally range from low to moderate in thrust and moderate to high in specific impulse. These increased specific impulses (up to 5000s for ion engines) offer the potential of significant propellant savings. The thrust levels of these devices make them suitable for space-based applications such as orbit transfer, space-station maneuvering and long-term space missions but not for launch vehicles. However, the launch vehicle costs of delivering the upper stage (with its propellant) to low earth orbit (LEO) are significantly reduced by obviating the need for lifting large quantities of chemical propellant into LEO. (In contrast, nuclear rockets offer high thrust and specific im pulse capability and have the potential to replace chemical rockets for many classes of space based missions. However, the lack of public support is a major hurdle in the development of these systems.) Several of the advanced propulsion concepts under consideration are classed under the term electric propulsion [3]. These include lasers, microwaves, arcjets, MPD and ion thrusters. The term electric propulsion is generally defined as the acceleration of gases for propulsion by electric heating and/or by electric and mag netic body forces. Further subclasses may be recognized depending on how the thrust is produced—electrothermal, electrostatic, or electromagnetic. Electrother mal thrusters [4] are similar to chemical rockets in that the hot gases are accelerated in a conventional rocket nozzle. The difference is that the heating is achieved by electric means instead of from chemical combustion of a fuel. Resistojets, arc jets, laser, solar and microwave thermal thrusters all fall under this category. On the other hand, electrostatic and electromagnetic thrusters, such as ion and MPD thrusters, use electric and magnetic fields, respectively, to accelerate the plasma Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3 and produce thrust. Each of the above propulsion systems is at different stages of development. Some of them such as arcjets and ion thrusters have been studied for several years and are well understood. Others such as solar and microwave thrusters are relatively new entrants to the field and are at the proof-of-concept stage. Ultimately, which of these systems will be selected for future space applications will depend on the individual advantages and disadvantages of the systems and how the major issues concerning them are resolved. The present thesis is concerned with assessing the potential of two advanced electrothermal propulsion systems—solar and microwave thrusters. The questions regarding solar and microwave propulsion are: Are these concepts feasible? W hat are the major feasibility issues? How can these systems be optimized? How do they compete with the other concepts being considered? To answer these questions, a fundamental understanding of the physics of the complex radiation-gasdynamic interactions is required. Toward this end, detailed physical models have been de veloped to study the absorption of solar and microwave radiation in gases. These models employ contemporary numerical methods to obtain solutions of the com bined radiation-gasdynamic fields. Along with experiments, the present analytical study should provide enhanced understanding as well as guidelines for future ex periments and analyses. In this chapter, the status of electrothermal thruster technology is summarized. First, the general characteristics of different electrothermal thrusters are discussed. This is followed by a review of various theoretical studies of laser propulsion; laser thermal propulsion has many similarities with solar and microwave propulsion and has been extensively studied both experimentally and theoretically. Solar and mi crowave propulsion are then discussed in the following sections. The overall details Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4 of the present research study are also discussed herein. Finally, a short summary of the present thesis is presented. 1.2 Electrothermal Thrusters Electrothermal thrusters employ electric means to heat a working fluid (often to a plasma state) and then expand it through a standard propulsive nozzle to produce thrust. The primary advantages of electrothermal propulsion over chemical propulsion are: (1) the propellant can be chosen on the basis of its molecular weight rather than its chemical potential (the typical choice is hydrogen or helium because of their low molecular weights); and (2) the peak temperatures can be considerably higher than in conventional chemical rockets. Both of these factors result in higher specific impulses. Different forms of electromagnetic radiation have been considered as energy source [4] ranging from direct currents, radio frequency and microwaves to infrared and visible wavelengths such as laser and solar beams. The primary issue in all these different methods is the coupling of the electromagnetic energy to the thermal (or kinetic) modes of the gas. Because of the wide range of the wavelengths involved, the problems associated with this coupling are unique for each radiation source. 1.2.1 R esistojets Historically, resistojets were the first electrothermal concept to be studied and have been used successfully on several space flights [5]. In resistojets, the fluid is heated by passing it through an electrically heated heat exchanger. Consequently, the maximum temperature to which the fluid can be heated is restricted by the tem perature limitations of the resistance element. Peak temperatures are typically Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5 in the 1500 to 2000K range and corresponding specific impulses are limited to about 300s. Despite their low performance capability, resistojets have the advantage of long lifetimes since they are electrodeless. They are being designed to use waste gases such as CO 2 , Hz0 , CH 4 , O 2 and H 2 as propellant on the Space Station Freedom [6]. 1.2.2 Thermal Arcjets Arcjets employ a dc electric arc between a cathode and an anode to heat the propellant. Electrons produced at the cathode by a combination of thermionic and field emission phenomena are accelerated toward the anode. These high energy electrons transfer their energy to the surrounding propellant by colliding with the heavy gas molecules. Arcjets have been studied extensively both experimentally [7] and theoretically [8] and the underlying physics are well understood. In general, arcjets promise better performance than resistojets but a serious disadvantage is having the electrodes exposed to the plasma [7]. This results in high heat losses (hence poor thermal efficiencies) and erosion of electrodes (hence short lifetimes). However, low power arcjets for satellite station-keeping are being developed for commercial applications [1]. 1.2.3 Laser Thermal Propulsion In view of the limitations of resistojets and arcjets, research attention, in re cent years, has focused on the development of electrodeless high-performance elec trotherm al thrusters. Laser, solar and microwave thrusters are all likely candidates in this class. Of these, laser thermal propulsion has received the greatest atten tion. Here, a high power laser beam is employed to initiate a plasma. The free Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6 electrons in the plasma absorb a major portion of the incident laser beam via in verse bremsstrahlung and transfer this energy to the gas molecules by collisions. The fraction of the incident energy that is absorbed by the plasma is termed the coupling efficiency. The hot plasma transfers this energy by thermal conduction, convection and radiation to the cold propellant gas which is then accelerated in the nozzle. The fraction of the incident energy that is retained in the propellant gas as thermal energy is defined as the thermal efficiency of the thruster. These definitions of coupling and thermal efficiency are employed in the present study of solar sustained and microwave sustained plasmas as well. Clearly, the feasiblity of any of the three technologies is contingent upon achieving high coupling and high thermal efficiencies in the thruster. A series of detailed experimental investigations have been conducted at the University of Tennessee Space Institute (UTSI) [9] in order to gain a quantitative understanding of the complex interactions between the laser beam and the flowing plasma. In these experiments, a 1.5 KW continuous (CW) CO 2 laser was used to sustain an argon plasma. Later experiments have studied plasmas sustained in hydrogen as well. The experiments were also combined with results from a theoretical analysis [10] in order to enhance understanding. Similar work has also been performed at the University of Illinois [11] using a CW CO 2 laser at powers up to 10 kW. Both these studies have demonstrated the feasibility of laser propulsion— the laser energy can be absorbed efficiently in a flowing gas; and the size, shape and location of the high temperature plasma can be controlled optically and fluid dynamically so as to keep the wall temperatures low. These studies also indicated plasma temperatures as high as 20,000K and corresponding specific impulses of up to 2000s. High coupling efficiencies (up to 97%) seem possible; however, thermal efficiencies are considerably lower because of radiation losses. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A comprehensive evaluation of laser propulsion has determined that lasers in the megawatt power range will be required to insure practical application of the concept. Current laser generation capability is restricted to powers in the kilowatt range. However, the development of higher power lasers is underway. Presently, pulsed free electron lasers (FEL) appear to be the most promising way to provide these power levels. The feasiblity of using these lasers to sustain plasmas is currently being investigated at Los Alamos National Laboratory [12]. 1.2.4 Solar Thermal Propulsion Solar thermal propulsion employs concentrated solar energy to heat the working fluid (hydrogen because of its low molecular weight). The primary advantage of solar propulsion is that the energy source is readily available and, unlike laser and microwave propulsion, a man-made source is not needed. However, solar radiation is a dilute form of energy. Moderately sized rockets need large collectors to intercept sufficient quantities of energy and, even after focusing, the maximum energy density remains quite low. This low intensity makes it difficult to couple the electromagnetic energy directly into the thermal modes of the gas. Both indirect and direct absorption concepts [13,14] have been proposed for coupling the solar energy into the working fluid. In indirect absorption schemes, the solar energy is first focused onto a heat exchanger surface which is then used to transfer the energy to the hydrogen gas. The coupling of the energy to the gas is thus readily ensured, albeit at the cost of losing some of the thermodynamic advantage of solar energy. The heat exchanger surface remains hotter than the propellant gas and material considerations restrict the surface temperature to values substantially less than the available thermodynamic maximum (which is 5800K, the effective temperature of the sun). Etheridge [15] and Shoji [16] have conducted Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the preliminary design of an indirectly heated solar thermal rocket that promises an Isp as high as 600s. The design is based on using rhenium coils as the heat exchanger surface, thus allowing metal temperatures of about 2700K. Current efforts in indirect absorption at the Air Force Astronautics Lab (AFAL) are directed toward fabricating a solar thruster for an Integrated Space Flight Test [14]. In the direct absorption thruster, the radiant energy is absorbed directly into the working fluid. Now, the peak temperature occurs in the fluid rather than on a material surface. This allows higher working temperatures and, hence, improved performance over the indirect absorption case. Nevertheless, even if the maximum thermodynamic temperature of 5800K could be reached, hydrogen, as well as most other candidate gases, would remain nearly transparent to solar radiation. The primary potential for direct coupling lies in the addition of small amount of seedant gases such as alkali metal vapors. Hughes Research Lab [17] has generated and sustained plasmas with a direct absorption process in a cesium vapor cell using simulated solar energy. Rault [18] also studied a direct solar radiation receiver using pure potassium vapor. Recently, the possiblility of using a suspension of small sub micron sized particles as the radiation absorber has also been considered [14]. Direct coupling involves a higher level of technological risk than indirect methods. For this reason, current research efforts (as in this thesis) are focused on understanding the complex physics of the radiation-gasdynamic interactions and future experimental research will depend on the results of these studies. 1.2.5 Microwave Propulsion Microwave propulsion is another promising electrodeless electrothermal concept [19]. Generation of microwave discharges have been studied in the USSR [20,21] and Japan [22] primarily in the area of fusion research. In the United States, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. experimental programs at Michigan State University [23-26], Penn State [27-30] and at the Aerojet Company [31] are currently studying microwave absorption for space propulsion applications. In microwave thrusters, the oscillating electromagnetic field accelerates free electrons in the plasma; these accelerated electrons collide with heavy gas atoms, transferring their kinetic energy to them. Again, thrust is obtained by expanding the hot gases in a rocket nozzle. Compared to other electrothermal concepts, microwaves have both advantages and disadvantages. The advantages include the fact that gas absorptivities are high er at microwave frequencies than at optical frequencies, thereby enabling absorption at lower power levels and lower gas temperatures than in laser devices. Further, the longer wavelengths of microwaves allow absorption in a resonant cavity instead of in a single pass through the gas, again enhancing the radiation gasdynamic coupling. Finally, microwave generation is very efficient (up to 85%), and power sources are smaller and less complex than high power lasers or solar concentrators. The primary disadvantage of microwaves is the increased difficulty in “focusing” the microwaves so that the plasma stays in the center of the flow rather than adjacent to the wall, thus preventing high wall temperatures. Several microwave absorption concepts are being considered for space propul sion. These include the resonant cavity, coaxial applicator and waveguide heated plasmas. The resonant cavity concept is the one that has been most extensively studied [23-25,27,28,31]. In this concept, the microwaves are fed by a waveguide to a cylindrical cavity which may be adjusted to tune a particular standing wave mode. The propellant gas flows through a quartz tube positioned along the axis of the microwave cavity. The discharge is typically formed in the propellant gas at a point of maximum electric field intensity on the cavity axis. As mentioned earlier, this arrangement allows multiple passes of the microwaves through the plasma, thus Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 enabling complete absorption of the energy. Indeed, in the above experiments, cou pling efficiencies greater than 99% have been obtained for helium. Measured peak temperatures are in the range of 10,000K and corresponding specific impulses are around 600s. The waveguide applicator is the newest absorption concept and is presently being studied at Penn State [29,30]. In this concept, the absorption chamber is simply a section of a microwave waveguide. Generally, the gas and the microwaves are introduced at the same end of the waveguide section. The microwave plasma generated in the waveguide is inherently unstable and tends to move toward the microwave source. When the gas flow velocity is equal to the propagation velocity of the plasma, the plasma may be held stationary. This situation is similar to combustion waves [32] and laser supported plasmas [9]. The plasma may also be stabilized by means of a bluff body [30] in a practical thruster. Measurements are currently underway to characterize the stabilized plasmas in rectangular and circular waveguide sections. The coaxial microwave thruster concept has been developed at Michigan State University [26]. Here, the plasma is formed and sustained at the tip of a central brass conductor. The primary disadvantages are possible erosion problems and also increased heat losses through the conductor. 1.3 Review of M odelling of Laser Absorption The objective of this research work is to develop comprehensive analytical mod els based on numerical techniques to study direct solar and microwave absorption in flowing gases. Being relatively new concepts, very little theoretical modelling has been attempted so far on either of these problems. However, as mentioned earli er, the related laser problem has been studied theoretically by several researchers. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11 Both one-dimensional and two-dimensional models for the gasdynamic and radia tion fields have been developed and validated with experimental data. It is thus useful to review the theoretical approaches to the laser propulsion problem before we discuss the unique aspects of the solar and microwave propulsion problems. The first analysis of coupling beamed laser energy to a flowing gas was accom plished by Raizer [33]. In his analysis of the subsonic propagation of an air plasma produced by a neodymium laser, he formulated a one-dimensional model and drew an analogy between laser-sustained plasmas and combustion waves. The presence of the laser sustained plasma in a gas gives rise to mechanisms (such as conduction and radiation) that preheat the cold incoming gases to ionization temperatures. The combined effect of these mechanisms in the one-dimensional approximation is to yield a unique propagation speed of the plasma which depends on the local in tensity of the laser beam, the pressure and the optical configuration. The plasma can be held stationary by maintaining gas flow precisely equal to this propagation speed; this gas flow rate thus constitutes an eigenvalue solution to the problem. The plasma wave constitutes a thin region of steep temperature rise and is referred to as a laser supported combustion (LSC) wave. Kemp and his co-workers [34,35] expanded Raizer’s analysis by using realis tic gas properties to study laser powered thrusters. Radiation losses from the hot plasma were included in the analysis by incorporating approximate expressions for the optically thick and thin limits of the radiation. They proposed two absorption schemes—one using pure hydrogen and the other using hydrogen mixed with trace amounts of other gases. Pure hydrogen plasmas showed typical plasma temperatures as high as 20,000K. At these temperatures, inverse bremsstrahlung is the dominant absorption mechanism in the gas. Radiation losses from the hydrogen plasma were Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12 found to be significant resulting in poor thermal efficiencies. In the alternate ab sorption scheme, hydrogen-seedant mixtures were selected so that absorptivity at lower temperatures was enhanced. Gases such as CO 2 and H 2 O provided low tem perature absorption through their molecular vibration levels while trace amounts of cesium were found to be effective at low temperatures because of cesium’s low ionization potential. For these seeded hydrogen plasmas, peak temperatures of about 6000K were predicted. However, experiments with lasers have concentrated on pure hydrogen and argon plasmas [9,11] and the problem of enhancing thermal efficiency has been addressed from within this context. The concept of hydrogen seeded mixtures, nevertheless, finds application in solar sustained plasmas because the maximum allowable temperature here is limited to the sun’s effective surface tem perature (5800K). In a practical laser thruster, the mass flow is set by the choked nozzle throat. Gulati and Merkle [36] studied a practical laser thruster including a convergingdiverging nozzle section. In such a laser device, they argued that two unrelated mechanisms strive to define the mass flow rate—the laser supported plasma admits a unique mass flow rate for a given intensity, whereas the choked nozzle also admits a single flow for given upstream conditions. The only path out of this dilemma is to use a convergent laser beam instead of a collimated beam; typically, the incoming beam is focused by means of a converging lens. The choked nozzle throat is then free to set the mass flow rate while the plasma moves axially along the varying intensity beam until it reaches a local intensity at which the eigenvalue condition precisely matches the choking condition at the throat. As a result of this study [36], an implicit numerical algorithm was developed and applied to study the laser flowfield interaction in a one-dimensional nozzle with a convergent laser beam. These one-dimensional analyses, although useful, rule out the im portant effects Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13 of beam convergence, radial temperature gradients and radial flow velocities. When a convergent laser beam is employed, the plasma is generally confined to the cen tral part of the flow away from the nozzle walls. The radiation gasdynamic field is then clearly two-dimensional. Further, the one-dimensional result of a unique propagation speed is completely misleading. This fact has been demonstrated in experiments [9,11] where laser plasmas have been sustained over a wide range of mass flow rates. The first two-dimensional model of a laser supported plasma was developed by Merkle et al. [37,38]. They solved the 2-D Navier-Stokes equations by using a contemporary time-marching technique. They calculated low tem perature plasmas (peak temperatures of 4000 to 4500K) in a hydrogen flow seeded with cesium. The calculations were performed for a rocket nozzle geometry but were limited to subsonic flows. Temperature dependent absorption coefficients were used but radiative losses from the plasma were not accounted for. In principle, it would have been straightforward to extend their approach to pure hydrogen plasmas (higher temperature plasmas) with radiation losses and choked nozzle flow. Jeng and Keefer [10] developed a two-dimensional numerical model for laser plasmas in pure argon and hydrogen. They included radiation losses from the plasma since this is significant at the higher temperatures encountered (15,000 to 20,000K). For the radiation loss, a “thick-thin” model similar to the one developed by Kemp et al. [35] was used. They carried out extensive comparisons with exper imental data [10]. At pressures greater than 2 atm., they found good agreement in temperature contours (within 10%) and in general features such as plasma size, shape and location. Transmitted laser power was, however, poorly predicted (more than 20% error) in some cases, indicating that the computed absorbed powers were significantly greater than the experimental absorbed powers. Laser beam refraction Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 and diffraction was thought to be responsible for this discrepancy. At low pressures, plasma size was under-predicted. The authors attribute this to the inadequacy of the treatment of optically thick radiation. (Indeed, analysis of radiation loss has been rather simplified in all laser plasma models developed so far.) Hitherto, all models have implicitly assumed the existence of local thermody namic equilibrium (LTE) in laser sustained plasmas. The situation may be pictured as follows. Free electrons in the plasma absorb high energy photons from the inci dent laser beam. These high energy electrons transfer their energy to the heavy gas particles and ions by colliding with them. If this collision frequency is high enough (as it should be at high pressures), all the energy is transferred and equilibrium (or LTE) is said to prevail. However, for the pressures of interest in laser plasmas (1-10 atm ), it is not clear if the LTE assumption is always valid. In the related area of arcjets [39,40], experiments and theory have confirmed the existence of a distinct electron temperature and heavy particle temperature in argon plasmas at atmospheric pressure, thus indicating non-equilibrium. No definite conclusions may be drawn from this since the physics of argon arc plasmas are considerably different from laser plasmas. Zerkle and his co-workers [11] are currently planning similar experiments to determine the extent of non-LTE conditions in laser plasmas. NonLTE considerations are likewise expected to be important in microwave plasmas as well. 1.4 Direct Solar Thermal Propulsion 1.4.1 Review of Modelling of Direct Solar Propulsion Compared to laser thermal propulsion, direct solar thermal propulsion has received very little research attention. Most of the previous work has concentrated Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 15 on indirect solar absorption schemes because they incorporate less technological risk than direct absorption methods. Indirect solar thrusters [14-16] have been designed and are currently undergoing extensive experimental testing. However, indirect methods are restricted to lower gas temperatures and hence lower specific impulses while direct absorption promises higher gas temperatures and, consequently, higher specific impulse. The practical application of direct solar thrusters depends on whether the promised higher performance can indeed be realized in practice. The overall objective of the present study is to determine this. Prior to the present research work, the only models available for solar-sustained plasmas were one-dimensional analyses [41-44] similar to the work of Raizer and others [33-35] for laser plasmas. As with laser plasmas, these analyses result in the prediction of a unique mass flux for a given beam intensity; therefore, the one-dimensional estimates are completely misleading. Further, the effects of beam convergence, radial temperature variation and radial radiation losses cannot be accounted for. Despite these serious limitations, these analyses reveal several inter esting aspects of direct solar absorption in gases. Direct solar absorption differs from laser energy absorption in two respects. First, laser plasmas operate at temperatures in the range of 15,000 to 20,000K. At these temperatures, ionization of hydrogen and other candidate gases is apprecia ble and absorption of laser energy is readily achieved by inverse bremsstrahlung (that is, absorption by free electrons). On the other hand, the maximum operating tem perature of solar plasmas is limited to 5800K, the effective tem perature of the surface of the sun. At these temperatures, hydrogen is virtually transparent to the solar radiation and other means need to be employed to ensure absorption. Further more, alkali metals, with their low ionization potentials, do not ionize appreciably at these temperatures which means that the inverse bremsstrahlung cross-sections Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 of even alkali metals may not be sufficient for solar absorption. The second important difference between laser and solar plasmas is that laser absorption occurs at a fixed wavelength (that is, the wavelength of the incident laser beam), while direct solar absorption occurs over a range of wavelengths (en compassing ultraviolet, visible and infrared regions of the spectrum). Thus, the direct absorbing species should possess an absorption coefficient th at is continuous over the entire solar spectrum. In view of these two points, attention was focused on molecular absorption processes which are appreciable at lower temperatures and continuous over wave length bands. Palmer [41] used pure cesium vapor and Mattick et al. [42,43] used potassium seeded hydrogen in their one-dimensional models of the solar plasma. They determined that these alkali metals were good absorbers of solar radiation. Dimer absorption at low temperatures and photoionization at higher temperatures were identified as the important absorption processes in these alkali metals. The one-dimensional models of Palmer [41], Mattick al. [42,43] and Rault [44] included non-gray radiation loss from the plasma. They predicted solar sustained plasmas between 1000 and 4000K. At these temperatures, the ionization of the absorbing species is quite negligible; the term solar sustained plasma is then a misnomer since the gas does not get heated to a plasma state. Their results also predicted high coupling efficiencies (up to 98 %); however, the high efficiencies are a consequence of the infinitely long absorption chamber assumed in the calculations. Recently, direct solar absorption of the incident solar beam by small submi cron size particles was considered by Thynell and Merkle [45]. They presented the first two-dimensional model to study the volumetric absorption of solar energy. A simplified two-dimensional energy equation (without diffusion terms) was coupled Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 17 with the P -l differential equation for the radiation and solved. The effects of a con verging beam and radiative losses were included for the first time. However, they employed a gray absorption coefficient which is suitable for absorption by particles but not for gaseous absorption. Their results indicated temperatures of 3000 to 3500K and high radial radiation losses if the wall reflectivity was small. Using a wall reflectivity of 1.0 and moderately high gas opacities, they predicted thermal efficiencies of up to 98% at optimal mass flow rates. 1.4.2 Synopsis of Present Solar Research In the present thesis, we address the conceptual feasibility of direct absorption solar propulsion by means of a comprehensive physical model. A schematic diagram of a typical solar propulsion system is shown in Fig. 1.1. The radiation from the sun is incident on the surface of a collecting mirror. This radiation is directed and focused into the absorption chamber of the rocket engine by means of a converging lens system. Hydrogen gas flowing through the absorption chamber absorbs the radiation from the solar beam. The hydrogen gas is then expanded in a supersonic nozzle to produce thrust. In Fig. 1.1, the incoming hydrogen gas is circulated around the nozzle walls thus regeneratively cooling them; this regenerative wallcooling is very important and will be discussed later. Direct absorption of concentrated solar radiation in a gas gives rise to a flowfield that is inherently two-dimensional; the radiation is primarily absorbed by a hot central core while cooler gas flows around it. The aforementioned one-dimensional analyses neglect radial temperature gradients and are therefore too simplified to make a realistic assessment of the capabilities of direct propulsion. The present research [46,47] uses a two-dimensional model which includes the effects of a con vergent beam, radial temperature variation and radial radiation losses. Further, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 18 Working Fluid Concen trated Solar Energy Thrust Transparent Window Figure 1.1 Schematic of a direct absorption solar thermal propulsion system. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 19 non-gray radiation that is characteristic of gaseous radiation is also accounted for [46]. The fluid dynamics are governed by the Navier-Stokes equations with source terms due to radiation in the energy equation. Physical modeling of solar absorption then involves the coupled solution of the Navier-Stokes equations and the equation of the thermal radiation transfer. An accurate modeling of radiation absorption and emission is very important. An exact analysis of the radiation problem is extremely difficult and costly. Further, the non-linear coupling with the fluid dynamics makes the exact analysis impossible even with current supercomputers. Usually, approximate methods such as the thickthin model [35] or the P -l approximation [45] are employed to make the problem tractable. However, the accuracy of these methods with regard to this problem is open to question. In order to estimate this, we have used an exact integral method to solve the radiative transfer equation for a gray gas in an axisymmetric geometry [47]. The results of the 6tudy show that the temperature contours and radiative fluxes from the exact and the approximate methods are in good general agreement. In assessing the potential of direct solar propulsion, we have considered two methods of solar energy absorption—absorption of by alkali-seeded hydrogen [46] and by small size particles [47]. For particulate absorption, the absorbing medium may be assumed to be spectrally independent (i.e., a gray gas); accordingly, the gray gas models described above may be used to evaluate this system. On the other hand, alkali vapors exhibit strong spectral behavior. The absorption coefficients of Na, K and Cs have been calculated as a function of wavelength and tem perature by making use of experimental and theoretical data available in the literature. Using these properties, a two-dimensional model has been developed for solar plasmas sustained in alkali-seeded hydrogen. Spectral and temperature dependence of the absorption coefficient are retained in the analysis. The non-gray radiation is solved Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 20 by a thick-thin model [35] which is calibrated using the P -l differential model. Full Navier-Stokes computations up to the choked nozzle throat are also performed in order to determine thrust and specific impulse of the thruster. The results of these analyses indicate that the direct absorption of solar energy is feasible and plasma temperatures in the range of 2500 to 3500K are predicted. However, the hot gases inside the absorption chamber radiate strongly so that a significant fraction of the absorbed solar energy is re-emitted to the walls. To recapture this energy for propulsive purposes, and to provide cooling to the walls, the working fluid is allowed to cool the nozzle walls regeneratively before it enters the chamber. Some of the re-radiated energy is also transmitted through the transparent window, and, thus, lost. As a result of this regenerative cooling, the absorption in the thruster is essentially a “hybrid” scheme in which energy transfer to the gas is achieved both by direct absorption from the solar beam, and by indirect heating from the walls. The indirect wall heating comes both from the gas re emission and from the remaining unabsorbed energy in the incoming beam after it has traversed the gas. With the regenerative heating, the results show that high thermal efficiencies are possible at optimum flow rates; and the performance of the thruster appears to be competitive with or superior to indirect engines. 1.5 Microwave Thermal Propulsion 1.5.1 Coupled Navier-Stokes and Maxwell Solutions In contrast to the amount of experimental research [23-31] in microwave dis charges, theoretical modelling of the problem has received very little attention. In particular, no two-dimensional models have been developed to describe the microwave-gasdynamic interaction. The major distinction here compared to laser Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 21 and solar plasmas is the necessity for solving the Maxwell equations for the mi crowaves in addition to the Navier-Stokes equations. A one-dimensional model has been developed for propagating waveguide-heated plasmas [29]. The analysis is sim ilar to the propagation of laser waves discussed in a previous section [33] and yields a unique propagation speed for the microwave sustained plasma. No theoretical models have been developed for resonant cavity plasmas prior to the present work. Most of the analyses involving coupled Navier-Stokes and Maxwell solutions have been confined to the related area of radio frequency (RF) discharges [48-50]. While RF discharges cannot be set up in a resonating cavity, they are similar to microwave discharges in that the electric and magnetic field solutions are given by the Maxwell equations. However, most of these analyses have considered on ly one-dimensional Maxwell solutions even though the gasdynamic analysis was two-dimensional [48,49]. Recently, Rhodes and Keefer [50] performed a full twodimensional analysis of the coupled electromagnetic-gasdynamic fields and com pared their results with available experimental data. There are two general approaches to the numerical solution of the Maxwell equations. One approach is to assume a harmonic time dependence of the form eIu;# and then solve the problem in frequency space. The analyses of RF plasmas referred to above have adopted this approach. This method is suitable if the so lution of the first mode is desired, but for higher modes, numerical convergence is seriously impaired. In waveguide-heated discharges (as in the case of the RF dis charges), the dominant mode solution is the only solution of interest and, hence, this convergence limitation does not affect the computation of these discharges. In the case of resonant cavity plasmas, it is possible to set up different standing wave modes simply by adjusting the length of the resonating cavity. In particular, higher order modes are frequently of interest. Then, the solution of the Maxwell equations Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 22 in the frequency domain is not useful. The second approach is to solve the unsteady Maxwell equations by marching in time (similar to time-marching procedures employed for the fluid dynamic equa tions) and averaging over the electromagnetic wave period in order to calculate the RMS values of the fields. This approach may be used effectively for any mode. For this reason, this method is employed in the analysis [51,52]. There are two possible disadvantages to this approach both stemming from the fact that an unsteady so lution is sought. First, the solution procedure adds an additional dimension, time, to an otherwise steady-state problem. The numerical solution of the Maxwell equa tions is however quite economical and the additional overhead is usually not very significant. Second, the unsteady procedure also prevents the implicit coupling of the Maxwell equations with the Navier-Stokes equations. This usually results in slower convergence characteristics of the numerical solution and is discussed later. 1.5.2 Synopsis of Present Microwave Research A comprehensive two-dimensional model has been developed [51,52] to inves tigate the physics and dynamics of microwave gasdynamic interactions in resonant cavity and propagating plasmas. The analysis closely follows companion experi ments being carried out at Penn State [27-28]. Key issues center around the loca tion and size of the absorption volume, absorption efficiency, and selection of the optimal coupling mode, discharge geometry and gas flow conditions for efficient and stable operation. The results of the model are used as an aid in understanding and interpreting experimental results and to extend experimental results to wider regimes (such as scale-up of size and power). In addition, the experiments are being used to validate the computational model. The schematic of the experimental set-up [27,28] which is used as the base Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 23 geometry for the resonant cavity calculations, is shown in Fig. 1.2. Helium gas flows from the top to bottom in the sphere-cylindrical tube while the microwave energy is fed into the cylindrical cavity by means of a coupling probe. In the experiments, the position of this probe and the length of the cavity (by usingithe sliding short) are adjusted in order to tune the cavity for a particular standing wave mode (TM ou and T M 0 1 2 have been used so far). The microwave energy interacts with the free electrons in the helium plasma to maintain the discharge. The microwave plasma discharge is typically formed near the centre of the quartz sphere as indicated. The frequency of the microwave source is 2.45 GHz. Powers up to 3 kW and gas pressures up to 5 atm are obtainable in this apparatus. Recent experiments with resonant cavity plasmas have employed a bluff body plasma stabilization device. Such a device has been employed by M artin et al. [53] for RF plasmas and by Whitehair et al. [24,25] for microwave plasmas. Balaam and Micci [54] are using it in their experiments in order to prevent instability of the plasma at higher powers, which is believed to be caused by the asymmetric location of the coupling probe (see Fig. 1.2). In the present thesis, we have modelled both free-floating and bluff-body sta bilized microwave plasmas. The results of the analyses indicate that the model is qualitatively realistic. Peak temperatures of about 10,000K and coupling efficien cies up to 99% are predicted. Experimental measurements have confirmed these values; moreover, plasma shape, size and location are also reasonably well predict ed. Threshold power, which is the minimum power at which a stable plasma may be sustained, is, however, overpredicted by the calculations. This discrepancy is prob ably due to non-equilibrium effects and is discussed later. More detailed validation will be performed as more detailed measurements become available. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 Sliding Short ' Plasm a Discharge Quartz Container Stationary Short Figure 1.2 Schematic of experimental set-up for resonant cavity plasmas (from Balaam and Micci [28]). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 25 1.6 Overview of the Thesis The rest of the thesis is organized with the details of the analyses presented in the first half and the results of the analyses presented in the second half. According ly, Chapters 2, 3 and 4 are devoted to the physical modelling of the electromagnetic, thermal radiation and fluid dynamic fields respectively. Chapters 5 and 6 present the results of the direct solar absorption and the microwave absorption problems respectively. Chapter 7 presents an overall assessment of the solar and microwave propulsion systems based on these analyses. The electromagnetic field solutions—necessary for the absorption of microwave energy—are discussed in Chapter 2. This involves the numerical solution of the Maxwell equations using time dependent techniques. First, the fundamental physics represented by the Maxwell equations are discussed. The numerical algorithm em ployed in the present work is then presented. Next, the method of characteristics procedure used to specify the boundary conditions on the Maxwell equations is ex plained. This method is superior to conventional boundary procedures and allows the specification of incident power instead of just the absorbed power. Also, the de termination of the electrical properties of equilibrium plasmas is discussed. Finally, typical Maxwell solutions for resonant cavity microwave plasmas are presented. Chapter 3 deals with the solution to the thermal radiative transfer problem. The Maxwell equations discussed in Chapter 2 are not useful in describing radiation over infrared and visible frequencies because of the extremely small length and time scales involved. In these regimes, the equation of thermal radiation transfer must be solved. For the thermal radiation solutions, two approximate methods and one exact method are discussed. Demonstrative solutions of the direct solar absorption problem using the three methods are presented and compared for accuracy. Further, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 26 the determination of radiative properties for alkali vapors are also dealt with. Chapter 4 describes the numerical methods employed in solving the NavierStokes equations for the fluid dynamics. Algorithms for the standard Navier-Stokes as well as special procedures for low Mach number computations are discussed. Flowfield solutions for the sphere-cylinder configuration and the bluff body con figuration employed in the microwave discharges are presented to demonstrate the performance of these algorithms as well as to get an initial feel for the nature of the flowfield in these geometries. Chapter 5 presents the results of the feasibility study of direct solar absorption in hydrogen gas. Two absorption schemes—one using small particles and the other using alkali vapors as the absorbing species—are studied. For the former, a gray gas analysis is performed, while for the latter, a non-gray analysis with a band model for the spectral dependence of radiative properties is used. Performance characteristics of three different thruster sizes are studied for a range of flow conditions. The investigation of microwave absorption in gases is the subject of Chapter 6. Coupled Maxwell and Navier-Stokes solutions are presented for different flow configurations. The effect of various parameters such as gas pressure, gas velocity and microwave power on the nature of the microwave discharges are investigated. The results from the computations are also compared with available experimental data. Chapter 7 draws various conclusions from the above analyses of direct solar absorption and microwave plasmas. The areas where the present analyses may be extended are identified and recommendations are offered for future analytical work. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 27 CHAPTER 2 SOLUTION OF THE MAXWELL EQUATIONS 2.1 Electromagnetic Field Theory Electromagnetic radiation fields are governed by the generalized Maxwell equa tions. While the analysis of these equations is suitable for the study of radiation in the microwave frequency range (viz., 0.01 to 10 GHz.), it is not so convenient to use these equations in the case of radiation in the optical and infrared frequencies. This is because the high frequencies and small wavelengths of optical and thermal radia tion introduces very small time and length scales which make the Maxwell solutions extremely tedious. Thus, in these regimes, the common practice is to employ the laws of geometric optics (or the equation of thermal radiation transfer) which are an asymptotic form of the Maxwell equations in the limit of very small wavelengths. In this chapter, we will focus on the numerical solution of the general form of the Maxwell equations which is applicable to the microwave absorption problem in res onant cavities and waveguides. The analysis of thermal radiation transfer which is applicable to solar absorption and plasma reradiation is dealt with in the next chapter. 2.1.1 The Maxwell Equations The Maxwell equations are written as the set of equations [55]: V-J + ^ = °. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.1) Here, J = the electric current density, H = the magnetic field intensity, E = the electric field intensity, D = the displacement electric field, B = the magnetic field intensity including the effects of magnetization and pe = the medium electric charge density. Equation (2.1) is simply a statement of the conservation of charge. Equation (2.2) is a generalization of Ampere’s Law and says that a time changing electric field produces a magnetic field. Corresponding to this, Eqn. (2.3), which is Faraday’s Law, says that a time varying magnetic field produces an electric field. Equations (2.4) and (2.5) are divergence relations that may be derived from the previous three equations. In these equations, B and D may be likenedto flow vectors.The electriccharge density pe acts asa flow source for the vector D while there are no flow sources associated with B. Various im portant relations exist between the electric and magnetic field vec tors which lead to some simplifications of the above equation set. D = eE, (2.6) B = pH. (2.7) where e is the permittivity of the medium and p is the permeability of the medium. In addition, the electric current density is usually related to the electric field inten sity by Ohm’s Law, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 29 J = <rE. (2.8) Here, a is the electrical conductivity of the medium and is generally a function of number of free electrons in the medium. Further discussion of this im portant property is reserved for a later section. Equation (2.8) suggests that the current density is parallel to and in phase with the electric field. This is not generally valid since alternating fields introduce phase lags. These effects are also discussed later. Using Eqns. (2.6) to (2.8), the Maxwell equation set may be expressed in terms of two fundamental vectors - the electric field intensity E and the magnetic field intensity H . In addition, the electric charge density may be assumed to be negligibly small (a valid assumption for weakly ionized plasmas). Then, the Maxwell equations for a general co-ordinate system may be written as: ATI V x H = <rE + 6 ^ - , at VxE = Boundary conditions the boundary represents , (2.9) (2.10) V • E = 0, (2.11) V *H = 0. (2.12) are needed to close the aboveset ofequations.Usually, an interface between twomedia. If this interface is a perfect conductor, no electric or magnetic field can exist within it and the following boundary conditions [55] apply. E t = 0, (2.13) H n = 0. (2.14) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 30 Here, E t is the component of the electric field tangential to the interface and H n is the normal component of the magnetic field. Since we are concerned with microwave fields in waveguides and resonant cavities, the boundaries are typically metallic and the above conditions are specified. Apart from these, the only other boundary condition of interest is the specification of incoming microwave power. This is discussed in a later section. 2.1.2 E lec tro m ag n e tic Force a n d E n erg y The microwave sustained plasma may be visualized as a volume enclosing a region of ionized gas. Free electrons in this volume are accelerated by the electric field drawing energy from it. These electrons then constitute an electric current which interacts with the electric and magnetic fields. A charged particle in an electric field experiences a force. Assuming the plasma to possess an electric charge density of /t>e, this force is given by F e = peE. Similarly, a charged particle moving in a magnetic field experiences a force which is given by F m = PeEind = Pe{y X B) = J X B, where v is the swarm velocity of the charged particles. The total electromagnet ic force (or Lorentz Force) exerted on the plasma by the combined electric and magnetic fields is then given by Fern = j x B + peE. (2.15) The power flow through any closed surface enclosing the plasma may be in terpreted from the surface integration of a power density flux vector V = E x H , Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 31 known as the Poynting vector [55,56]. The integral form of the theorem of Poynting may be expressed as / T , d s = I i / t H B + E •D ]dV + J ^ 3 >EdV. (2.16) The term on the left is the ingoing power flux as indicated earlier. The first term on the right denotes, at any instant, the time rate of increase of the total stored electromagnetic energy within the volume V enclosed by S. The | E • D represents the stored energy per unit volume associated with the electric field while the | H - B represents the stored energy associated with the magnetic field. The last term is the total power dissipated (or generated) within V at any instant. Substituting J = <rE identifies this last term as an Ohmic power loss term. Equation (2.16) has an important significance to the problem of microwave absorption. In the microwave plasma, the electrons accelerated by the incident electromagnetic (EM) field impart their energy to the heavy gas molecules and ions by collisions. This is then the power loss term in Eqn. (2.16). This power has to be replenished by the ingoing power flux to maintain the plasma. Generally, we are interested in Eqn. (2.16) in the time-averaged sense since the time scales of fluid motion are much slower than the time scales of the electromagnetic field. For a steady harmonic EM field, it may be readily shown that the time average of the rate of change of stored electromagnetic energy (that is, the first term on the right of Eqn. 2.16) is zero. Thus, the time-averaged power flux entering the closed surface must equal the power dissipated within the interior volume. If the incident microwave power exceeds the power absorbed in the plasma, then the fraction of the incident power that is not absorbed is reflected back such that the net power flux across the surface stays equal to the net power absorbed. In the event that Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 32 there is no power absorption (when the plasma is extinguished), the net ingoing time-averaged power flux has to be zero. This means that the incident microwave power is completely reflected back to the source. In the presence of the plasma, depending on the dynamics of the electron gas molecule collisions, a certain fraction of the incident power is absorbed and the rest is reflected back. This percent of the incident power that is absorbed is called the coupling efficiency of the absorption process. 2.1.3 T h e E lec tro m ag n e tic W ave E q u a tio n The Maxwell equations represent the general solution of any electromagnet ic field. For instance, these equations may be employed along with appropriate boundary conditions to solve for the electric and magnetic fields inside a waveguide or resonating cavity. The electromagnetic wave equation is useful in providing a physical picture of electromagnetic waves in such devices. Taking the curl of Eqns. (2.9) and (2.10), and using the following vector identity V x V x E = V(V • E) - V 2E we can derive the wave equations for E and H in a general co-ordinate system. 0E 02E VB = ^ j r + V p - . W B - + (2.1T) There are three fundamental modes of electromagnetic wave propagation transverse electromagnetic (TEM), transverse electric (TE) and transverse magnetic (TM). In a given situation, we are interested in only one of these modes. As a demonstration of the wave equation solution, we consider a TEM wave. TEM waves Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 33 have no field component in the direction of propagation. In other words, the electric and magnetic fields are purely transverse to the direction of propagation. If TEM waves are propagated along the x-axis in a rectangular co-ordinate frame, E x = 0 and H x = 0. Also, assuming that the z-axis is oriented parallel to the magnetic field and the y-axis parallel to the electric field, we get H y = 0 and E t = 0. Finally, we assume that the fields are uniform in the transverse plane and that the medium is lossless (<r = 0). Then, Eqn. (2.17) becomes dZEy a* = dEy d’ H , SB, - S T = ^ -gr- <2-18) The solution of the electric field may be written as : E y == JB0+ eiw(*/c,“ f) + E 0~ e - iu,(x/ci+t). (2.18a) This solution represents two sets of waves - forward travelling and backward trav elling. The velocityof these waves is given by cj = -j==, the velocity of light. In practice, the forward travelling wave represents the incident wave and the backward travelling wave the reflected wave component as discussed in the previous section. In a lossless semi-infinite waveguide, the reflected component is zero, in which case, there is only the forward travelling wave. 2.1.4 Microwave Waveguides TEM waves cannot be supported by a single conductor waveguide (they can only be carried by a two-conductor coaxial waveguide). TE waves and TM waves are then commonly used in waveguides. In a TE wave, the electric field is transverse to the propagation direction. If the direction of propagation is along the x-axis, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 34 E x = 0. For TM waves, the magnetic field is transverse to the propagation direction and H x = 0. In general all the other electric and magnetic field components will be non-zero. As an example, we consider the propagation of TM modes in a cylindrical waveguide. Using a cylidrical co-ordinate system, xyQ (the axial, radial and tan gential directions respectively), H x — 0 and the E x satisfies the wave equation w ritten in cylindrical co-ordinates. For a lossless, homogeneous medium and a har monic time variation (e,UJ<), the wave equation may be solved using separation of variables [57]. E x = E oJ n{kcy)cos{n9)ei(u>t- 0x\ (2.19) Here, J n{kcy) represents the nth order Bessel function where n is the tangential mode number. Applying the boundary condition at y=R, that is at the waveguide wall, we get J„(kcR) = 0. Jn(kcR) is an oscillating function with an infinite number of roots. Depending on the particular radial mode desired, the appropriate root is selected. For the m th radial mode, the m th root is selected. If the m th root of J„(kcR) = 0 is qnm, then the wave number for the radial mode, kc — qnm/RFinally, /? is called the propagation constant of the waveguide and for a given mode, it is given by (3 = y/u>2[i€ — kc2. Since has to be real for propagation of microwaves in the waveguide, this relation puts a lower limit on the value of the frequency of the microwaves. Thus, for a given mode, wave propagation occurs only if the signal frequency is greater than this cut-off frequency. The cut-off frequency is a function of the mode number, waveguide dimensions and the material properties of the medium filling the waveguide. From the above discussion, it is clear that Eqn. (2.19) represents an infinite set Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 35 of solutions since m and n can take on values from 0 to oo. Each solution represents a mode of transm itting electromagnetic energy down the circular wave-guide. The particular mode is designated as T M nm where the subscripts m and n indicate the number of half-wave variations of the field components in the y and 9 directions respectively. The other field components may now be readily written out in terms of the known E x component. From Eqns. (2.9) and (2.10), we get the following relations. ifi d E x k 2 dy ’ Hy = iw e 9 E X 09 ’ -s e dEt Ilkc2 09 ' CM Eg — iwe 0 E X k 2 Or ( 2 .20 ) For waveguide-heated plasmas, the choice of propagation mode influences the size, location and shape of the plasma and hence the absorption characteristics of the device. Figure 2.1 shows some typical mode patterns for TM waves in a circular waveguide. Consider the TMom modes which have no 9 variation of the fields. This implies th at the fields are axially symmetric and, further, that the Eg and H v modes are zero. The only non-zero components of the axisymmetric TM electromagnetic field are E r , E y and Hg. Similarly, the axisymmetric TE modes (TEom) have non zero Eg, Hx and H y components. In this research study, we focus our attention on these axisymmetric modes. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction TEu Wave Type Field distributions in cross-sectional plane, a t plane of maximum transverse fields Field distributions along guide ™01 ^ ^ ^ ^ ^ ^ lis trib u tio n s F below along th is plane •/-r— •*•A***'•ill *i» tf jU jf c a U B @# riiiliiife i:: « 0 - -« k-.-- - Hz’ Hr ’ V Er« E0 ™02 ^^^^^U stributions below along this plane |aTC3j i s m K T^)K T^l>: Ez*, Er*, E„, 0 Hr ’, H0a Ez« Er ’ H0 u i s i S Field components present ™n TE01 i Ez’ Er ’ He V « r' Ee prohibited without perm ission. « Figure 2.1 Field distribution for TM and TE modes in a circular waveguide (from Gandhi [57]). CO o> 37 2.1.5 Microwave Resonators A microwave resonator is a cavity in which electromagnetic waves are enclosed. The wall of the cavity is usually a good conducting material. When the frequency of the microwaves is tuned to the resonant frequency of the cavity, the waves form a standing wave pattern. In a lossless medium, the waves oscillate endlessly without losing energy and the energy oscillates between the electric and magnetic fields in the cavity. In the presence of a lossy substance such as a plasma, the energy is dissipated by ohmic heating as the microwaves bounce to and fro. The cavity has to be replenished with this lost energy and this is done by coupling a waveguide carrying microwave energy to the cavity (as shown in Fig. 1.2). The resonant cavity may be viewed as a waveguide with shorted ends and may be analyzed in the same manner as waveguides are. For a given cavity, a unique resonant frequency exists for each mode. These resonant modes of the cavity are denoted as T E nmp and T M nmp where m, n and p are the number of half-wave variations in the tangential, radial and axial directions respectively (for a cylindrical cavity). The resonant frequency corresponding to a particular mode may be calculated as follows. Since the propagation constant (32 = u r2fie — kc2, f vfr* + r 27r^//Ie Substituting for /3 and kc, we get Alternately, for a given frequency, the cavity length when resonance occurs may be determined (for a given diameter of cavity). This corresponds to the physical Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 38 situation in experiments with resonant cavity plasmas [27,28], where the cavity length is adjusted by means of a sliding short in order to tune the cavity to a particular resonant mode. (Figure 1.2 shows this schematically.) In analyzing wave propagation in waveguides only the forward travelling wave solution is considered. For resonant cavities, since waves are bouncing back and forth, both forward and backward travelling wave solutions must be included. The axial component of the electric field for TM waves in a lossless cavity filled with a homogeneous medium is given by E x = E 0 Jn(kcy)coa(n6 )(eik*x + e~ik*x)/2 or, E x = E 0 J n(kcy)coa(nd)coa(kxx). (2.22) Here, kx has been used instead of /3 and is given by kx = pir/l. The other relations given by Eqn. (2.20) still hold and the remaining components may thus be deter mined. These solutions represent the general T M nmp mode. Some typical mode solutions are shown in Fig. 2.2. Again, the TMomp modes have an axisymmetric field pattern and E x, E y and Hg are the only non-zero components. Just as in the case of waveguides, the mode patterns determine the location, size and shape of the plasma. Experiments have employed the TM on and the T M qu modes frequently [27,28]. Consequently, these are the modes used in the calculations as well. 2.2 Time-Dependent Solutions o f the Maxwell Equations 2.2.1 Introduction to Maxwell Solution Procedures The exact solutions given in the previous section may be obtained only for ho mogeneous media. In the presence of a high temperature plasma, the gas properties Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (a) TM01() mode (b) TM0U mode (c) TMm mode (Useful for rectangular waveguide mounted fre quency meters) (d) TM012 mode mode (i) TM131 mode Figure 2.2 Field distribution for TM waves in a cylindrical resonant cavity (from Gandhi [57]). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 40 are strong functions of temperature and, consequently, they vary considerably over the volume of interest. Furthermore, for microwave plasmas, it has been shown [58] that the modes are generally not pure but are of hybrid character. In addition to the modes that are related to those of the empty (lossless) cavity, additional resonances that do not exist in the empty cavity are found to occur. In general, the solution of the electric field then has to be obtained by the superposition of several mode solutions. Computational solutions are then the only means of obtaining general Maxwell solutions. Different approaches have been used ranging from variational and integral approaches to finite element and finite difference methods [59,60]. The logical choice here is finite differences because of their robustness and efficiency and also because it is the method that we employ to solve the fluid equations as well. Previous research with radio-frequency plasmas [48-50] have assumed a har monic time dependence of the electromagnetic fields and thus removed the timedependence of the Maxwell equations. As mentioned earlier, this method is efficient only if the solution to the lowest order mode is desired. The radio frequency analy ses have considered only the lowest order mode and, for that reason, these analyses have also been restricted to the one-dimensional Maxwell equations. For higher order modes, the convergence of these schemes deteriorate [59]. To offset this, conjugate gradient techniques employing a bi-harmonic operator [61,62] have been developed. Results from using these approaches seem encouraging; however, the complexity of the formulation somewhat offsets these advantages. Time-dependent methods to solve the unsteady Maxwell equations were de veloped in the mid-sixties [63]; these methods have since been popularly termed as FDTD (Finite Difference Time Domain) schemes. They have been applied to solving a wide variety of microwave problems [64-66] such as waveguide propaga tion and scattering by complex obstacles. However, they have not been used in Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 41 the analysis of microwave discharges in gases. In fact, there are very few instances in the literature of the analysis of coupled thermal and electrical fields using the FDTD procedure. One example is the computation of the EM fields and tempera tures within a microwave irradiated eye [66]. The powers employed here were of the order of a few milliwatts and temperature variations were only a few Kelvin. The electromagnetic field was then uncoupled from the temperature field which means that the EM field was first computed and then the energy equation was solved using the computed field distribution. The advantage of using the time-dependent procedure is that higher order modes may be computed with the same ease as lower order modes. As mentioned earlier, most of the experiments with resonant cavity plasmas have employed the TM U12 and the TM on modes [27,28]. Clearly, it is desirable that the physical model developed have the capability of handling such higher order modes. A second advan tage with this procedure is that hybrid modes and mixed modes are automatically included in the analysis. Previous analyses [58] have demonstrated the existence of hybrid modes in the presence of a lossy plasma. In principle, it is possible to ex actly replicate the electromagnetic fields in the cavity. However, since we choose to solve the axisymmetric Maxwell equations, non-axisymmetric modes are excluded. Inclusion of these modes is straightforward but that would mean extension of the analysis to three dimensions. The main disadvantage with the time-dependent approach is the need to solve an unsteady electric field for an otherwise steady-state problem. In addition, the electric field solution is non-linearly coupled to the solution of the fluids equations. Then the question arises as to how this coupling between an unsteady electric field and a steady flow field may be achieved. Physically, the time scales of fluid motion are much larger than the time-scales of the oscillating EM field. Thus, the fluid sees Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 42 the electric field as a time-averaged field. This RMS value of the oscillating field is determined by averaging over an EM wave period. The time-averaged Lorentz forces and power loss terms are calculated and used as the source terms in the momentum and energy equations respectively. Exactly how this is achieved is discussed along with the fluid dynamic model in a later chapter. A standard FDTD procedure for solving the unsteady Maxwell equations is described by Yee [63]. The Yee algorithm is a second-order accurate (in space and time) explicit Crank-Nicolson [67,68] scheme. Recently, the hyperbolic nature of the Maxwell equations has led to the development of CFD-based solution methods. Shankar et al. [69] have used a characteristic-based finite volume upwind algorithm similar to Euler equation solvers [67] while Goorjian [70] has developed an implicit approximate factorization scheme [67,68]. While these new algorithms possess some attractive properties, there does not seem to be any specific advantage to using characteristic splitting of the electromagnetic flux vectors except in specifying the microwave inflow boundary condition. Secondly, explicit schemes are easier to code than implicit schemes. For these reasons, we have chosen to employ the explicit procedure of Yee along with a characteristic-based boundary condition procedure. We describe the algorithm in the following sections while the boundary condition procedure is discussed in a later section. 2.2.2 Vector Representation of the Maxwell Equations As mentioned earlier, our interest is confined to axially symmetric modes (such as the TM U12 and the TMun ) where there are no tangential variations of the fields. For this axisymmetric case, the Maxwell equations simplify further such that Eg, H x and Hy are the only non-zero components for TE waves and H$, E x and E y are the only non-zero components for TM waves. Accordingly, the axisymmetric Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 43 Maxwell equations may be written in the following compact vector form. 9Q 9E 9F , W + lfc + a r * - ' 2-23> For TM waves, / Hey\ Q = \ E ty , \ EyyJ f ( l / f i ) E yy \ E=[ 0 , \(iA )# « y / / - ( l / n ) E xy \ F = - (1 /e)Hey , \ 0 / ( (1 ! n ) E x \ H = -(<r/e)Exy . \-l<T/e)Eyy J We may write the Jacobians A = 9 E / 9 Q and B = 9F /9Q . If the eigenvalues of A and B are real then the system of equations represented by Eqn. (2.17) is hyperbolic in character. Then the equations may be likened to the Euler equations of fluid mechanics. In addition, these eigenvalues represent the characteristic speeds or the propagation speeds of the electromagnetic waves. For instance, the Jacobian A for TM waves is : 0 0 0 The eigenvalues of A are 0, +^== and —-j=*. The last two terms are once again the velocity of light in the medium. 2.2.3 FD TD Solution Procedure From Eqn. (2.23), it is clear that the magnetic field depends only on spatial derivatives of the electric field while the electric field, in turn, depends only on spatial derivatives of the magnetic field. Accordingly, the flux vectors may be split up as Q = Qh + Qe, E — Eh + E e and F = Fk + Fe. The vectors are defined Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 44 so that those with subscript h contain only magnetic field components while those with subscript e contain only electric field components. Qe = E xy E P= U A )Hey 0 / —( l / y ) E xy \ Fe - I 0 \ 0 \ Ft = I -(1 , 0 / ■ Then, Eqn. (2.23) may now be written as: f + f % + + f = o, * l i t =s ' <2,24i) The Yee algorithm may now be represented by the following two step procedure. Taking central differences in time with Eqn. (2.24a) at time level n and Eqn. (2.24b) at time level n + | , we get, Q r M r ! = - 2A < ( f + f ) ” ~ * (1 + ^ A t)(Q ”+1 - Q") = _ 2 A < [ ( ^ + ^ - ) (2.25a) j* - f X ' - H n] (2.256) Here, we note that H = —<r/eQe. The first step (Eqn. 2.25a) advances the magnetic field solution in time while the second step (Eqn. 2.25b) updates the electric field solution. The spatial derivatives are handled using central differences which are second-order accurate. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 45 Since the above algorithm is an explicit time-marching scheme, stability of the algorithm becomes an issue. In general, the stability characteristics of a hyperbolic system of equations axe well understood [67,68]. Stability considerations reveal a conditional stability criterion due to the CFL condition [67]. The CFL number is traditionally defined as ci8 t / 8 £, where c/ is the characteristic speed of the system (or the velocity of light), 8 t is the time step size and 8 £ is the minimum grid cell dimension. The CFL condition limits the maximum allowable CFL to 1 for stability of the time-marching procedure. In other words, there is a limit to the size of the time step that may be used in the calculations. Further, the finer the grid the more severe this limitation. This point has to be borne in mind especially when local grid stretching is desired. A second issue arises from the odd-even splitting character of central difference schemes. In the extreme case, two completely different solutions may result from the scheme - one from the odd numbered grid lines and one from the even numbered grid lines. Upwind differencing, staggered gridding and the use of artificial smoothing operators are the common ways of eliminating this problem. In the present work, as in the original Yee formulation [63], staggered grids have been used; hence, no additional artificial smoothing is required. This means that the different electric and magnetic field components are evaluated at appropriate points within the grid cell. Figure 2.3 shows these locations within a typical grid cell. 2.3 Characteristic Boundary Conditions In the experiments with resonant cavity plasmas, microwaves are supplied to the cavity by means of a coupling probe as shown in Fig. 1.2 [27]. This location of the coupling probe is asymmetric with respect to the cavity axis. In order to keep the axisymmetric formulation of the computational model, the microwave input is Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 46 i,jH 1/2 i + 1/2, j + 1/2 Figure 2.3 Staggered grid showing locations of electric and magnetic field components within a grid cell. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 47 modelled by a coaxial waveguide located in the bottom of the cavity. This is depicted in Fig. 2.4, where i2,- and R 0 stand for the inner and outer radii of the coaxial waveguide respectively. Note that the microwaves are introduced into the cavity at the bottom end since the coupling probe in the experiments is located near this end. The coaxial waveguide, which carries microwaves in the TEM mode, may be included in the computational model simply by extending the computational domain to include a portion of the waveguide. At the power input end of the waveguide, the incident microwave power has to be specified as a boundary condition. The RF plasma studies [48-50] have invoked an awkward iterative procedure to specify the E-field at the boundary. W ith this procedure, it is only possible to con trol the amount of power absorbed in the plasma and it is impossible to distinguish incident power from reflected power. In the present study, a physically realistic boundary condition procedure involving characteristic variables was developed in order to enable the specification of incident power. Using characteristic boundary conditions, the power reflected from the cavity may be determined from the solu tion of the Maxwell equations in the cavity. This procedure thus precisely describes the interaction between the microwaves and the plasma in the cavity. Clearly, the situation in waveguide-heated plasmas is similar and the same procedure may be used for this case as well. The characteristic form of the Maxwell equations may be derived in the same way as for the Euler equations [67]. To obtain the characteristics in the x-direction, we start with Eqn. (2.23) and multiply through by the x-directional modal matrix (the matrix of left eigenvectors of the Jacobian A). <2-26> where M x~l is the modal matrix for the ®-direction and for axisymmetric TM waves Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 48 Figure 2.4 Geometry of microwave cavity showing the inlet waveguide and the co-ordinate system. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A in Eqn. (2.26) is the diagonal matrix containing the eigenvalues of A and Qx is the characteristic vector in the x-direction that is defined as Qx = M x~l Q and is given by: / qz = E xy i(H<i+ \ f i E»)y The characteristic vector, Qx, contains the forward and backward travelling waves or the Riemann invariants of the system. For the co-ordinate system shown in Fig. 2.4, the backward travelling wave that corresponds to the negative eigenvalue (-ci) is the incident microwave. This is given by q~ = \y(He - E y/ Z 0) and has to be specified as a boundary condition. The forward travelling wave represented by the quantity q+ = \y(He + E v/ Z 0) corresponds to the positive eigenvalue (+cj) and is the reflected component. In both these relations for the characteristics, Z 0 = y / E and is the characteristic impedance of the waveguide. This definition of Z 0 is analogous to the impedance in transmission line theory [57]. The co-axial waveguide that supplies microwaves to the cavity supports TEM waves. The Maxwell equations for TEM waves reduce to the wave equations given in Eqn. (2.18). The solution of these equations are the superposition of the for ward travelling and backward travelling wave solutions given in Eqn (2.18a). From transmission line theory [57], we can show that the forward characteristic is equal to the forward travelling electric field component and the backward characteristic is the backward travelling electric field component. Z 0 q+ = E 0 +eiui, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 Z 0 q~ = E„~eiut. Then, if E y = E yewt and H$ = Hee'ut, we get: , £>+ = 5»(6 + VZ„)Z„, - E,IZ„)Z,. The instantaneous incident power is given by the incident Poynting Vector which is: Zo Usually, it is the average incident power that we wish to specify rather than the instantaneous power. To do this, the above equation is averaged over one wave period. For a co-axial waveguide, this is given by: Tinc = -£-(E 0 - f l n { R 0 / R i ). "0 (2.27) Thus, incident power may be specified by fixing E 0~ at the boundary according to Eqn (2.27). It is clear that instead of specifying the electric field at the bound ary, the correct boundary condition involves the specification of the appropriate characteristic variable, which is a combination of electric field and magnetic field components. The reflected power is likewise given by an equation similar to Eqn. (2.27) except that E 0+ replaces E 0~ . This reflected power is thus the outgoing characteristic and is calculated from the solution of the Maxwell equations in the cavity. The power absorbed by the plasma is the difference between the specified in cident power and the calculated reflected power. Once the reflected power is known, the coupling efficiency which is defined as the ratio of absorbed power to incident power may be determined. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 51 2.4 Electrical Conductivity o f the Plasm a 2.4.1 Electrical Conductivity Expressions The electrical conductivity of the plasma plays a crucial role in determining the magnitude of the ohmic power dissipation of the microwave field; thus, it controls the magnitude of power absorbed by the plasma. Simple relations estimating the electrical conductivity in direct current and high frequency fields are derived by Jahn [71]. The motion of electrons under the influence of an electric field may be described by the relation: ^ 7 = — E - vcv at me (2.28) which may be regarded as the equation of motion of a ficititious average electron. In Eqn. (2.28), q is the electronic charge, m e the mass of an electron, v is the swarm velocity of the electrons and vc the collision frequency of the electron with heavy particles.Then, the first term on the right is the acceleration of the electron by the applied electric field and the second term is the deceleration due to collisions. The electric conduction current J arises due to the motion of electrons and is defined as J = n egv, where n e is the electron number density in the plasma. Equation (28) may now be re-written as: ^ = — B_^J. at me (2.29) For a steady electric field, the time derivative in the above equation drops out when steady-state is achieved and we get for the conduction current j = m euc = ^E . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (2.30) 52 The <TdC is called the dc (direct current) conductivity. This expression shows that the electrical conductivity is proportional to the number of free electrons in the plasma and inversely proportional to the collision frequency. For a high frequency field, the current and electric fields may be assmed to vary harmonically with time (e,wt). Substituting this in Eqn. (2.29), we still get J = <rE except now the electrical conductivity is complex. This means there is a phase difference between the driving electric field and the current field due to the inertial response of the particles and the frequency of collisions. The real and imaginary parts of this complex ac electrical conductivity are given as [71]: n eq2 r <Tr _ vc ~ m e to;2 + i/2. neq2 r a> m e lw2 + v*\' The dc limit is recovered as u; tends to 0. This equation for the ac conductivity is useful only if we are solving the electromagnetic field in frequency space. However, since we are solving the time-dependent equations, the complex ac conductivity cannot be directly introduced into the equations. In order to include the effects of a high frequency field, it is necessary to use Eqn. (2.29) instead of the Ohm’s Law relationship (Eqn. 2.8) assumed so far. This would allow phase differences between the electric field and the conduction current to be included automatically. In other words, Eqns. (2.9) to (2.12) are not valid and we have to go back to Eqns. (2.1) to (2.5) to derive the modified form of the Maxwell equations th at is valid for alternating current fields. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 53 2.4.2 Modified Maxwell Equations Introducing the dc conductivity as a variable, Eqn. (2.29) may be rewritten as: = (o-dCE - J)i/C. Note th at the above equation reduces to Eqn. (2.8) when J j = 0, which is the dc limit. We employ this equation in addition to the standard Maxwell equations to represent the microwave field. Accordingly, the new set of governing equations for axisymmetric TM modes may be written in our conventional vector form. 8 Q IF + / H 9y \ E vV Q= E xy J vy \ J xy / f - ( l / f i ) E xy \ 0 F = -(1 /e)Hgy , 0 \ 0 / 8 E 8 F + <*=*’ (2.32) / ( l / l i ) E yy \ (1 M H e y E = 0 0 0 / 6 = -(1 l v ) E x - ( 1/<0 J I,y \ -(l/e)/*y (<TdcEjy —Jy)vcy {(ffdcEx ~ Jx)vcy) This enhanced set of equations is more general than Eqn. (2.23) in th at the phases of the electric field and the conduction current are not assumed to be the same. In particular, for the alternating current field, two additional equations are used to determine the conduction current instead of the algebraic Ohm’s law relationship th at is valid for direct current fields. The solution of Eqn. (2.32) proceeds in much the same manner as described in Section 2.2.3 with the current components being treated on the same time level and at the same spatial locations as the corresponding electric field components. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 54 2.4.3 Equilibrium Ionization The dc conductivity given by Eqn. (2.30) depends on the electron number density of the plasma. Assuming local thermodynamic equilibrium between the electrons and the heavy particles, the ionization of gaseous helium (which is the working fluid of interest here) is described by the celebrated Saha Equation [72]. For singly ionized helium, the Saha equation gives: 7ie7i; /2ttm ekT e \ = 2Qe ( — . _ exp(—I p/ k T e). 3/ 2 (2.33) Te is the electron temperature which is approximated as being equal to the heavy particle temperature. The ion number density n; is equal to the electron number density n e for a neutral plasma. The neutral atom number density n„ may be assumed to be equal to the number density of the helium gas because the degree of ionization is very small. The n„ is simply a function of the local pressure and temperature of the discharge. Substituting in all the constants in Eqn. (2.33), the Saha equation may be used to determine the equilibrium electron density as a function of the local temperature and pressure of the plasma. n e = 9.829 x 10lon Q1/2Te3/4esp(-1.425 x 10s/T e). For a temperature of 10,000K and pressure of 1 atm , n e = 5.3 x 1O10 m ~3. The degree of ionization, a = n e/n„ = 7 x 10-5 which is only few thousandths of a percent. The assumptions of negligible charge density that we made in Section 2.1.1 are then certainly valid. 2.4.4 Collision Frequency Electrons in the plasma collide with neutral gas atoms, ions and other electrons and, in doing so, transfer their kinetic energy. At the ionization levels that we are Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 55 interested in, electron-neutral atom collisions predominate and all other interactions may be neglected. The collision frequency of a single electron may be related to the atomic-scale momentum-transfer cross-sections [71]. According to this, uc = n aQve. (2.34) The Q is the momentum-transfer cross-section for atom-electron collisions. From the calculations of Frost and Phelps [73], Q is found to be reasonably constant over the range of electron energies of interest. We have used Q = 5.2 x 10" 20 m2. The ve is the mean scalar speed of the electrons relative to the heavy particles. This should not be confused with the swarm velocity v of the electrons. Assuming a Maxwellian distribution of velocities, fm Z Ve = \ ------ • V TTTTle Thus, the collision frequency may be determined as a function of local temperature and pressure of the plasma. Note that vc is proportional to the number of heavy particles - that is, the greater the density of heavy particles, the more collisions an electron will suffer. Thus, at higher pressures when the number density increases, the collision frequency increases. Increase in temperature results in two effects. It decreases the atomic number density but increases the translational energy of the electrons. The total effect is to decrease collision frequency as Te *. 2.4.5 Electrical Conductivity Calculations Using Eqns. (2.30), (2.33) and (2.34), the dc conductivity may be determined as a function of temperature and pressure. Figure 2.5 shows the conductivity plotted against electron temperature for pressures ranging from 0.01 to 100 atm . It may be seen th at the conductivity increases exponentially with temperature indicating the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 56 0.01 0.1 1 atm 100 10 O E "t! b J3 -3 1 -6 0 5000 10000 15000 20000 TEMPERATURE, K Figure 2.5 Direct current conductivity of a helium plasma at different temperatures and pressures. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 57 ‘avalanche’ effect characteristic of ionization. The conductivities become apprecia ble (about 1 mho/m) only at temperatures over about 8000K. Figure 2.5 also shows th at the conductivity increases as the pressure decreases; this is a consequence of smaller collision frequencies at lower pressures. These results agree closely with the theoretical results of Lick and Emmons [74]. We have formulated the Maxwell equations for the microwave field in terms of the dc conductivity. However, the physical parameter of interest is the oc conduc tivity. For this reason, we consider the behaviour of the ac conductivity though we do not use this parameter explicitly in our calculations. The real part of the ac conductivity (Eqn. 2.31) is shown as a function of temperature and pressure in Fig. 2.6 for a field frequency of 2.45 Ghz. We can see that pressure has a very different effect on ac conductivity compared to dc conductivity. Increasing pressure from 0.01 atm increases the conductivity. At a pressure of about 1 atm the conductivity reaches a maximum and then starts decreasing. At low pressures, the collision fre quency is much smaller than the field frequency. Increasing pressures increases the electron number density and hence the conductivity of the gas. At higher pressures, the collision frequency becomes comparable to the field frequency and increasing pressure starts reducing the conductivity just as in the dc case. This behaviour has an im portant effect on the absorption characteristics of microwave plasmas and is discussed in Chapter 6. 2.5 Demonstrative Solutions of the Maxwell Equations Various essential features of the Maxwell equation solution procedure have been discussed. In this section, some sample results are presented. To decouple the analysis from the fluid dynamic solution, the plasma profile was specified arbitrarily. In other words, the electrical conductivity of the gas was specified as a known Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 58 _ 1 atm • ' 100 0.1 0.01 atm -3 -6 5000 10000 15000 20000 TEMPERATURE, K Figure 2.6 Real part of ac conductivity for helium plasmas at different temperatures and pressures. Field frequency = 2.45 GHz. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 59 function of spatial location. The assumed plasma profile in the cavity is shown in Fig. 2.7. The plasma core is assumed to have a uniform electrical conductivity while, surrounding the plasma, the electrical conductivity is taken to be zero. At the surface of the plasma, the electrical conductivity grows exponentially from zero at the outside to its assigned value on the inside. The electrical conductivity of the core was varied for different cases in order to evaluate the effect of the magnitude of the lossy term on the Maxwell solutions for the cavity. In order to compare the solutions for the no-loss case with the exact solution, a case was computed with zero conductivity. The result from the exact solution for the T M 0 1 2 mode is shown in Fig. 2.8. The axial and radial electric field components as well as the intensity of the E-field squared (|E |2, or the electromagnetic energy field) are shown. This no-loss case was computed numerically as well. Figure 2.9 shows the temporal development of the electric field at a given location in the cavity for this case. The electric field does not settle into a harmonic ‘steady’ state. This is because it is difficult to precisely match the length of the cavity to the desired resonant length. Thus, the waves in the cavity do not form an exact standing wave pattern. Figure 2.10 shows the corresponding spatial variation of the axial and radial electric fields and the E-field squared. In spite of the difficulty in generating an exact standing wave pattern in the cavity, comparison of Figs. 2.8 and 2.10 show th at the overall shapes of the field contours are in good qualitative agreement. The only disagreement is in the region where microwaves enter the cavity. This is to be expected since the exact solution does not include the coupling of the coaxial waveguide to the cavity. From the no-loss results of the TMoiz mode, three points of maximum intensity - two at the ends and one at the center - are seen to occur along the cavity axis. Gas breakdown resulting in the formation of the plasma can initiate at any one of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2.7 Assumed plasma profile at the center of the cylindrical cavity. Plotted lines represent log10<Tdc. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 61 (b )£ , 0 (c) |E |5 Figure 2.8 Analytical solution for the T M ^ mode in a no-loss cylindrical cavity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 62 O oo K ?" -----i----------------- 1------------------------------------11----------------- 1----------------- 1----------------- 1----- o o© n H------ u— ’.0 0 0 .6 2 5 1 .2 5 0 1 .8 7 5 2 .5 0 0 3 .1 2 5 3 .7 5 0 4 .3 7 5 x10'*-8 TIME, s Figure 2.9 Temporal development of the electric field at a point in the no-loss cavity. Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission. “1 5 .0 0 0 M 0 (C) |E| Figure 2.10 Spatial distribution of computed electric fields inside the no-loss cavity. TMui2 mode. / = 2.45 GHz. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 64 these three nodes. Experiments have demonstrated that the plasma is most stable at the middle node. The plasma profile in Fig. 2.7 assumes the plasma to be located near the middle node as well. Figure 2.11 shows the time development of the electric field for a case where the electrical conductivity of the plasma core was 1 mho/m. In contrast to Fig. 2.9, the field appears settle into a harmonic steady state after an initial transient response. This result suggests that increased absorption of the microwave energy by the plasma changes the resonance characteristics of the cavity. However, there is still evidence of a low frequency response, which may be due to the low level of power absorption by the plasma. Figure 2.12 shows the spatial distribution of the electric field components and the E-field squared. Note the major difference in the shape of these contours compared with the no-loss case. In the region of the plasma, there are steep gradients of the E-field as the energy is lost to the plasma. Most of the energy is absorbed by the surface of the plasma and, within the plasma, the E-field is nearly zero. This inability of the E-field to penetrate the plasma is called the skin effect. Outside the plasma, the basic shape of the TM 012 mode is seen to be retained. A second case was computed with the electrical conductivity of the core taken as 100 m ho/m . Temporal and spatial solutions for this increased conductivity are shown in Figs. 2.13 and 2.14. The temporal development in Fig. 2.13 shows a slightly faster relaxation of the solution to a harmonic steady state compared to the previous case (Fig. 2.11). The spatial contours in Fig. 2.14 are qualitatively similar to the solution in Fig. 2.12 except that the plasma is bigger. Also, the electric fields remote from the plasma are much weaker because of the stronger absorption of the plasma. Still, the overall mode shape is evident in these regions. For the above case, the incident microwave power was 3 kW. Of this, 1.17 kW Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. - 7 .0 0 0 000 i .6 2 5 i 1 .2 5 0 T 1 .8 7 5 i 2 .5 0 0 i 3 .1 2 5 j - ........................j .. 3 .7 5 0 4 .3 7 5 | 5 .0 0 0 x10"-8 TIME, s Figure 2.11 Temporal development of the electric field at a point in the lossy cavity. Peak value of <Tdc = 1 mho/m. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 66 (a) E , (b)B, (c) |E|> Figure 2.12 Spatial distribution of the electric field in the lossy cavity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -1.400 - .9 3 3 D8/:+^^A .0 0 0 £25 1 .2 5 0 1 .8 7 5 2 .5 0 0 3 .1 2 5 3 .7 5 0 4 .3 7 5 ?000 x10"-8 TIME, s Figure 2.13 Temporal development of the electric field at a point in the lossy cavity. Peak value of a jc = 100 m hos/m . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 68 (a) E x (b) Ey (c) |E |2 Figure 2.14 Spatial distribution of the electric field in the lossy cavity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 69 was absorbed by the plasma and 1.83 kW was reflected back to the source. This gives a coupling efficiency of 39%. In general, the cavity length may be adjusted to minimize the reflected power component and thus increase the coupling efficiency. This is because the presence of the plasma detunes the cavity by changing the resonant cavity length. Adjusting the cavity length to the new resonant length (by trial and error), allows multiple passes of the microwaves through the plasma and hence better absorption. This is, indeed, the fundamental advantage of the resonant cavity arrangement for microwave absorption. Figure 2.15 shows the change in the coupling efficiency with slight variations of cavity length. The length in the figure has been normalized with respect to the no-loss resonant length. The coupling efficiency is seen to peak when the cavity length is about 5-10% shorter than the no-loss resonant length. These results show that the time-dependent Maxwell algorithm can be used ef fectively to solve for the microwave field in resonant cavities. The EM field solution is, however, coupled to the solution of the gas dynamics by the strong tem perature dependencies of the electrical conductivity and the collision frequency. The tem per ature field is in turn coupled to the microwave field because of the microwave heat source. Thus, the Maxwell solution procedure has to be coupled iteratively to the fluid dynamic model for the determination of the flow and electric fields. Results from the fully coupled model are discussed in Chapter 6. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 70 COUPLING EFFICIENCY o On LO_ h- OlO lO_ CM .92 .94 .96 .98 1.00 1.02 L /L 0 Figure 2.15 Variation of coupling efficiency with cavity length. No-loss resonant length of the cavity, L 0 = 0.148 m. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 71 CHAPTER 3 SOLUTION OF THERMAL RADIATION TR A NSFER 3.1 Direct Solar Absorption in Gases Solar radiation is electromagnetic radiation in the range of relatively small wavelengths. Figure 3.1 shows Planck’s black body intensity distribution at the sun’s surface temperature of 5800 K. It can be seen that most of the radiation is over the visible wavelength range (0.4 to 0.8 fim) with small percentages in the ultraviolet and infrared portions of the spectrum. As mentioned earlier, these small wavelengths make the Maxwell equation analysis impractical, thus necessitating the need to employ the equation of thermal radiation transfer. There are two aspects to the solar absorption problem - one, the direct absorption of solar radiation by the gas, and two, the reradiation from the hot gas to the surrounding cold gas and the chamber walls. This chapter is devoted to both these aspects. The analysis is primarily directed toward the solar absorption problem but the same methods may be used to describe the radiation from microwave plasmas as well. In a solar thruster (Fig. 1.1), concentrated solar radiation is employed to heat a gas which is subsequently expanded in a supersonic nozzle to generate thrust. The hot gases in the absorption chamber radiate strongly - emitting a significant portion of the energy. A portion of this reradiation is absorbed back by the gas and the remainder is incident on the chamber walls and the inlet window. Further, energy in the beam that is not directly absorbed by the gas is also incident on the chamber walls. The energy incident on the walls may be used to regeneratively Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 72 INTENSITY (H7(m*.m)) 3 2 1 0 0.0 0.2 0 .4 0.6 0.8 W AVELENGTH 1.0 1.2 1 .4 1.6 (mm) Figure 3.1 Black body intensity at the surface temperature of the sun. Tsun = 5800 K. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 73 heat the working fluid before it enters the chamber thereby utilizing the energy productively. Hydrogen gas is an attractive choice of working fluid because of its low molec ular weight. However, the low intensities of solar radiation introduce difficulties in the coupling of the solar energy to the gas. In fact, hydrogen is virtually trans parent to radiation in the visible and the infrared regions at the temperatures of interest. Different methods of direct absorption such as seeding with alkali metal vapors and suspending a swarm of small particles are being considered. Of these, alkali vapors seem more promising because it is relatively easy to seed hydrogen gas with small quantities of alkali metals. Introduction of small particles such as soot is more difficult. Further, particles may increase the weight of the propulsion system and may also cause wall erosion problems. The radiative properties of alkali vapors are important to ensure the absorption of solar radiation. The main absorption mechanism in alkali vapours is transitions between electronic levels of atoms and dimers - that is, bound-bound transitions [4144]. The alkali metal dimers (eg., Na.2 , C 3 2 ) contribute to most of the absorption and are known to be strong in the visible and infrared regions. Molecular absorption may operate at temperatures as low as room temperature up to the tem perature at which the molecule dissociates. At elevated temperatures, dimer absorption is likely to be less important and other absorption mechanisms need to be considered. Photoionization, or bound-free transitions, is a likely candidate given the low ionization potentials of alkali metal atoms. In fact, for alkali vapours, high-energy photons (in the ultraviolet region) are effective even at low temperatures. Photoionization bands in the visible region become important around 2000K. A third mechanism is inverse bremsstrahlung (or free-free transitions) which involves the interaction of electrons with ions or neutral particles. However, alkali metals (even heavier Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 74 atoms like Cs) ionize appreciably only at temperatures above 4000K; therefore, in verse bremsstrahlung is not likely to be an important absorption channel for solar plasmas. The following picture then emerges for the solar absorption process. The gas is first heated regeneratively by wall heating to about 1000K. This elevated tempera ture is necessary to maintain the alkali metal in the vapor state. The hot inlet gas is then continuously heated up by the solar energy due to strong dimer absorption of the alkali molecules. At about 2000 to 3000K, the photoionization bands start becoming im portant. Due to dissociation of molecules, the dimer absorption drops down to values comparable with the photoionization absorption. At still higher tem peratures (about 4000K), inverse bremsstrahlung may contribute to the absorption process, especially in the infrared. In this chapter, we first consider the determination of the radiative properties of different alkali metals such as Na, K and Cs. The absorption coefficients of alkali metals are found to be strongly dependent on frequency and temperature. In con trast, the absorption coefficients of small particles are approximately uniform with wavelength. A gray gas assumption is therefore appropriate for particle absorption but not for alkali vapors. Thus, the radiative properties of the absorbing species influence the form and solution of the radiation equation. Accordingly, models of different levels of complexity have to be employed to represent each of these radi ation problems. Following the discussion on radiative properties, the equation of radiative transfer is discussed and various methods of solution of both the gray as well as the non-gray problem are described. Finally, some representative solutions from these models are given. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 75 3.2 Radiative Properties Modelling 3.2.1 Dimer and Atomic Line Absorption The internal energy of a gas is contained in discrete vibrational, rotational and electronic energy states of its atoms or molecules. The absorption of a photon can cause a transition of the atoms or molecules to a higher energy state. Because only discrete energy states are involved in these transitions, only photons of certain energies can be absorbed (or emitted). If the energies of the upper and lower discrete levels are E j and 22,-, then only photons of energy E p = E j —Ei can be absorbed. The energy of a photon is related to its frequency through Bohr’s condition Ep = Ej —Ei — hu (3*1) where h is the Planck’s constant. Consequently, the discrete bound-bound transitions result in the absorption of photons of very definite frequencies resulting in line spectra. The absorption coefficient of the gas is related to the rate of these transitions by means of the Einstein coefficient B , j , which is defined as the probability per unit time and volume of a transition occuring from state i to state j per unit incident energy flux. For gaseous atoms, this relation is [75]: kv = BijNihu{ 1 - e~h^ kT) = k*( 1 - e~hu/ kT). (3.2) The factor (1 —e~hv/ kT) accounts for induced emission. This is because the ab sorption of energy by a gas stimulates some of the gas atoms or molecules to emit energy. Equation (3.2) predicts that very little energy would be absorbed from the inci dent spectrum by any given absorption line since only photons of a single frequency Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 76 are absorbed by that line. However, other effects cause the line to be broadened and, consequently, each line has a finite width about v. Some of the important line-broadening mechanisms are natural broadening, Doppler broadening, collision broadening and Stark broadening [75]. In diatomic molecules, in addition to electronic transitions, two additional modes of motion, which do not occur in atoms, are possible. First, the molecule can rotate as a whole about an axis passing through the center of gravity and per pendicular to the line joining the nuclei (the internuclear axis) and, second, the atoms of the molecule can vibrate relative to each other along the internuclear axis. The transitions between vibrational and rotational energy states in addition to the electronic transitions result in very closely spaced spectral lines. In fact, in cer tain spectral regions, the absorption lines are so closely spaced that the individual lines cannot be resolved. The lines often overlap and merge to form continuous absorption spectra or band spectra. All alkali metal vapours have two prominent dimer bands. These are the socalled X —►A (due to E+ —» n u transitions) and X —> B (due to E+ —* E+ transitions) bands. In addition to these, alkali metals possess a strong atomic resonance doublet. Typically, the X —* B dimer band lieB on the blue side of the resonance line and the X —* A band lies on its red side. It has been observed that the resonance doublet, though atomic in origin, is broadened following the selection rules of the molecule [76]. Accordingly, the quasi-static line-broadening theory can be applied to explain the line-widths of the doublet. In the present work, the resonance lines have been represented by Lorentzians as described in Ref. [76] and are dependent on the square of the atomic number Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 77 concentrations. As for the dimer bands, they are calculated by using A?(jimer — 0"absJVd imer (3.3) where <rabs is the absorption cross-section of the dimer and JVdimer is the dimer con centration. Data for the absorption cross-sections were obtained from the literature [77-80]. The dimer concentration was calculated by considering the equilibrium reaction for the molecule: 2[Aik] & [Alk2]. (3.4) [Alk2] is the number density of alkali molecules and is given by: [Alk2] = JVdimer = {Aik}' RT exp[-A G °{T)/R T] P0A \ (3.5) where A is the Avogadro number, R is the universal gas constant and the exponen tial term containing the Gibbs free energy term G° is the equilibrium constant, K p. Combining Eqns. (3.3), (3.4) and (3.5), we get RT ^dim er — °"abs [-A/fe] PoA exp[-A G°(T)/RT) (3.6) The data for the equilibrium constant K p as a function of tem perature was obtained from the J A N N A F tables [81]. It should also be noted that induced emission has been neglected. This assumption is justified for the wavelengths and temperatures of interest. The maximum error introduced is about 10 percent but is significantly less over most of the wavength-temperature domain. 3.2.2 Photoionization The absorption of a photon with sufficient energy by a gaseous atom can cause ionization or the ejection of an electron. The resulting ion and electron are free Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 78 to take on any kinetic energy; hence, the transition is bound-free. The bound-free absorption coefficient is a continuous function of the photon frequency v as long as the photon energy hu is sufficiently large to cause ionization. Alkali metals have relatively low ionization potentials (especially the heavier atoms like Cs) because of their free valence electrons in the outermost shell. The absorption coefficient is calculated using fcphoto = c abaN exp~I^ kiT ( 3 .7 ) where I p is the ionization potential, N is the number density of the alkali atoms and <rabs is the frequency-dependent absorption cross-section for bound-free absorption. The cross-sections for ground-state photoionization were obtained from Ref. and the excited-state cross-sections from Ref. [82] [83]. 3.2.3 Inverse Bremsstrahlung The inverse bremsstrahlung mechanism involves the interaction of electrons with ions and neutral particles. A free electron, passing near an ion, interacts with its electric field resulting in a free-free transition. The electron can absorb a photon thereby gaining kinetic energy. Since the initial and final free energies can have any value, a continuous absorption spectrum is produced. By the same token, inverse bremsstrahlung can also be produced if an electron passes very close to a neutral atom since there is an electric field very close to an atom. Since the mechanism requires electrons, it can only operate at temperatures where free electrons are present - that is, after the gas has sufficiently ionized. In the case of alkali metals, inverse bremsstrahlung becomes an effective absorber at temperatures of 4000 to 6000K. Also, free-free cross-sections vary as 1/t/3 and will be more im portant in the infrared and far infrared regions of the spectrum. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 79 The free-free absorption coefficient, according to reference [84], is given as 64flr4 m N e 10z4 exp(—*j) brem " 3\/3 2h'cv* ^ where xi — Ip/k^T and Iv is the ionization potential. N is the neutral atom concentration. Figure 3.2 shows the bremsstrahlung absorption coefficient for Na, K and Cs at a temperature of 4000K. It is appreciable in the longer wavelengths (longer than 1.2 pm) for Cs but is somewhat less important for K and Na at the temperatures of interest here. Since only a small proportion of the solar energy (about 10%) is contained in this wavelength region, the bremsstrahlung contribution is not im portant and is neglected in the present work. 3.2.4 Complete Absorption Coefficient Based on Eqns. (3.6) and (3.7), the absorption coefficient of three alkali metals - Na, K and Cs - were constructed. Figures 3.3, 3.4 and 3.5 show the absorp tion coefficients of these alkali metals carried in hydrogen gas (hydrogen does not contribute to the absorption) as a function of wavelength for three temperatures of 1500K, 2500K and 4000K. The properties are calculated using a total pressure of 3 atm and an alkali mole fraction of 0.05 (i.e., a partial pressure of 0.15 atm). These values are representative of working conditions in the absorption chamber of a solar thruster. Consider the absorption spectrum of Na as an example (Fig. 3.3). Here, the photoionization bands extend up to 0.4 pm, while the dimer bands extend from 0.4 to 0.8 pm; the resonance line occurs at 0.57 pm. Ground-state photoionization is active even at low temperatures but excited-state cross-sections become appreciable only at higher temperatures. Dimer bands occur on either side of the resonance line; the dimer cross-sections are very high at 1500K and drop Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 80 A B S O R P T IO N C O E F F IC IE N T (m - 1 ) 10 10 Na 10 -2 10 - 3 10 - 4 10 - 5 0 .5 2.0 W A V E L E N G T H [a m ) Figure 3.2 Inverse bremsstrahlung absorption coefficient for Na, K and Cs. T = 4000 K, P0 = 3 atm, yca — 0.05. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 81 10 C O E F F IC IE N T ( m - *) 10 - A B S O R P T IO N 10- 15 0 0 K 2S00K 4000K 10 -1 0.2 W A V E L E N G T H (m m ) Figure 3.3 Absorption coefficient of Na vs. wavelength. jP„ = 3 atm, = 0.05. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 82 10 1500K C O E F F IC IE N T ( m ~ ') 10 2500K A B S O R P T IO N 4000K 10 0.6 W A V E L E N G T H (m m ) Figure 3.4 Absorption coefficient of K vs. wavelength. P0 — 3 atm, t / k = 0.05. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 83 10- ABSORPTION COEFFICIENT (m --1500K - 2500K — ItOOOK 0. 0 0. 2 0.4 0.6 0. 0 1.0 1.2 WAVELENGTH (pm) Figure 3.5 Absorption coefficient of Cs vs. wavelength. P0 = 3 atm , yca = 0.05. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 .4 84 rapidly with temperature. The absorption spectrum of K and Cs are very similar but they cover a larger wavelength range. The absorption coefficients are dependent on the number density of the al kali atom. Increasing this number density will increase the absorption coefficient (fcphoto increases as N while fcdimer increases as N 2). The alkali concentration may be increased by increasing the partial pressure of the alkali. There are two ways of doing this - by increasing the percentage of seed in the hydrogen and by increasing the total pressure. Figure 3.6 shows the mean absorption coefficient of a H i/C s mixture as a function of temperature. To study the effect of increasing pressure and mole fraction, curves are shown for two total pressures (1 atm and 10 atm) and two alkali mole fractions (®cs = 0.05 and 0.1). It is clear that at lower chamber pressures, higher alkali mass fractions are required to maintain high absorptivity levels; however, higher percentages of alkali metal in hydrogen will increase the average molecular weight of the working fluid. At higher chamber pressures, high absorptivity levels may be attained with relatively small alkali mole fractions; how ever, gas pressures are likely to be limited by structural considerations of the thrust chamber. Clearly, there is a design trade-off involved between performance and structural considerations; we will discuss this aspect in greater detail in Chapter 5. Since the strong resonance lines and the dimer bands of Na, K and Cs occur over different wavelength intervals, it would seem advantageous to use a mixture of the three alkali metals as seed. This would provide a high absorptivity level over a larger portion of the solar spectrum. Figure 3.7 shows the absorption coefficient of such a mixture. The gas was chosen as a mixture of 95% hydrogen and 5% Na, K and Cs (i.e., 1.67 % of each species). The total pressure was again taken to be 3 atm. The figure shows that the dimer bands now extend from 0.4 to 1.4 (im while the photoionization bands cover the ultraviolet and a portion of the visible spectrum. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 85 10 10 3 x c s = 0 .0 5 10 2 p = l.a tm A B S O R P T IO N C O E F F IC IE N T ( m ” 1) p = lO .a tm i c« = 0.1 p = l.a tm x c « = 0 .0 5 10 I 2000 2500 3000 3500 4000 4500 5000 Figure 3.6 Mean absorption coefficient of Cs vs. temperature for different pressures and alkali mole fractions. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 86 ABSORPTION COEFFICIENT I (m lOOOK 2000K AOOOK 0 .0 0.2 0.4 0.6 0.0 1.0 1.2 WAVELENGTH ( ym) Figure 3.7 Absorption coefficient of a 5% mixture of Na, K and Cs (1.67 % of each alkali) in hydrogen. P 0 = 3 atm. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1.4 87 Note that again the dimer absorption coefficients decrease with temperature while the photoionization absorption coefficients increase with temperature. 3.2.5 Band M odel for Non-Gray Gases It is evident that the absorption coefficient of alkali vapors is strongly wavelength-dependent. Traditional gray gas models which use a mean absorption coefficient will be inaccurate. In order to consider the spectral behavior of the ab sorption coefficient in the radiation model, a band model is used. The band model approximates the absorption coefficient as a constant within each band: N B fcA(A) = X > ( * ) (3.9) »=i where N B is the number of bands of width AA; and fc; is the average absorption coefficient for the ith band and is zero outside the band. The different bands used for the alkali mixture shown in Fig. 3.7 are listed in the Appendix. 3.3 Equation o f Radiative Transfer 3.3.1 The Thermal Energy Equation The thermal radiation problem is strongly coupled to the thermal modes of the gas. High temperatures and high absorptivities result in greater radiative transfer which, in turn, governs the temperature of the gas. The radiation transfer term appears in the thermal energy equation as a gradient of the radiative flux (similar to the thermal diffusion term which is the gradient of the conductive heat flux). In the analysis of beamed energy interaction, however, it is customary to divide the radiation problem into two parts - (1) the direct absorption from the incident Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 88 beam, and (2) the reradiation by the gas. Accordingly, the energy equation has two source terms due to radiation energy transfer: pCp% + PCpu J ^ = V - K CV T + q"’ - V • qR (3.10) where V is written in cylindrical co-ordinates and K c is the thermal conductivity. The q1" term is the radiant heat absorption from the solar beam per unit volume and V •qR is the gradient of the reradiative flux vector (which includes the effects of emission and re-absorption by the gas). Equation (3.10) can be solved independent of the fluid dynamic equations if the mass flux pu is known. In a solar thruster, the mass flux is typically specified by the choked nozzle throat and its determination involves a full fluid dynamic analysis of the absorption chamber and the thruster nozzle. However, to simplify the analysis, we decouple the thermal energy equation from the momentum equations by assuming quasi one-dimensional flow and constant pressure. The problem of specifying the correct value for the mass flux is addressed in Chapter 5. 3.3.2 T h e R a d ia tiv e T ran sfe r E q u a tio n In the absence of scattering, radiation travelling along a path is attenuated by absorption and is enhanced by both spontaneous and induced emission. The equation of radiative transport [75,84-87] along a path s is given as dh{d^ X) + A) = A), (3.11) where k \ is the wavelength (or frequency) dependent absorption coefficient, I \ is the spectral intensity of the radiation and I\,\ is the black-body function. The radiative flux vector qR is defined by q j = / h fld w , J 4tt (3.12) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 89 where each component of is obtained by taking the dot product with the unit vector in that direction. The conservation of energy requires the gradient of the radiative flux vector, which is one of the radiation source terms in the energy equation. This may be derived by integrating equation (3.11) over all solid angles. V - j J = * »(< »*»- G » ) . (3.13) where G \ is the incident radiation that is given by G \ = f I xdw. (3.14) JiTC W ith those definitions in mind, we proceed to formulate the equation of ra diative transfer for a cylindrical medium. The intensity vector I \ is a function of x, y, 6, <j>and A; here, x and y are the axial and radial co-ordinates, 6 is the polar angle and <f> is the difference in the azimuthal angles of the intensity vector and the radial position vector. Starting from Eqn. (3.11), we arrive at the following equation of transfer in axisymmetric cylindrical co-ordinates [85]. 'iToxT + ' 2 Ioy T “ “y TiJ 0<p + = kxhx (3-15) where I \ = l\(x,y,0,<f>, A) and l\ = cos 9, h = sin 6 cos <j>and Is = sin 6 sin <f>. An exact solution to Eqn. (3.15) may be obtained by representing the angular distribution of intensity by a series of spherical harmonics called the P -n approxi mation [86]. However, this requires the solution of a very large number of coupled differential equations and is not practical. Various methods have been proposed to solve Eqn. (3.15) [87]. Methods such as the zonal method, the discretized intensity method or the Monte Carlo methods are accurate; however, their computation al requirements skyrocket when spectrally-dependent and temperature-dependent properties have to be included. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 90 As mentioned earlier, the two radiative source terms in the energy equation are treated independently. The direct absorption term is treated by using a com prehensive ray-tracing technique [37], which is described in the following section. The gradient of the reradiative flux may be treated in a variety of ways. We look at three such methods - (1) an exact integral technique [88], which is time consuming and is suitable only for gray problems with constant absorption co-efficients; (2) the P -l Differential method [87], which has the advantage of being a differential equa tion model and can handle spatial variation of the absorption co-efficient; however, the method becomes cumbersome for non-gray gases since a system of differential equations must be solved; and (3) the ‘thick-thin’ model [34,35] which involves no differential equations and is hence very economical for non-gray solutions. We use the exact technique primarily as a benchmark with which the other approximate models may be validated. The radiative transfer term is solved using the P -l model for particle absorption and the thick-thin model for the alkali vapor absorption. 3.3.3 Direct Solar Absorption For direct solar absorption from the incident beam, a ray-tracing technique following the laws of geometric optics is employed to split the incoming solar beam into a series of individual rays. Due to the finite size of the sun’s image, the rays need to be traced from the surface of the collector. Figure 3.8 shows this ray model with 180 constant area rays. The larger the number of individual rays, the closer we are to an exact representation of the incoming beam and the more accurately we can represent the absorption pattern in the working fluid. The shape of the beam is set by the focal spot size and the /-num ber of the converging lens (or collector mirror). Downstream of the focal point, a ‘reflected’ ray, corresponding to radiation originating in the oppposite half of the beam, enters the computational domain. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 91 'tx'KKimau SSSSSWCMI ^ -sX ? 5 ^ S S ! ^> ^5«O S «! £S*J®WSSS&; t>\v.x\:scvxv> s^OSN!S>SSS! £§S59^N J^f^ S S S .< X ^ : 555^< ^S ^; £5S>5>Ss!V!J9« S§S!S?2S5»; nmmnmwwrsj’ ; IIIIUIH»»!i iiiiiiiiir h uvmnmvMVA \ummwTMVjr <a\ umuvt/.rjv.rsj mmYMYSSSSA \mmv//sssj:»w;, H 19//M YjVaV1VA\ W M V YSSM W 'S s: wawsmw/wA rmwmxY/Mm. MMV’/xvy/sssm mwtm.w&mm. UWXWMZXZtW. wmsssxssfSKasttt vrmv.v&xv.'ys. ;!7?5sesss«3??s»! 'xmwtsvw&sxtx ’m w m & zxsm 'mai&z&xzjryt■&tK4*C-X!f?}6;yZ\ Figure 3.8 Solar ray model superimposed on a fluid dynamic grid. 240 rays are shown with a 51 x 51 grid. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 92 Along each ray, the intensity is determined by Beer’s law, = -M m , (3.16) where m stands for the ruth ray. Note that Eqn. (3.16) is the same as Eqn. (3.11) without medium emission. For the non-gray problem, the spectral dependence of the absorption coefficient is represented by several bands within each of which the absorption coefficient is constant. Further, the incident intensity is also assumed constant across each wavelength band. Thus, a series of equations like (3.16) need to be solved (one for each band) for each m th ray. Equation (3.16) has an exact solution which is given in exponential form, fm (i) = /o,m exp[- / kda], J0 However, exponentials are expensive to compute and a numerical integration procedure is employed. For each ray, grid points are placed midway (at point M ) between the points of intersection with the constant £ faces of the fluid dynamic cell (see Fig. 3.9). The local intensity of a particular ray is taken as constant over each fluid dynamic cell. This ‘average’ intensity for each cell is taken as the intensity calculated at the mid-point M . This calculation is done by marching along the ray starting at the inlet. For the ( t,j) th cell, j=1 j=l which gives J. _ (3 1T) Equation (3.17) is employed successively starting from I\ which is specified at the inlet of the flowfield. The heat absorbed from one ray by any particular fluid Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 93 Fluid d y n a m ic g rid cell Figure 3.9 Intersection of a single solar ray with the fluid dynamic grid cells. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 94 dynamic cell (i , j ) is determined from the local intensity the path length A 3{j and the absorption coefficient k{j. The total heat absorbed by each cell is obtained by summing over all rays that pass through that cell. The absorption term in Eqn. (3.10) is then given by, NR = V L. .A Tm (Vol. of cell) (3.18) v ' where N R is the total number of rays. This procedure guarantees conservation of energy and is easy to use. 3.4 Solution Techniques for Radiative Transfer 3.4.1 Exact Integral M ethod The ray-tracing technique accounts only for absorption of the solar beam. Re radiation from the hot gas, reabsorption by the gas, radiative heating of the wall and losses to the surroundings are represented by the V*q^ term in the energy equation and this gradient term has to be modelled. Equation (3.15) represents the governing radiative transport equation in cylindrical co-ordinates. Recently, Thynell [88,89] has shown that Eqn. (3.15) may be transformed into an equivalent integral form of the Fredholm type. The radiation problem then reduces to one in which the dependent variable is the radiosity of the peripheral surface The resulting integral equation for the radiosity involves surface-to-surface and gas-tosurface kernel evaluations. The radiosity is then expressed as a standard power series and the unknown coefficients of this series are determined using a colloca tion strategy based on the zeroes of the Chebyshev polynomials. The surface and volume integrals are evaluated by using Gaussian quadrature. Once the radiosities are known, the net radiant heat fluxes and incident radiation may be computed. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 95 Equation (3.13) may then be applied to determine the gradient of the radiative flux vector which is the required source term in the energy equation. The mathematical details of the formulation are given in Refs. 88 and 89. In principle, wavelength dependence may be accounted for accurately by in corporating a Gaussian quadrature scheme. This would involve multiple solutions of the integral equation for different wavelengths. The procedure would then be come computationally very expensive. Secondly, the present formulation does not account for spatial variation of the absorption coefficient. Again, in principle, this should be straightforward to include but the evaluation of the integrals would be come tedious and very time-consuming. Thus, this integral solution procedure may be used only to solve gray gas problems with a constant absorption coefficient. Even for this problem, if highly accurate integration involving a large number of quadra ture points is desired, the scheme is expensive. Due to these limitations, the exact solution procedure was used only as a means of evaluating other more approximate schemes over different optical thicknesses. The results of this study forms the last section of this chapter. 3.4.2 The P -l Differential Equation M odel Another general method to solve Eqn. (3.15) is to expand I \ in an infinite series, in terms of spherical harmonic functions. For the simplest approximation, which is the P - l, the infinite series is truncated after the first order terms [86,87]. By taking the zeroth and first order moments of equation (3.15), the P - l differen tial approximation for a gray absorbing-emitting (non-scattering) medium may be derived as: d 1 dG 1 d y dG , . S i h H + y f y k f y = 3 k { G - i *h ) ' . . (3'19) V • qfl = k(4irlb - G). (3.20) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 96 For a given temperature field, Eqn. (3.19) may be solved numerically for the incident radiation G(x,y). This may be substituted in equation (3.19) to obtain the gradient of the radiant flux in the energy equation (3.10). Thus, the energy equation for the temperature and the P - l equation for the incident radiation may be solved iteratively to obtain a converged temperature-radiation field. The solution of Eqn. (3.19) requires the specification of four boundary con ditions - i.e., at the inlet and outlet for the axial direction and at the centerline and wall for the radial direction. Marshak’s projection scheme [90] is normally used to specify the boundary conditions. The boundary condition at the wall (y = R), which is assumed to be opaque, gray and diffusely emitting and reflecting, may be derived as: + |( 1 + U ) 6G(-g’R) - 4irewh ( T w), (3.21) where eu, is the wall emissivity, ( w is the wall reflectivity and Tw is the wall tem perature. Similar boundary conditions may be derived for the inlet and outlet and are given in Ref. 90. The inlet is generally assumed to be purely transmitting and the outlet to be purely reflecting. The latter condition is employed to account for the radiation that is incident on the outlet surface from downstream of the compu tational domain (i.e., from the nozzle portion). For the centerline, the symmetry condition is employed. The governing energy and P - l radiation equations may then be written in the coupled vector form: 8Q dE J t + dF 8 E V 8FV aT + „ ( 3 ' 22) where Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 97 KeVdT\ ( tU3LQL\ V c ). fc 07 / V / i±sJL2A-\ Key o r F « = ( Cis a \ 3fcy(G - 47rJ6) ) . k 6y / / The first element in each of the vectors is the energy equation while the second element is from the P - l equation. The numerical procedure employed to solve Eqn. (3.22) is similar to the time-marching technique employed for the fluid dynamic equations. In fact, the unsteady terms of Eqn. (3.22) have been included only for the purpose of marching the solutions in time until a steady state is reached. The details of the computational algorithm are given in the next chapter. When the absorptive properties are highly frequency dependent, the radiation transport equation must be solved on a detailed spectral basis. We may then con sider the spectral versions of Eqns. (3.19) and (3.20) and integrate these over each wavelength band. The details of this procedure are outlined in Ref. 90. The re sulting P -l model then involves the simultaneous solution of a series of differential equations - one for each wavelength band considered. When the number of bands is larger than two or three, the P -l model becomes very time-consuming and is not economical for non-gray problems. In the present work, the P -l model was employed only to study absorption of solar radiation by small particles. Assuming the medium to be gray, it was possible to conduct extensive parameter surveys on the absorption and reradiation characteristics of small particles. The results of this study are included in Chapter 5. 3.4.3 The Thick-Thin Model We have seen that the radiative properties of alkali vapors are highly frequency dependent. Consequently, it is necessary to employ a model that describes this Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 98 spectral variation and is, at the same time, economical. A closer look at the radiative properties of the working fluid suggests that a simpler model may be applicable. The absorption coefficient is, in general, in two parts — photoionization and dimer bands. The molecular (or dimer) bands are very strong at low temperatures and these bands cover a substantial fraction of the wavelengths of interest. It is expected that the gas is optically thick over these wavelength bands for most of the flow-field. Thus, a simple diffusion approximation, which is accurate in the optically thick limit, may be sufficient to model this contribution. In the ultra-violet region and part of the visible wavelengths, photoionization bands occur. The absorption coefficients in this region are quite small except at elevated temperatures. So, these bands are likely to be optically thin over most of the flow-field and their contribution may be described by an optically thin (or purely emitting) model. The thick-thin model [34,35] ensures that the predictions are accurate when the optically thick and thin limits are approached. The band model described in Eqn. (3.9) is used to describe the spectral variation of the absorption coefficient. In a region where a radiation band is optically thin, the optically thin approximation is used to determine the V-q^ term. In a region where a radiation band is optically thick, the radiation flux is described by the Rosseland diffusionapproximation. Criteria are established to determine whether a particular band is optically thick or thin. In order to prevent abrupt transition between the two models, a gradual transition is enabled by means of an interpolation function. For optically thin bands, emission dominates over absorption. This gives V • q f = ATrkxI bx. (3.23) Integrating this over the particular wavelength band,we get (V -2,')khi l = 4 r t pA . (3-24) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 99 where kp{ is the Planck mean absorption coefficient which is given as _ <rT4 — ’ where F; is the black-body fractional for the ith band. For optically thick bands, the diffusion approximation [75] is used, 2? = ~ T Z W u - (3 '25 ) Integrating this over each wavelength band, we get (v -S ?)“ ? = - y v -(ffrV T ), ( 3 .2 6 ) where K r is the Rosseland diffusion coefficient and is given by K' - L h w » - h L w » At any point in the flow-field, a given band is taken to be either thick or thin, so the radiation term is given by either Eqn. (3.24) or (3.26). The criterion for transition from thick to thin is provided by the product of the absorption coefficient and the radius R of the chamber. To avoid an abrupt transition between the thick and thin limits, aninterpolation function is employed. The form of this interpolation is obviously im portant to the overall accuracy of the model.The P -l model was employed to calibrate the results of the thick-thin model for intermediate opacities. On the basis of this study, an appropriate form of the interpolation function was determined. k{R < 0.2 thin kiR > 5.0 thick (3.27) 0.2 < k, R < 5.0 interpolate by Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 100 /(*«*) 1 + (5 - k i R Y Then, (V • qRU A, = II " /(A>. JZ)](V • , B) S " + /(* ,fi)(V ■2B)SSk- (3.28) The complete gradient source term in the energy equation is evaluated simply by summing up Eqn. (3.28) for all wavelength bands. The thick-thin method does % not involve the solution of any additional differential equations for the radiation. Then, any number of wavelength bands may be considered with very little overhead to provide an accurate spectral representation of the absorption co-efficient. This is the main advantage in using the thick-thin model. 3.5 D em o n stra tiv e S olutions In this section, we consider the results of the radiation problem in the solar thruster using the ray-tracing technique for solar absorption and the exact, P -l and thick-thin models for reradiative transfer. The accuracy of the two approximate methods are assessed using the exact solutions over a range of optical dimensions. These calculations also provide a preliminary picture of the thermal field within the thruster. The geometry employed in these calculations is shown in Fig. 3.10. The calculations are restricted to the absorption region and the downstream nozzle region is not included. The boundary conditions on the radiation problem are also shown on Fig. 3.10. The upstream boundary which is the len6 window is assumed to be completely transparent to thermal radiation. The peripheral walls are chosen to be diffusely gray and the downstream surface to be perfectly black. For these sample calculations, the concentration ratio of the focusing lens was taken as 104 : 1 and the power in the inlet solar beam was 1 MW. The corresponding Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 101 T ra n s p a re n t S urface O p aq u e W all Black W all Figure 3.10 Geometry used in the computations with the boundary conditions. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 102 focal spot size of the incident beam is 0.2 m diameter. To accomodate the focal spot and to allow for beam convergence within the absorption chamber, the window diameter was taken as 1 m. The inlet temperature of the gas was 300 K. A grayconstant absorption coefficient was chosen to simplify the calculations. Figures 3.11-3.13 show the contours of the gas temperature for absorption coefficients of 0.4, 2.0 and 6.0 m ~ x respectively. Results from all three models are included in each figure. The outline of the incident solar beam is also indicated on the contours. In Fig. 3.11, k = 0.4 m -1 or the gas is optically thin. Much of the solar radiation is transmitted through the gas and only a small fraction of the energy is absorbed. Most of this absorption occurs in the region of the focus because the intensity of the radiation is largest here. The peak temperature of the gas is only 2450K. The three models show a remarkable degree of agreement. This is not too surprising since at these low opacities, the gas reradiates very weakly. The tem perature field is then basically absorption-driven and since the absorption is modelled identically in the three schemes, we would expect the solutions to be close and they are. Figure 3.12 shows the results for an intermediate opacity (fc=2.0 m -1 ). The absorption from the solar beam is nearly 100% and because the absorption path length (inverse of the absorption coefficient) is 0.5 m, most of the absorption again occurs in the region of the focus. The peak gas temperature is now about 3600K and the reradiation from the plasma is significant. Still, the tem perature contours are predicted quite well by the P -l and thick-thin models showing that the radiative fluxes are fairly well represented by these models. The maximum discrepancy in the temperatures is only about 5%. Note that the thick-thin model results for intermediate opacities have been calibrated to match the P -l solution. Figure 3.13 shows the temperature contours for a higher optical thickness of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 103 EXACT 700 K 105Q K P-1 THICK-THIN Figure 3.11 Comparison of exact, P -l and thick-thin solutions. k = 0.4 m " 1, P = 1 MW, C R = 104 :1. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 104 EXACT 700 K 1050 K 3500 K P-1 3500 K THICK-THIN 3500 K Figure 3.12 Comparison of exact, P -l and thick-thin solutions. k = 2 m ~ \ P = l MW, C R = 104 :1. Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission. 105 EXACT 2100 K P-1 2100 K THICK-THIN 2100 K Figure 3.13 Comparison of exact, P-1 and thick-thin solutions. k = 6 m - \ P = 1 MW, C R = 104 :1. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 106 Xs=6.0 m - 1 . Now, the absorption path lengths are shorter and a large percentage of the radiation is absorbed near the inlet window itself. Consequently, the gas heats up upstream of the focal region. The peak temperature of the gas is only about 2200 K. This is because the solar beam in the focal region is quite weak (as a result of the aborption) and also because the radiation losses from the plasma is very high resulting in a more uniform temperature distribution in the gas. The results also show that the P-1 and thick-thin models underpredict the peak gas tem perature by about 10%. The maximum discrepancy occurs around the focal spot; this is a consequence of the importance of reradiative transfer for this gas opacity. Still, the overall shape of the temperature contours are in good qualitative agreement. These results demonstrate that both the P-1 and thick-thin models give very good predictions of the gas temperatures and the radiative fluxes in solarsustained plasmas. Figure 3.14 shows the distribution of the radiative heat flux along the wall for two optical thicknesses (fc = 1 and 3 m - 1 ). Results are shown from the P-1 model and the exact method. Overall trends are fairly well predicted by the P -1 model. For the higher optical thickness, the P-1 model underpredicts the wall heat flux by about 20%. This may be a consequence of the lower gas temperatures evident in Fig. 3.13. Furthermore, the P-1 model is notoriously inaccurate near boundaries [86]. The wall heat fluxes from the thick-thin model are not shown because of the spatial distribution of the heat flux from thin radiation was not determined. This involves calculation of volume-to-surface integrals which are quite complicated; because there is little to be gained from knowing the distribution accurately, the calculations were not performed. On the whole, the above results are quite encouraging. The P-1 model and the thick-thin model succeed in representing the temperature fields well. The wall Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 107 0.004 Exact 0.003 k=2 m * 0.002 0.001 k=6 m 0 .0 0 0 0. 0 0.75 0.375 AXIAL LENGTH, m Figure 3.14 Distribution of radiative flux along the peripheral wall. Non-dimensional heat flux, g* = Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 108 heat fluxes are also reasonably well-predicted. The approximate results show the same overall trends with varying optical thickness of medium as the exact solution does. The approximate models have been applied to evaluating solar thruster con figurations over a range of flow conditions and power sizes. The results of these parametric studies are presented in Chapter 5. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 109 CH APTER 4 THE FLUID DYNAM IC MODEL 4.1 Solution o f the Navier-Stokes Equations 4.1.1 The Navier-Stokes Equations The fluid dynamics in the thruster chamber is governed by the Navier-Stokes equations. These equations may be written in axisymmetric co-ordinates in the following vector form: dQ 8E dt "** dx dF dy ^ +^ + B ,+ B „ Ox oy (4.1) where ( puv \ (pu2 + p)y E = puvy \ (e + p)uy } Q= Ev = F = 0 w ft + w if Fv = 0 \ -f£ (H P - 3^7 + _ §u f£ 0 \ ■ fw fc + i m t t \ ( pvy \ puvy (pv2 + p)y (e -f- p)vy J K*y f Hz / 0 _ . \ -P9V -fiJvSe + p J xH$ ~q"'y — (V • qR) y / Here, p, u, v, p and T represent the density, x and y velocity components, the pressure and temperature respectively, e is the total energy per unit volume Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 110 and is given by, e = pCvT + 1/2p{u2 + v2). The source terms are written as two vectors Hi and Hz. The vector Hi contains the standard fluid dynamic source terms because of the axisymmetric co-ordinates. The vector Hz contains the source terms arising from the interactions with electromagnetic radiation. The energy equation contains the q'" term which is the heat source from electromagnetic radiation. For solar radiation, it is given by Eqn. (3.18) and for microwaves, it is given by J • E. Further, the V ■q R term standing for radiation energy transport is also included in the energy equation. In the momentum equations, the source terms arising from buoyancy and the Lorentz forces are included. These sources are used only in the microwave analysis and not in the solar heating cases. The viscous dissipation terms in the energy equation have been dropped because because our primary interest is the fluid flow in the absorption chamber where the Mach number is low and these terms are negligible. However, we are also interested in nozzle flowfield solutions in order to predict thrust and specific impulses of these devices. For such calculations, these terms must be included as well. For the purposes of the numerical algorithm, it is convenient to rewrite the viscous part of Eqn. (4.1) in terms of a new dependent variable Qv. Then, we get: dx dy dx xx dx dx xy dy dy yx dx dy SO . yy dy where Qv = { p , u , v , T ) T and, 0 ivy 0 0 0 0 \ 0 0 1 My 0 ’ 0 K cy J II St H 0$ /° 0 Rxx — 0 \0 0 /° 0 0 0 m \0 0 0 -f^ y 0 0 0 0 5 0/ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I ll /O 0 0 fly Ryy — 0 0 \0 0 0 Ryx = I o -§M y o ,0 0 0 o 0, 0 0 I ny 0 0 \ 0 0 K cy l In order to accomodate complex flow-field geometries, it is convenient to ex press Eqn. (4.1) in a generalized co-ordinate system. In other words, the physical co-ordinates (sc, y) are mapped onto a rectangular equally spaced computational do main (Z, 77). This transformation is achieved by appropriate functions, £ = £(®,y) and 77 = 77(3 , 2/). The transformed set of equations is given by [67,68]: 8Q dE 8F where Q - Q/J, J — txtfy E = ( U E + £yF ) / J , r d D ~ tfxZy F = (77mE + 77yF ) / J d d D d e _ d a + dZ ^ d r i + dr) dZ + drj d dr}' The tildes on the flux vectors should not be confused with the Maxwell equations (Eqn. 2.10). Henceforth, the tildes on the transformed variables are dropped. The transformed viscous matrices are given by: R& = E-xx d- ZxZy{Rxy d* Ryx) d- Zy Ryy]l Rfr — ^ x'HxR xx 4" ZxVyRxy d~ ZyVxRyx d* ^yVyRyyl/ Rif$ = [ZxfJxRxx d- ZyVxRxy d- ZxVyRyx d* ZyVyRyy\/ RflV = [Vx Rxx d- VxVyiRxy d* Ryx) d* tfy Ryy}/ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 112 4.1.2 Time-Iterative Solution Procedure The Navier-Stokes equations represent a mixed hyperbolic- elliptic set of equa tions. The convective or inviscid flux terms form the standard Euler equations which are hyperbolic in time. The eigenvalues (u , u , u + c and u —c, where c is the speed of sound) of the Jacobians of the flux vectors E and F represent the charac teristic speeds at which disturbances are propagated. Early CFD algorithms for the Euler equations [91,92] used these physical characteristics to drive out disturbances and relax the flow field to a steady-state solution. Since then these time-marching algorithms have progressed rapidly and have become the preferred method for com puting a wide variety of flows. For instance, the same general philosophy has been extended to the Navier-Stokes equations [93,94] even though the meaning of the characteristics speeds loses its physical significance (because of the ellipticity in troduced by the viscous terms). In addition, preconditioning schemes where the physical speeds of the system are altered in order to accelerate convergence have also been developed [95,96]. In this section, the algorithm for the solution of the standard Navier-Stokes equations is presented. Special procedures for low Mach number calculations are discussed in the following section. The Euler Implicit scheme [67] is applied to march Eqn. (4.3) in time to a steady state. The discretization is first-order backward differencing in time and second-order central differencing in space. The left hand side operator of the result ing delta form is approximately factored following the Douglas and Gunn procedure [97]. This gives rise to two block tridiagonal operators that may be solved sequen tially using the Thomas algorithm [67]. The approximately factored algorithm may then be written as [67]: dA d_ 8i dV T - i dJB _ d_ 9V dii drj nn dr) AQ = -nn Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. (4.4) and the Jacobians A, B, V and D are defined as: A-?E BQ’ B = Ja i BQ’ B - — BQ’ BQ' Note that the viscous cross-derivatives have been treated explicitly in order to pre serve the block tridiagonal structure of each of the factored operators. In general, fourth-order artificial dissipation is added in order to prevent odd-even splitting [67], In Eqn. (4.4), only the source vector Hi containing the fluid dynamic source terms is treated implicitly. The source vector Hi containing the electromagnetic interaction terms are treated explicitly. Consequently, the fluid dynamic solution is said to be explicitly coupled to the electromagnetic field. In other words, the EM field is allowed to lag behind the fluid field. Implicit coupling would mean that both the EM field and the flowfield are advanced in time simultaneously. Generally, implicit coupling or the implicit treatment of these im portant source terms is pre ferred in order to enhance the convergence characteristics of the flowfield solution. However, the solution procedure of the solar and microwave fields render the implic it treatm ent of these terms difficult and often impossible. The specific difficulties with the solar and microwave absorption terms are however quite different and are discussed below. For direct solar absorption, where the ray tracing technique described in Chap ter 3 is used, the heat source at each grid location depends on the solution of the radiative transfer equation along each ray. This means that the heat source at a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 114 point is coupled to the solution at several points upstream of that point. Implicit coupling of the radiation source to the fluid dynamic solution would give rise to a full m atrix operator on the left hand side of Eqn. (4.4) rather than the desired five banded structure. Solution of this new matrix operator would require the same level of computational effort as the direct inversion of the Navier-Stokes operator— which is not practical and is what we seeked to avoid in the first place. For the microwave absorption problem, the situation is peculiar. The Maxwell equations are solved in physical time and the resulting time-harmonic field is averaged over a wave period to determine the average heating source and Lorentz forces. On the other hand, the fluid equations are marched in time iteratively to a steady-state. Thus, the two solution procedures, though similar in principle, cannot be coupled implicitly. The one exception is the radiative transfer term in the energy equation. When approximate methods such as the P-1 or the thick-thin model are used, this source term may be treated implicitly. However, when integral methods are used to determine this term, implicit coupling is again no longer possible. The proper specification of boundary conditions is extremely im portant. For inviscid flows, the Method of Characteristics (MOC) procedure suggested by Rai and Chaussee [98] may be implemented. In regions where viscous effects dominate, this procedure is not correct and appropriate viscous boundary conditions [67] must be used. For the present study, the flow inlet boundary was treated with the inviscid boundary procedure (i.e., MOC) while the wall and outflow boundaries, where viscous effects are likely to be im portant, were treated with viscous procedures. For subsonic flows, the acoustic waves ensure that three eigenvalues are al ways positive while one is always negative. This means that we must specify three boundary conditions at the upstream end and one at the downstream end. The remaining conditions that determine AQ at these boundaries (four conditions are Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 115 needed at each boundary) must come from the equations of motion. At the up stream end where MOC is used, this is done by premultiplying Eqn. (4.4) by the modal m atrix containing the left eigenvectors of the Jacobian A 0 'as decribed in Chapter 2) and then multiplying this result by a selection m atrix that selects those characteristic equations that represent information flowing toward the bound ary. At the downstream end, three flow variables—the two velocity components and temperature—are extapolated from inside the flow field. The boundary conditions specified at the upstream end are stagnation temperature, flow velocity and entropy. At the downstream end, the single boundary condition specified is constant static pressure. This pressure fixes the pressure level in the system. At the wall surface, no-slip viscous boundary conditions are imposed. These are u = v = 0 and T = T,vau- To derive a condition on the pressure, a convenient one is to use py = 0 which follows from the classic boundary layer assumption. Finally, for the centerline, symmetry conditions are applied in lieu of a boundary condition. 4.2 Low M ach N u m b e r F o rm u latio n Conventional time-marching techniques are found to be inefficient at low Mach numbers. This is because of the disparity between the acoustic speeds and the particle speeds which leads to a stiffess of the equation set. Typically, this stiff ness results in the rapid deterioration of convergence at Mach numbers less than about 0.1. Preconditioning methods [95,96] have been devised for low Mach num ber computations whereby the temporal derivatives axe modified so as to yield a well-conditioned set of eigenvalues. In physical terms, these methods remove the rapidly propagating acoustic waves from the equation system and replace them by a new set of ‘pseudo-acoustic’ waves that travel at a speed that is of the same order Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 116 as the particle velocity. These artificial acoustic waves then become the agents that propagate errors out of the flow-field and drive the solution to a steady state. An alternate approach to this problem was suggested by Merkle and Choi [99,100] who used a perturbation expansion in Mach number squared to obtain a new set of equations that are valid at arbitrarily low Mach numbers. In principle, the method is closely related to the preconditioning method. Traditional time-iterative schemes were applied to solve the low Mach number system of equations for viscous flows and convergence rates that are independent of either Mach number or Reynolds number were obtained [97,98]. This latter approach has been demonstrated to be quite robust and accurate and for this reason, we have chosen it for our analysis. According to this procedure [100,101], all the flow variables are expanded as: P = Po + tPi + ••• u = uo + eui + ...etc. where e = 7 Moo2. Upon substituting these expressions in the Navier-Stokes equa tions, we may write out a sequence of equations for like powers of e. The equation to order e_1 says that the zero order pressure, P0 is a constant. The low Mach number viscous equations to order unity in axisymmetric co-ordinates may then be written as: dQ dt >— L j dE dx dF 8y 1------------ d „ dQv "T d dQv . d dQv d p dQ v r 7T" dx “xy dy T dy**yx dx + a i R n ~ W + H (4.5) where 1 0 0 P 0 V 0 P T 0 0 u fpiy/P} °\ 0 0 P' , Q= \ uy vy Ty J Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 117 ( puy (pu2 + epi)y E = puvy \ puTy / - pvy puvy (pv2 + ep1)y pvTy 0 \ - |& ( H CPl - 3 ^ + f p f e - f Vf j V 0 / The viscous terms and the source terms (Hz) introduced by the interaction with the solar and microwave radiation are the same as in Eqn. (4.1). The eigenvalues of this system of equations are A[P-1^4] = u } u } (u + c ')/2, (u —c')/2 where A is the Jacobian (8E/dQ)und the ‘pseudo-acoustic’ speed c' is given by c' = \ / u 2 + 4e/3, e= Here, in order to get well-conditioned eigenvalues /3 is scaled as R T r. 4.3 Dem onstrative Flow Field Solutions 4.3.1 Sphere-Cylinder Configuration The flow field th at we investigate for the microwave problem duplicates the sphere-cylinder arrangement being used in the Penn State experiments. Flow enters a 2.5 cm cylindrical tube from the top (see Fig. 4.1) and expands into a 10 cm diameter sphere. The flow leaves the sphere through a similar cylindrical tube attached to the bottom of the sphere. The microwave cavity surrounding the flow tube is also shown with the no-loss field pattern for the TMoiz mode. The microwave discharge is typically maintained in the middle of the sphere. Entry gas velocities Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 118 MICROWAVE CAVITY SPHERECYLINDRICAL TUBE INLET WAVEGUIDE Figure 4.1 Fluid dynamic grid for sphere-cylindrical tube with the microwave cavity. Right side of the cavity shows the no-loss field pattern of the TM 012 mode. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 119 range from 0.02 to 2 m /s and discharge pressures from 0.4 to 4 atm. For helium at room temperature these conditions correspond to Reynolds numbers ranging from 0.8 to 800 and to Mach numbers from 10-5 to 10- 3 . To get an initial feel for the fluid dynamics in this flow configuration, some cold flow and specified heat addition cases were computed. These cases were computed with a moderate 81 x 41 grid size. A typical convergence plot of the low Mach number computations of the flow in the sphere-cylinder geometry is shown in Fig. 4.2. The Reynolds number of the flow was 80 and there was no heat addition. The figure shows the change in the solution vector (i.e., vector Q) as a function of iteration steps. The convergence is seen to be approximately linear on the semi-logarithmic plot. There is about 6 orders of magnitude reduction in error over 500 iteration steps. The Mach number of the flow is of the order of 10-4 . The low Mach number algorithm is thus seen to be very effective. Similar convergence rates are obtained if the Mach number and/or Reynolds numbers are varied. At significantly higher Reynolds numbers (above 500), convergence deteriorates but this is a consequence of the grid geometry rather than the numerical algorithm. Figure 4.3 shows cold flow results for Reynolds numbers of 8, 80 and 400. At the lower speed, the streamlines diverge upon entering the sphere and the flow slows down. As the velocity is raised, the sphere becomes the site of a rather large toroidal recirculation cell. This flow separation becomes even stronger at Re=400. Also note that the recirculation region is forced toward the right as the flow velocity is increased. This severity of the recirculation region at higher Reynolds numbers is responsible for the slowing down of convergence referred to earlier. To see the effect of heat addition on the flowfield, a heat source was specified as a function of x and y in the energy equation. Figures 4.4 to 4.6 show these heated Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 120 1 -2 ,-4 b/ bv -6 -8 -10 100 200 300 400 500 ITERATIONS Figure 4.2 Convergence of the low Mach number algorithm for cold flow in the sphere-cylider geometry. 81 x 41 grid. Re=80. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 121 Figure 4.3 Streamlines in plasma tube for cold flow at various Reynolds numbers. Flow is from left to right. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 122 HEAT-ADD1TI ON PROFILE 1000 K TEMPERATURE 500 K 7000 K STREAMLINES Figure 4.4 Effect of specified heat addition on flow through spherecylinder configuration. Re=8, P=150 W. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 123 HEAT-ADDITION PROFILE 800 K TEMPERATURE 400 K STREAMLINES Figure 4.5 Effect of specified heat addition on flow through spherecylinder configuration. Re=80, P=300 W. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 124 HEAT-ADD 1TI ON PROFILE 800 K 1200 K TEMPERATURE STREAMLINES Figure 4.6 Effect of specified heat addition on flow through spherecylinder configuration. Re=80, P=300 W. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 125 gas results. In all cases, the top figure shows the distribution of heat addition; the middle plot shows the temperature contours; and the bottom plot shows the streamlines. In Fig. 4.4 where the Re=8 and the total heat flux is 150 W, the gas temperature is seen to reach 7000K. The streamlines show the flow decelerating in the region of heat addition and the flow is pushed toward the sphere surface. Figure 4.5 shows similar results for a Reynolds number of 80. Here, the peak gas temperature is only 4400K even though the heat addition is increased to 300W. This is due to the higher gas velocities. The velocity of the gas also pushes the isotherms downstream resulting in a asymmetric flowfield. Consequently, the recirculating bubble is moved rearward (compared with Fig. 4.3). The results of heat addition in a region upstream of the sphere center is shown in Fig. 4.6. This results in a dramatic change in the flowfield. The recirculating region is almost completely eliminated and the solution resembles the calculation at the lower Reynolds number. 4.3.2 Bluff Body Configuration Recent experiments with resonant cavity plasmas have used a bluff body in the flowfield in order to stabilize the plasma at higher powers. Accordingly, the computational model was also adapted to study bluff body stabilized plasmas. In this section, we look at some cold flow solutions for this new geometry. The grid employed in these calculations is shown in Fig. 4.7. The length of the flow tube is 0.2 m and its diameter is 0.1 m. The diameter of the bluff body is 0.04 m. The grid size used was 161 x 81. The inlet gas velocity and gas pressure were varied and the calculations were performed for a series of flow Reynolds numbers. The Mach numbers of the flow were once again of the order of 10-3 to 10-5 . Therefore, the low Mach number algorithm was employed in the calculations. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 126 llllllllllllllllllllllllllllllltlllllllll llllllllllI lllllllllllI f llllllB I I I I U I I ll I tl ll ll l l l l t l l l l l l l l l l l l l l l l l U l l l l I I I I I llllllllllllllllllllllllltlllttllllltlll IIIIIIIIIII1 IIIIIIIIIIIIIIIIIIIIIIIIIIH iiniiiiiiaiiiiiiiiiiiiiuiiiiiiiiiiuii iiiniiiiiiiitiiiiiiiiiiiiiiiiiiiiiiiiii lllllllf1 II1 IIIII1 I1 IIIllllllllltlllllflS III11III1IIIIIIIIIIIIIIIIIIII1IIIIIIIIII m m iim m tiiiiu im tiiiim m iiii m u tin m n u m im tittiiiiiiiiiiiii iH u iu u n tu tin m u u m itiiiiiiiti \u m n u \u m \im im iim u iiiiu n m u u n iiiiu u iitim m im iim iiit u u m m m titiiv m iiitiiiiiiim iin U M im m u u m m iu im m im iim u m u iiiiiiim iu im im iiiiiitiiiii u m m m m m m itiu u iiiiim m ii m u u u u iim u im iiim iiiiiiiiiiii u m m im u m m iu tm iiu im im i u m m m m m m tm im m im iin i ittmmiuiM iimimmiiumiiim tm m u um m m m im iim iiiin ii m u u m m um m iim m iim m m ummntmm m tiiu im itiim m i MICROWAVE CAVITY t in t iiiiiiiM iiiiiiiit iiin in iiiiiiiii IHlUIIIIIIH IIUM tltlH IM H IU IIIIII U im ilU lllllllltlllllllllllllM ltllll U lU U IIU IIIItllU IIIU IIIIIItlH llllt BLUFF BODY . . ................................ . mi mmm ..... . i ii n i ii ii ii ii ii ii ii ii ii im i ii ii m i ii ii ii ii ii n m ii n i ii im im i ii ii ii ii M ii mu iiiiiiiiiiii>ijiiiiiiiiiimiiiiiiuiitiiiiiiiiiiiiiiiiiiiiiiuiiiiiiMii;ii(ii mliiiiiiiiiiiiriiiiiiiiiiiiiiiuiiiiiiiiiiiiiiilli'iiiim iinmiiiiiiiiiii'i i ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii ii in ii ii ii ii ii ii ii ii ii ii ii in m i iu i i i i i i h ........... . ......... h i n i m t u m i i n i i i i i i i i i i i H i i i n i 1:11:11111111 n 11 1.1 in m i in m in i n u n 11111 n t iim iiiiiiuiiim iiiiiiiiiiiimiiiiiiiiiiiiiiiiiiiiiiiiiHiMiiiiiiiiiiiiiiiiiini I I I l 1 111 I II 11 i l l I I i l l I .................................. ................................................... H I 111 I I II M i l iiiiiiimiitii 11 mi i’i ii m u m m in i m m 11111111 inn in mm 1 1 m inim um 1 mi ........ 1 1 1 1 mu............................ mm mm I mini mm 1 1 1 1 I I 1 1 1 1 1 1 1 I 1 wmrnw 11niiiiiiiuiniiiiii 1 11111111111111 1 m i ................ 1 1 n i n 1 1 11 11 11 11 i i i t i i i i i n i i l i n i iiiitilliilliiiin 111............ 1 1 1 i i ’i u i n .1111.1 ittiiii ii n i ii iu i I ................. i l l II i l l i l l K l l i l l l i l l l l l I I * I I I II 111 II I I I I l l l l t l J l l t l l l l l l l l 1 1 111 m i 111 n t 11 n u t i i i i i i i i i i i i m . 1 1 1 1 1 1 1 1 1 1 ii. 11 iiiniiiiiiiiiniiii iiuiiiMHimii 1 1 1 111 11 111 hi 111 iiiiiiiniiniiiiti I 1 1 1 1 1 11111 in in 11 iiiiiiiiiiiinniin iiiiiiiiin 1 1 1 1111 i n n 11111111111111111111111 mu 1 1 1 1 1 1 ri in 1111111111111111111111 mini 1 1 i n n 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 11 1 1 1 . 1 1 1 i n n 1 ill 11 ill ill Mil m i||(1 Mltl|ll1 |l1 ll)l1 l|MIMI||lllll|HliniMlinMI|l||llirilMMIIIIIM I I mi 1 i 1 1 i 1 1111 1 i n 11 111111 u i i i i i i i t i i i n i i i i i 1 1 I I I II I I I II I I I I U I I I II I II I I I III Figure 4.7 Fluid dynamic grid for axisymmetric bluff body geometry. Flow is from top to bottom. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 127 Figure 4.8 shows the convergence history for a series of cases. The Reynolds number of flow (based on inlet conditions and tube diameter) is given alongside each of the curves. The first case was computed with a uniform flow as initial condition. Subsequent cases were computed by using the flowfield from the previous case as the intial condition. Convergence rate of the calculations is quite slow for the lowest Reynolds number shown (Re=5). There is barely a two orders of magnitude reduction in the error over 1000 iteration steps. The reason for this poor convergence will be evident in Fig. 4.9 which shows the flowfield solutions for different cases. In Fig. 4.8, the convergence rate steadily improves as the Reynolds number increases. For the highest Reynolds number shown (Re=1000), the error in the solution decreases by almost six orders of magnitude over 1000 iteration steps. Figure 4.9 shows the streamline contours for four different flow Reynolds num bers (Re=5, 60, 350, 1000). At a Reynolds number of 5, the flow is very viscous and almost creeps along the surface of the bluff body. The recirculating bubble in the wake of the bluff body is very small and the flow makes a 90 degree turn. This drastic switch in flow direction as well as the predominant viscous effects make this a difficult case to compute—thus explaining the poor convergence in Fig. 4.8. Increasing the Reynolds number steadily increases the size and the strength of the reirculating region in the wake of the bluff body. At the highest Reynolds number shown (Re=1000), the length of the recirculating bubble is approximately 1.6Z?6 where Db is the diameter of the bluff body. Also, as seen from Fig. 4.8, the nu merical algorithm performs very efficiently for higher Reynolds numbers. This is because of smaller viscous effects at these Reynolds numbers. It also appears that the recirculation itself does not cause any serious convergence problems. These initial cold flow results for the sphere cylinder and bluff body configu rations demonstrate the performance of the low Mach number Navier Stokes code Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without perm ission. O’ 't o o •J -10 2000 4000 6000 300 0 ITERATIONS Figure 4.8 Convergence of low Mach number computations for the bluff body geometry at different Reynolds numbers. 10000 o n II v Pi © II V Pi Figure 4.9 CO Bluff body computations for different Reynolds numbers. Re = 5, 60, 350, 1000. 129 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 130 over a wide range of flow conditions and for different geometries. The solutions also demonstrate that the apparently simple geometries are rich in complexity. The microwave heat addition is likely to have a profound effect on the nature of the flowfield. Clearly, prediction of these flowfields can be very illuminating and may aid in experimental investigations as well. Investigation of the coupled microwave gasdynamic field is considered in Chapter 6. More detailed parameter studies in cluding variation of pressure, gas velocity and microwave power level are presented. Validation of results with available experimental data is also carried out. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 131 CH APTER 5 FEASIBILITY STUDY OF DIRECT SOLAR PROPULSION 5.1 Direct Solar Absorption 5.1.1 Sizing o f the Solar Thruster The complete solar propulsion system consists of a collector to collect and focus the solar energy and a rocket engine in which the working fluid is heated and expanded through a supersonic nozzle. The size of the collector is dictated by the power requirements of the engine which are, in turn, determined by the thrust size of the engine. For a given optical focusing system, the size of the focal spot, or the size of the sun’s image in the absorption chamber, may be determined. The focal spot size prescribes the minimum diameter of the absorption chamber. We have used a simple procedure based on ideal vacuum expansion in order to size the collector and the engine. These one-dimensional relations are given below. Momentum balance for the whole system gives: T h = mVe Energy balance for the heat addition chamber: P = mCpTo Energy balance for the jet nozzle: Ve = y/2CjFo- Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 132 Here, T h is the thrust, Ve is the exit velocity of the jet, m is the mass flow rate, P is the power of the system and T0 is the average stagnation temperature of the solar plasma. Assuming the plasma temperature, the exit velocity may be determined. Fixing a thrust size then gives the mass flow rate and, hence, the power requirements of the thruster. Figure 5.1 summarizes the power requirements for different thrust levels and gas temperatures. Once the power requirements are known, the size of the collector may be fixed from the magnitude of the sun’s irradiation at the earth, 1315 W / m 2. Fixing a concentration ratio (CR) for the converging lens, the focal spot size (or the size of the sun’s image in the thrust chamber) may also be determined. The collector size and focal spot size are shown on Fig. 5.2 for a CR=1 X 104:1. The focal spot size is an important parameter in d e te r m in in g the engine diam eter. Very simply, the engine diameter must be large enough to contain the focal spot of the concentrated beam so that all the energy in the beam may be adm it ted into the chamber. It might seem that the focal spot can be made arbitrarily small by concentrating the sun’s rays adequately but this is not the case. The maximum intensity of the concentrated solar radiation is limited to the value of intensity of radiation at the sun’s surface and, consequently, the maximum possible concentration ratio is about 4 x 104:1. Corresponding to this concentration ratio, there is a minimum focal spot size that can be achieved by the focusing appara tus. This limitation in the CR is also responsible for the thermodynamic limitation in the maximum temperature that may be attained by the gas in the chamber; as mentioned earlier, this maximum temperature is 5800 K which is the surface tem perature of the sun. In practice, CR’s are limited to between 1 x 104 : 1 and 2 x 104 : 1 because there are practical difficulties with ensuring high quality optical properties of the reflecting surface [14]. In Fig. 5.2, we have employed a CR=1 x 104:1. For a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 133 6000 = 5000ft = 3000ft PO W ER (K W ) 4 0 0 0 ------ = 2000 ft 2000 0 200 400 THRUST 600 800 (N) Figure 5.1 Power requirements of a solar thruster for different thrust levels and gas temperatures. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1000 134 u . o —— To = 5 0 0 0 i f To = 3000/f COLLECTOR DIAMETER (10* m) FOCAL DIAMETER (m) To = 2000/f 0 . 4 ------ 0.2 200 400 THRUST 600 800 (N) Figure 5.2 Collector and focal spot sizes as a function of thrust size. C R = 1 x 104 :1 . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1000 135 N (40 lbs), the power requirement is 1 MW (Fig. 5.1) and the corresponding focal spot diameter is 0.2 m (Fig. 5.2). This sets the minimum window diameter of the engine at 0.2 m. For larger thrust and power levels, the focal spot is correspondingly bigger and so is the solar engine. Table 5.1 lists the power requirements and engine sizes for different engine thrust levels. These estimates are for an average plasma temperature of 3000 K, which is a reasonable estimate as our results will show. Note that employing a higher CR reduces the engine size; for this reason, we have employed a CR = 2 x 104:1 in most of our calculations. Still, the solar propulsion system is relatively larger than chemical, laser or microwave systems. The large size is a direct result of the dilute nature of solar radiation and is one of the drawbacks of solar propulsion. 5.1.2 Issues in Direct Solar Absorption There are three major issues facing direct solar propulsion—the efficient cou pling of the solar radiation to the thermal energy of the gas, the performance po tential of the direct solar thruster, and structural considerations introduced by the relatively large sizes of these thrusters. The primary issue concerning the feasibility of direct solar propulsion is the efficient transfer of energy from the solar beam to the working fluid. The major difficulty with coupling the solar energy to the thermal modes of the gas arises from the relatively weak intensities of solar radiation even after focusing. In previous chapters, we have discussed two possible ways by which absorption of solar radia tion may be achieved; these are absorption by small submicron size particles and absorption by alkali vapors. In both cases, hydrogen is selected as the carrier gas by virtue of its low molecular weight. In this chapter, both these absorption schemes are studied for practical solar thruster configurations. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 136 Table 5.1 Sizing of the solar thruster for different thrust levels and concentration ratios. D c - Collector diameter D f \ - Focal spot diameter for C R — 1 x 104 :1 D f 2 - Focal spot diameter for C R = 2 x 104 :1 Th, N P, MW D c, m Df i , m D f2, m 100 0.5 20 0.2 200 1 500 1000 2.5 5 30 50 0.3 0.5 0.7 0.15 0.21 2000 10 70 100 1.0 0.35 0.50 0.70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 137 A further issue that is related to the transfer of energy to the gas is the radiation loss from the solar-sustained plasma. In other words, we are concerned not only with how much of the incident radiation may be absorbed directly by the gas but also with how much of this energy is lost due to reradiation from the hot gas to the cooler surroundings. As mentioned in Chapter 1, the reradiative flux to the wall as well as incident radiation that is transmitted through the gas may be used to regeneratively heat the inlet gas. The regenerative heating is an indirect means of heating the propellant; thus, it seems that a “hybrid” concept that embodies both direct and indirect heating is necessary for high overall efficiency of energy transfer. Accordingly, we define two efficiencies—the collection efficiency, rjc, and the overall efficiency, t/0. The collection efficiency is the percent of the incident power that is directly absorbed and retained by the gas after reradiation. The overall thermal efficiency is the net fraction of the total incident power that is coupled, both directly and indirectly, to the gas. Note that the term coupling (or absorption) efficiency is sometimes used to define the percent of the incident power th at is directly absorbed by the gas. This is not so useful because it does not include radiation losses from the plasma; however, it is commonly applied in the study of laser and microwave plasmas. The second major issue in solar propulsion is its performance capability. Per formance is usually measured in terms of the specific impulse (Isp) of the propulsion system. Indirect solar propulsion promises specific impulses in the range of 500-600s [15,16]. These estimates make indirect systems very attractive especially since these systems are relatively simple. On the other hand, direct systems present a signifi cant technological challenge. For this reason, direct systems should be capable of providing better performance than indirect systems. Because of the significantly higher gas temperatures that are possible with direct absorption, direct systems do Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 138 promise better performance. But the question—whether this promised high per formance is attainable in practice—still remains. In fact, further development of the direct concept depends, to a large extent, on how well it can compete with the indirect concept. Finally, the relatively large size of the solar engine presents major technological challenges. Large sizes introduce structural difficulties of the thruster chamber as well as the inlet window of the thruster. In practice, these structural considera tions are likely to limit the pressure in the absorption chamber to fairly moderate levels—in the range of 1 to 3 atm. The pressure limitation strongly influences the performance characteristics of the direct solar thruster as our results will demon strate. In fact, size, pressure and structural weight considerations are likely to be the major questions facing the direct solar propulsion designer. In the following sections, the performance potential of both particle-absorption and alkali-absorption solar propulsion systems are assessed. Since the solar en ergy absorption process is inherently two-dimensional, a two-dimensional analysis of the propulsion systems is perfomed to provide realistic performance estimates. Parametric studies of these engines are performed; in particular, the effects of vary ing the gas flow rate and chamber pressure are investigated. The calculations are carried out for two CR’s (1 x 104 : 1 and 2 x 104 : 1) and three engine sizes (1, 5 and 10 MW). Further, the results of our two-dimensional model are contrast ed with previous one-dimensional estimates to demonstrate the inadequacy of the one-dimensional models. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 139 5.2 Direct Solar Propulsion Using Particle Absorbers The performance of direct solar thruster using small submicron size particles as the absorbing species is assessed by using the P -l differential model for a gray medi um. The absorption co-efficient of particles varies inversely with wavelength and the gray gas assumption is appropriate for the analysis. Accordingly, the calcula tions are performed by using the simplified radiation model described in Chapter 3; the thruster geometry employed is an axisymmetric cylinder as shown in Fig. 3.10. The calculations are restricted to the absorption region of the thruster and the choked nozzle throat region is not included. As mentioned earlier, radial velocities are neglected and the axial velocity component is assumed to be uniform. Further, regenerative wall-heating of the incoming propellant is not carried out. Scatter ing effects are also neglected; this assumption is generally valid for submicron size particles [87]. In all the calculations presented in this section, the power in the incident solar beam was taken to be 1 MW (thrust of 200 N) with a CR of 104 :1 of the collection system. The corresponding focal spot size is 0.2 m and the focal point was located near the downstream exit plane of the absorption chamber. The window diameter was fixed at 1 m to allow for beam convergence within the absorption chamber. The absorption coefficient of particles may be varied simply by changing the particle number density. Thus, any desired optical thickness of gas may be obtained by appropriately adjusting the fraction of particles carried by the hydrogen gas. Results of temperature contours and wall heat flux for different gas opacities were presented in Chapter 3. On the basis of those results, it appears that practical thrusters would operate at intermediate gas opacities. In this section, we examine the effect of gas velocity and wall emissivity on gas temperature and thermal efficiency. We Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 140 also present thrust and specific impulse results for the thruster by assuming ideal expansion of the exit gas in a nozzle. Figure 5.3 shows the effect of gas velocity on the temperature of the gas. The outline of the incident solar beam is indicated on these contour plots. Results for two different gas velocities are shown. In both cases, the absorption coefficient is 2 m -1 . Absorption (or coupling) is a function only of the absorption coefficient and is nearly 100%. In Fig. 5.3(a), the inlet gas velocity was 0.1 m /s corresponding to a Boltzmann number [47] of 0.68. The maximum temperature corresponding to this case was 4025 K. Fig. 5.3(b) shows the results correponding to a higher gas velocity of 0.4 m /s or a Boltzmann number of 2.74. At the higher gas velocity, the gas temperatures are significantly lower. The peak temperature now is only 3000 K. This reduction in temperatures is because of the associated increase in the mass flow rate of the gas in the thruster. The effect of gas velocity on the wall heat flux is shown in Fig. 5.4. Results from the P -l model as well as from the exact analysis of Chapter 3 are included on the same plot for comparison. The results from the two models are in good qualitative agreement with the P -l model consistently underpredicting the heat flux. The heat flux for the lower gas velocity is seen to be substantially higher than th a t for the higher gas velocity. This is of course expected as the gas temperatures are higher for the lower gas velocity as seen in Fig. 5.3. Naturally, it follows that at higher gas velocities, the reradiation loss becomes smaller while convective heat transfer becomes more important. This result points to the fact th at better thermal efficiencies may be expected at higher gas velocities. Figure 5.4 also shows the effect of wall emissivity. A perfect black wall (e = 1) absorbs all the energy that is incident on it. As the wall emissivity is reduced below 1, the wall starts reflecting (assumed to be diffuse) a portion of the radiation Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 141 t = 0.2 0.3 (a) SOLAR BEAM 0.3 zZ 0.8 (b) Figure 5.3 Temperature contours at two different inlet gas velocities (a) 0.1 m /s, (b) 0.4 m /s. t = ^ . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0. 008 Exact Bo=0.685 £=1 / 0. 006 Q Bo=0.685 '£=0.5 0. 004 B o=1.37 £=1 0.002 B o * 1.37 s=0.5 0 .0 0 0 o.o 0.375 AXIAL LENGTH, 0.75 m Figure 5.4 Radiative flux on the wall for different gas velocities and wall emissivities. q* = . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 143 incident on it. This reflection causes the net heat flux to the walls to decrease in proportion to the decrease in wall emissivity. The overall performance of the thruster is considered in Figs. 5.5 and 5.6. The collection efficiency, peak gas temperature, operating thrust and specific impulse are estimated over a wide range of mass flow rates (or flow velocities) and absorption coefficients. The predictions of thrust and specific impulses are first-order approxi mations since the nozzle calculations were not performed. Average exit conditions are determined by assuming ideal isentropic nozzle expansion. In Fig. 5.5, the peak gas temperature and collection efficiency (which is the same as the overall efficiency because there is no regenerative wall-heating) are de picted as a function of the mass flow rate through the thruster. The results are shown for two absorption coefficients — k = 2 m -1 and 6 m -1 . As seen before, higher mass flow rates cause a decrease in the peak gas temperature. This is be cause the same power (1 MW) is being coupled to larger gas quantities resulting in lower temperature rises. The higher optical thickness medium also shows lower gas temperatures. This is because smaller absorption path lengths result in heating of the gas upstream of the focal region. Lower temperatures mean smaller radiation losses and hence higher thermal efficiencies. Consequently, thermal efficiencies in crease with increasing mass flow rate and with increasing optical thickness. In fact, thermal efficiencies up to 100% seem possible with optimum mass flow rates. Note that there is a trade-off between efficiency and mass flow rate. High efficiencies re sult in more effective utilization of the incident solar energy but at the cost of using higher mass flow rates. Lower mass flow rates result in higher gas temperatures and hence higher specific impulses or better thrust per unit mass flow rate. Figure 5.6 shows the dependence of thrust and specific impulse on mass flow rate. Increasing the mass flow rate increases the thrust almost proportionately but Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 144 i1 J- i i I1... i i i i II i i i i II it Ii 1.00 ^- 0. 75 -jm ax f 0. 50 ^ o.25 - 0.00 1 -1 2m — *k -= f-! -1 — — k = 6m - t Figure 5.5 Collection efficiency and peak gas temperature as a function Boltzmann number. <ma!C = | Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 145 1600 1200 600 -4 0 0 - Th (N) (sec) 800- -200 400 0 1 Bo 2 3 Figure 5.6 Thrust and specific impulse of the direct thruster using particles for absorption. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 146 decreases the specific impulse of the thruster. In general, the results show that a 1 MW solar thruster is capable of delivering thrusts of a few hundred Newtons at specific impulses of 400 to 600 s. 5.3 Direct Solar Propulsion with Alkali Vapors 5.3.1 Details o f the Analysis The radiative properties of alkali metals exhibit strong spectral as well as tem perature dependence. It is, therefore, necessary to employ a detailed non-gray analysis of the radiative transfer problem in the solar thruster. Accordingly, the radiative properties are represented in the wavelength domain by a band model and the radiative transfer term is evaluated by using the thick-thin model. The details of these models were presented in Chapter 3. The thick-thin model is calibrated with the P -l differential model in order to ensure accuracy for intermediate optical dimensions. The demonstrative results in Chapter 3 show that the thick-thin model provides fairly high accuracy in its predictions. A practical solar thruster configuration is employed as the base geometry in these calculations (see Fig. 1.1). The grid geometry used is shown in Fig. 5.7; in this figure, the incident solar rays are superimposed on the flud dynamic grid. The walls of the thruster converge with the incident beam thus reducing the overall size of the thruster. The nozzle throat is also included in the computational domain so that the mass flow may be set by the choked throat rather than by an ad-hoc hypothesis. In order to simplify the computations, however, the flowfield is split up into two regions—the absorption region and the nozzle region. In the absorption region, where the wall convergence is slight, the transverse component of the contravariant Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 147 Figure 5.7 Solar thruster geometry with ray model superimposed on fluid dynamic grid. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission 148 velocity is neglected. In the more rapidly converging portion of the nozzle, the full fluid dynamic equations are solved. The propellant working fluid was chosen as a mixture of 95% (by volume) hydrogen gas and 5% alkali vapors (composed of equal mole fractions of Na, K and Cs). The absorption spectrum of the three alkalis cover about 80% of the solar spectrum and thus provide reasonable absorption over most of the spectral range. The total pressure in the absorption chamber was taken as 3 atm because of possible structural limitations. Under these conditions, the molecular weight of the propellant is about 5. As a parametric study, the performance of the solar thruster was evaluated for a higher pressure of 10 atm as well. Reradiation loss from the hot gas is an important issue and the necessity of us- ing regenerative wall-heating to provide higher overall energy conversion efficiencies was discussed earlier. In all the calculations presented in this section, the regen erative heating of the inlet propellant gases is included. By regenerative heating, the propellant mixture is preheated to a temperature of about 1000 K; this ensures th at the alkalis in the mixture are in the vapor state. We consider three engine sizes corresponding to thrust levels of 100, 1000 and 2000 N. The corresponding power levels are 1, 5 and 10 MW with collector sizes of 25, 70 and 100 m. These collector sizes are realistically deployable from a practical viewpoint. In the calculations, the concentration ratio was selected as 2 x 104 : 1; this is considerably below the theoretical maximum of 4 x 104 : 1 and is likewise attainable with realistic designs [14]. The focal spot sizes corresponding to these designs are 0.22, 0.5 and 0.7 m respectively. To allow for beam convergence inside the engine, the window diameters are selected to be 0.4, 1 and 1.5 m respectively. The results for the three engine sizes are presented in the following sections. General flowfield characteristics are first discussed followed by a detailed consideration of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 149 the efficiencies of energy conversion and the performance numbers of the thrusters. 5.3.2 Flow field C h a ra c te ristic s Temperature fields and Mach number contours are presented in Figs. 5.8 5.11 to demonstrate the overall characteristics of the absorption region. Figure 5.8 shows results for the 5 MW engine with a throat radius of 11 mm (m = 0.032 kg/s), while Fig. 5.9 is for the same power level but a larger mass flow (m = 0.128 kg/s) corresponding to a throat radius of 22 mm. Figures 5.10 and 5.11 show results for the 1 MW size ( throat radius 9 mm, rh = 0.02 kg/s) and the 10 MW size (throat radius 33 mm, m = 0.288 kg/s), respectively. All cases are for a chamber pressure of 3 atm. The temperature contours for the lower mass flow 5 MW engine (Fig. 5.8) exhibit a steep temperature gradient near the entrance where the absorption from the dimer bands is strong. Further from the inlet, as the gas gets hotter, the dimer bands dissociate and absorption occurs primarily because of the photoionization bands. The peak gas temperature of 3295 K occurs on the axis near the focal region. The gas in the outer periphery is considerably cooler than this because of the lower solar intensity in this region. Downstream of the focal region, the hot central core gases lose heat by radiation and mixing until the average tem perature at the nozzle throat is 2273 K. These temperature contours demonstrate the strong twodimensional character of the absorption region. The corresponding Mach number plots demonstrate the low speeds for this mass flow rate and the acceleration through sonic speed at the throat. Increasing the mass flow at this same power level (Fig. 5.9) reduces the temper ature gradient near the entrance and pushes the peak tem perature farther down stream. The peak temperature is now slightly lower at 3095 K but the average Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 150 0.5 m 2000K Temperature Contours 2500K 3000K 1m 3250K 0.0 l.E -2 • E-3 Mach Number Contours Figure 5.8 Temperature and Mach number contours for P = 5 MW, m = 0.032 kg/s. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 151 0.5 m { 20001? Temperature Contours 2500K 1m 3000K 0.0 • E-3 \1 .5 E -3 Mach Number Contours Figure 5.9 Temperature and Mach number contours for P = 5 MW, m = 0.128 kg/s. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 152 0.2 m 1500K 2000K 2500K 0.0 0.4 m Figure 5.10 Temperature contours for P = 1 MW, m = 0.021 kg/s. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 153 0.75 m 2000K 300 0.0 1.4 m Figure 5.11 Temperature contours for P = 10 MW, m = 0.288 kg/s. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 154 tem perature at the throat, 2271 K, is nearly the same. The velocities are now significantly higher. Though the temperature field is not significantly different for this higher flow rate, there is a substantial increase in the overall efficiency as is discussed later. The effects of size scale-up on the temperature contours are given in Figs. 5.10 and 5.11 for the smaller (1 MW) and the larger (10 MW) power levels respectively. The mass flows and throat areas chosen for these two cases were determined by maximizing the net amount of power transferred to the gas making them analogous to the high flow rate 5 MW results shown in Fig. 5.9. The effect of increasing the size (optical depth) of the engine has little effect on the overall shape of the tem perature contours. However, the larger sizes result in more net energy transfer to the fluid as is noted later. 5.3.3 Efficiency o f Direct Absorption We have defined two efficiencies - collection efficiency and overall thermal effi ciencies. In this section, we determine these two efficiencies for the thruster config urations described above. These results for the three power sizes are summarized in Table 5.2. Refering to the results in Table 5.2, we first consider the low mass flow rate (m = 0.032 kg/s) condition for the intermediate (5 MW) power case. About 4.3 MW of the incident beam lie in the spectral range of the alkali mixture; 2.7 MW of this energy are absorbed directly by the gas, the remainder along with the radiation th at lies outside the spectral range is incident on the walls. A significant portion of the absorbed energy is later reradiated to the walls (1.13 MW) or lost through the inlet (0.9 MW) resulting in a net collection efficiency of only 14%. At this low flow rate, the working fluid can absorb (by regenerative heating) only a small portion of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 155 Table 5.2 Energy balances for the three thruster sizes (P = l,5 and 10 MW) for different mass flow rates. P (MW) - Incident power in solar beam Qaba (MW) - Power directly absorbed by gas Qwall (MW) - Reradiation to the walls. Qiose (MW) - Reradiation loss through inlet window 77c - Collection efficiency Tj 0 - Overall thermal efficiency p m, kg/s Qaba Qwall Qloaa Vc Vo 1 0.325 0.115 0.055 0.064 0.251 0.087 0.350 0.160 0.095 1 0.0025 0.01 0.021 0.410 0.035 0.010 0.454 0.692 5 0.032 2.700 1.130 0.970 0.137 0.204 5 0.064 0.128 0.865 0.355 0.750 0.225 0.268 0.541 0.414 5 2.770 2.920 10 10 10 0.036 0.144 5.790 2.705 1.740 0.570 2.430 0.082 0.117 1.410 0.380 0.341 0.653 0.481 0.928 1 0.288 5.880 6.170 0.35 0.805 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 156 the wall flux while the rest must be rejected as waste heat. The resulting overall thermal efficiency at this low mass flow rate is therefore very low at 20%. Increasing the mass flow rate in the thruster to 0.128 kg/s (by increasing the nozzle throat area) while fixing the input power improves collection efficiency be cause of reduced reradiation losses. But, the increased mass flow rate provides a much higher overall thermal efficiency because it enables more effective use of the wall flux through regenerative heating. At this condition, the amount of energy absorbed by the gas is marginally higher at 2.9 MW, but the faster moving gas re duces the reradiated energy to 0.35 MW, with a similar reduction in losses through the inlet to 0.22 MW. This increases the collection efficiency to 54 %. This in creased mass flow also closely balances the heat load on the walls providing better utilization of the energy by regenerative preheating. This results in a significant improvement in the overall efficiency to 80%. Further increases in the mass flow above this level would reduce the energy load to the walls even further so that there would not be enough energy to maintain the incoming gas mixture at 1000 K. Thus, for all three power levels, we limited our maximum mass flow rate to the value at which the radiant flux on the walls would bring the inlet mixture to 1000 K. For the 1 MW engine size with the mass flow rate chosen to balance the radi ation loading on the wall, the decreased optical depth causes a smaller percentage of the incident power to be absorbed directly by the gas compared to the 5 MW engine. The collection efficiency at the optimum mass flow rate (0.021 kg/s) is reduced to 45 % (compared to 54% for the 5 MW case) while the overall efficiency is 69% (compared to 80%). The energy balance for the 10 MW engine with the mass flow rate optimized for maximum efficiency provides slightly better performance than the 5 MW engine. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 157 The longer absorption path lengths in the working fluid of this larger engine allow a larger percentage of the beam to be absorbed directly compared to the two smaller engines. As a result, both the collection efficiency (65%) and the overall efficiency (above 90%) are higher for this case (m = 0.288 kg/s). One would certainly expect the overall efficiency of direct absorption to increase as their size increases, and these results give an indication of the rate at which they improve. 5.3.4 Comparison to One-Dimensional Predictions In this section, we contrast our two-dimensional predictions with previous one dimensional estimates for direct solar thermal propulsion [44]. The findings of the one-dimensional analyses have been presented in terms of the variation of the collection efficiency, 7]c, versus the peak temperature. For comparison purposes, we also present our collection efficiencies in this mannner. One difference is that we use the average tem perature rather than the peak temperature. For the two-dimensional problem, the average temperature, which is the more fundamental measure of rocket performance, is considerably lower than the peak temperature. Table 5.3 lists out the peak and average temperatures for different cases. By contrast, for the one dimensional case, the peak and average temperatures are identical. The results of the two-dimensional analyses for the 5 MW engine (taken from Tables 5.2 and 5.3) are shown on Fig. 5.12 along with the previously published one-dimensional results. The qualitative difference between the two sets of results is immediately apparent from this figure. The one-dimensional results lie essentially horizontal while the two-dimensional results run vertical. The peak tem perature for the one-dimensional solution varies over a wide range from 0 to 3500 K; each tem perature corresponds to a unique mas6 flux or eigenvalue. Increasing the mass flux Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 158 Table 5.3 Performance of the direct solar thruster for different sizes, different sizes, mass flow rates and chamber pressure. (Rt is the radius of the nozzle throat). Pressure=3 atm P, MW m, kg/s i2«, mm 2ma*i K TavyK Th, N I spy 6 1 0.0025 3 3046 2555 21 543 1 0.01 6.5 2974 2498 84 537 1 0.021 9 2827 2337 164 520 5 0.016 0.064 8 16 3295 3262 2273 2324 130 524 524 0.128 22 3095 2271 1035 522 10 0.036 11 3320 2192 287 515 10 10 0.144 23 2260 1177 527 0.288 33 3280 3127 2231 2325 521 5 5 528 Pressure=10 atm 5 5 0.092 0.137 10 12.5 3169 3050 2185 2137 731 1075 650 641 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 159 1-D Results .0 EFFICIENCY CR=2 x 104 CR=1 x 10 0.8 2-D Results 0.6 COLLECTION CR=2 x 104 .4 C R -1 x 10 2 0.0 1000 2000 3000 4000 AVERAGE GAS TEMPERATURE Figure 5.12 Collection efficiency vs. average gas temperature. Solid lines show 1-D results (from Rault [44]). Dashed lines show present 2-D results. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 160 decreases the temperature and increases the collection efficiency (by reducing rera diation losses). By contrast, the two-dimensional solution shows th at the average temperatures lie in a very narrow band around 2000 to 2500 K; increasing the mass flow rates results in higher collection efficiencies as noted earlier. We also note that the peak temperature computed for the two-dimensional problem is commensurate with the peak (average) temperatures predicted by the one-dimensional analyses, but the average temperature is lower. The two-dimensional collection efficiencies are also considerably below those obtained for the one-dimensional problem. Improved efficiencies may be obtained by using larger engine sizes as noted earlier. The above result demonstrates the inapplicability of one-dimensional models to the study of the solar propulsion problem; temperatures and collection efficiencies are poorly predicted and completely misleading trends are indicated. This is because the one-dimensional analysis omits much of the basic physics of the problem. The overall efficiency, rfc,, is a more appropriate measure of the effectiveness of solar propulsion system. Figure 5.13 shows the variation of the rj0 with thrust. The thrust levels for the three power levels are listed on Table 5.3. Because average tem perature is essentially independent of mass flow (as noted in Fig. 5.12), thrust is a direct function of the mass flow. Thus, in Fig. 5.13, the overall efficiency increases with thrust until the optimum mass flow condition is reached; clearly, the maximum thrust for each thruster size represents the most effective energy utilization. 5.3.5 Solar Thruster Performance Estimates Hitherto, we have discussed the efficiency of energy transfer from the solar radiation to the propellant. We now turn our attention to the second major issue in the feasibility of the direct solar concept — its performance potential. The performance of the thruster has been estimated based on ideal vacuum expansion Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. OVERALL EFFICIENCY 161 0 . 8 ---- 0 . 6 ----- 1 MW 10 MW 5 MW 0.0 THRUST(N) Figure 5.13 Overall efficiency vs. thrust for three engine sizes P = 1, 5 ,1 0 MW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 162 from the throat. Table 5.3 shows the thrust and specific impulse of the three thruster sizes while . 5.14 shows specific impulse plotted against thrust. For the 1 MW engine, thrust capability ranges from 20 to 160 N; at 5 MW, it ranges from 130 to 1000 N; and at 10 MW, it ranges from 300 to 2300 N. Because the average temperature varies only weakly with input power or mass flow rate, the specific impulse (Isp) is nearly the same (about 500 s) for all cases. This comparatively low Isp is a consequence of the high molecular weight (5) of the working fluid. The molecular weight may be decreased by reducing the alkali fraction in the working fluid; however, this will result in poorer absorption charac teristics of the working fluid. Another way of reducing the molecular weight of the propellant is by increasing the chamber pressure and by reducing the alkali mole fraction in such a manner that the number density of alkali atoms (and, hence, the absorption characteristics) is unchanged. For example, by increasing the cham ber pressure from 3 to 10 atm, the total alkali mole fraction can be reduced to 1.5% (0.5% of each alkali) without changing the absorption characteristics and the molecular weight may be reduced to 3. Computations for the 5 MW thruster under these conditions show that the efficiencies and thrust are essentially unchanged but the specific impulse is significantly higher at about 650s (see Fig. 5.14). Further increases in chamber pressure will provide similar improvement in performance. The specific impulse predictions for the 3 atm condition are approximately equal to those for the indirect solar propulsion [15,16], while those for the 10 atm condition are significantly better. In this regard, it should be noted that our twodimensional predictions for direct solar propulsion are considerably lower than previ ous one-dimensional estimates, and it is likely that more in-depth multi-dimensional analyses of indirect solar propulsion may likewise reduce the estimated performance Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 163 800 SPECIFIC IMPULSE (s) 10 atm 6 0 0 ------ 3 atm 5 MW 10 MW 400 200 THRUST (N) Figure 5.14 Performance of the direct solar thruster - specific impulse vs. thrust level for three power levels, P — 1, 5, 10 MW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 164 of that configuration. It is, however, clear that successful implementation of the di rect concept will involve optimization between thruster size, propellant mass flow rate and chamber pressure. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 165 CHAPTER 0 INVESTIGATION OF MICROWAVE ABSO RPTIO N IN GASES 6.1 Resonant Cavity Microwave Plasmas In this chapter, we switch our attention to the second problem of interest, which is microwave propulsion. Just as in the case of solar propulsion the ultimate objective with microwave propulsion is to address the feasiblity issues of the con cept. In the present study, we focus our attention on the absorption of microwave energy in flowing gases and reserve the detailed evaluation of practical thruster configurations for future study. A physical model involving the coupled solution of the Maxwell equations for the microwaves and the Navier-Stokes equations for the fluid is employed in the study. Experimental programs evaluating different absorp tion schemes [23-31] such as resonant cavity plasmas, waveguide-heated plasmas and co-axial applicator plasmas are presently underway. The present analysis is for resonant cavity plasmas and the study is closely coupled to a companion experimen tal program [27,28]. The experimental data is used to validate the computational results and, at the same time, the numerical results are used to provide a better physical understanding of the experiments. Once the physical model is completely validated, the results may be extended to larger sizes and power levels. Practical geometric configurations may also be evaluated. Furthermore, it is quite straight forward to extend the model to study waveguide-heated plasmas as well. The analytical model was closely coupled to the companion experimental pro gram [27,28]. Figure 1.2 shows a schematic of the experimental set-up. The working Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 166 fluid (helium and nitrogen have been used in the experiments so far) flows from top to bottom in the sphere-cylindrical tube while the microwave energy is fed into the surrounding cylindrical cavity by means of a coupling probe. The frequency of the microwave source is 2.45 GHz. Powers up to 3 kW and gas pressures up to 5 atm are obtainable in the experimental apparatus. The experiments using the T M ou mode for the microwaves in the cavity have determined the stability boundaries of microwave-sustained helium plasmas for a range of pressures (0.4 to 3 atm) and flow rates (up to 2 m /s and 3 x 10-4 kg/s). Coupling efficiencies (the percent of the input microwave power that is coupled to the plasma) of up to 80% and overall thermal efficiencies (including radiative and diffusive losses) of up to 40% have been obtained. The reason for the relatively low coupling efficiencies (other experiments [24,25] have reported coupling efficiencies up to 99%) is the inclusion of cavity wall losses in these measurements. Electron temperature measurements have yielded peak plasma temperatures in the range of 10,200 to 10,900 K which are relatively insensitive to variations in flow conditions. Most of these data were limited to pow ers of about 500 W because, at higher powers, the plasma tended to arc with the quartz walls of the tube. This instability of the plasma is thought to be due to the asymmetric location of the coupling probe. Further, even at lower powers, tuning the cavity to minimize reflected power resulted in the same difficulty. This is the second reason for the relatively low coupling efficiencies obtained. In order to stabilize the plasma at high power levels, a bluff body was intro duced into the flow. Recent experiments [101] have demonstrated th at the bluff body is an excellent plasma stabilization device and plasmas at powers up to 2 kW could be sustained. Plasma instability problems were no longer experienced when the cavity was tuned and coupling efficiencies up to 100% were measured for pressures of 1 to 3 atm. These experiments with the bluff body also included an Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 167 orifice plate in order to enable thrust and specific impulse measurements. Further, the TM on mode was employed in order to locate the plasma near the orifice. This lower order mode also has the advantage of reducing the physical size of the de vice. Currently, detailed temperature measurements including radial temperature profiles of the plasma are being conducted. Nozzle exit velocity measurements are also planned to determine thruster performance. As mentioned earlier, the calculations followed the experimental set-ups as closely as possible. Both the sphere-cylinder configuration and the bluff-body con figuration have been modelled. Representative grids for these geometries were pre sented in Chapter 4 (Figs. 4.1 and 4.7). In addition to the experimental configura tions, a straight duct was also studied in order to assess the effect of flow geometry on the discharge characteristics. Most of the calculations employed a fluid dynamic grid of 101 x 61. Calculations with finer grids, such as 201 x 101 gave the same basic solutions. For the Maxwell solutions of the microwave field in the cavity, a standard rectangular grid was employed. Information transfer between the microwave grid and the fluids grid was achieved by a two-dimensional linear interpolation proce dure. To resolve the strong gradients in the electromagnetic field in the region of the plasma, a reasonably fine grid was employed. Most of the calculations shown here used a 121 x 121 grid. A coarser 61 x 61 grid was employed in some of the calculations especially those involving the TM on mode. The electromagnetic field solutions were somewhat more sensitive to grid size especially with regard to nu merical convergence. In the following sections, various results showing the characteristics of mi crowave plasmas are discussed. First, we present some basic computational results for different geometries such as the straight duct, the sphere-cylinder and the bluff Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 168 body configuration. Next, the validation of the model by comparison with experi ments is discussed. In the final section, we undertake some parametric studies such as the effect of microwave power, gas pressure and gas velocity on the plasma char acteristics. In all cases, the gas employed was helium. The equilibrium viscosity and thermal conductivity were obtained from Ref. 74. The electrical conductivity was calculated as discussed in Chapter 2. The effects of buoyancy and Lorentz forces were included; however, radiation losses from the plasma were not accounted for. 6.2 Basic Computational Results 6.2.1 Microwave Plasmas in a Straight Duct The radius of the duct was taken as 0.05 m which is the outer radius of the spherical bulb in the experimental geometry; and the length of the duct was taken as 0.2 m. The flow conditions employed are similar to experimental conditions [27,28]. The T M 0 1 2 mode was employed for the microwaves in the surrounding cavity. Corresponding to this mode, for a cavity radius of 0.089 m, the resonant length is 0.148 m. A gas pressure of 1 atm , an inlet temperature of 1000 K and an inlet gas velocity of 0.11 m /s were employed for the representative case. The corresponding mass flow rate is 4.2 x 10-5 kg/s. The high inlet gas tem perature was chosen to improve convergence of the numerical scheme and does not materially affect the temperature profile of the plasma. The rate of convergence of the solution are shown in Figs. 6.1a and 6.1b. In Fig. 6.1a, a relatively coarse EM grid of 41 x 41 was employed. The convergence rate is seen to be linear on the semi-logarithmic plot for the first two or three orders of magnitude reduction in error. After that, the solution stops converging. This Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 169 < az> Ui ec UJ £ o -1_6 ---- -3 0 100 200 300 400 500 600 ITERATION STEPS (a) MORII OF RESIDUALS o — r— L-2 - 6 ----- -3 - 10 0 200 400 600 300 1000 ITERATION STEPS (b) Figure 6.1 Convergence of the coupled Navier-Stokes and Maxwell equations. Fluids grid 101 x 51, EM grid - (a) 41 x 41, (b) 121 x 121. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 170 difficulty in convergence is due to the energy equation and may be attributed to in sufficient resolution of the electric field in the region of the plasma. To demonstrate this, Fig. 6.1b shows the convergence for the same case except that now a finer EM grid of 121 x 121 was used. The convergence is now seen to be consistently linear down to about five orders of magnitude reduction in error. This result shows that the steep gradients in the electric field require a reasonably fine grid resolution. Contours of temperature and velocity in the straight duct along with stream lines for the representative case described above are shown in Fig. 6.2. The outline of the cylindrical cavity containing the microwaves and the co-axial waveguide by which the microwaves are introduced into the cavity are also shown. The microwave input power was 2 kW. The resulting solution indicates that 1.84 kW of the input power are absorbed, corresponding to a coupling (or absorption) efficiency of 92%. The remainder of the power is reflected back to the source. The length of the cav ity may be adjusted to maximize the efficiency as described in Chapter 2. In the present calculation, however, the cavity length was not optimized. The tem perature contours in Fig. 6.2 show a peak temperature of 9300K. The undistorted TM 012 electric field has three nodes any of which could serve as the site for the plasma. In the present calculations, the plasma is seen to form at the lower node. The velocity profiles in Fig. 6.2 show the input profile which corresponds to fully developed flow rapidly deforms as it approaches the plasma. The centerline velocity initially decelerates and then accelerates to a maximum at the center of the plasma. The corresponding streamlines show the flow turns slightly outward as it passes the plasma, indicating that the microwave plasma causes the traditional blockage effect that is observed in laser plasmas [37]. This blockage is a result of reduced gas densities at the higher temperatures encountered in the plasma. The axial and radial components of the electric field are shown in Fig. 6.3. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 171 FLOW IPO K [OOO 0.2 8500 IT* 9000 K / ^ 0 .5 m/s TEMPERATURE V ELOCITY STREAMLINES Figure 6.2 Representative flowfield solutions for a floating plasma in a straight duct. P0 = 1 atm , m = 4.2 x 10“ 5 kg/s, Re = 13, Pinc = 2 kW, Pref = 0.16 kW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 172 NORMALIZED AX IAL ELECTRIC RADIAL E L E C T R IC F IE LD F IE L D Figure 6.3 Distorted electric fields for the straight duct solution in the TM 0i2 mode. The electric field is normalized by E ref = 25000 V /m . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 173 As seen in Chapter 2, the presence of the conducting gas, distorts the electric field significantly from its no-loss solution. The steep gradients in the plasma referred to earlier are seen at the region of the lower node. The location of the peak temperature of the plasma in these plots is marked by a darkened square. The electric field is seen to be a minimum inside the plasma. This is because of the traditional skin effect whereby most of the microwave energy is absorbed at the surface of the plasma itself and very little energy actually penetrates into the plasma. We also note that even though the TM 012 mode solution is distorted beyond recognition in the region of the plasma, the general features of the mode solution are still evident at regions remote from the plasma. For the TMui 2 mode, there are three likely sites at which the plasma may be sustained. In this particular case, the plasma forms at the bottom node. Ex perimental results with the sphere-cylinder geometry consistently show the plasma forming at the middle node. Even when the plasma formed at one of the other locations, it was observed to be unstable and the plasma always reverted to the middle node. There may be several factors involved in the selection of the node by the plasma such as the fluid dynamics in the discharge tube, the nature of the mi crowave energy field and the initial conditions of the computation. One possibility is th at the plasma tends to form at the end where the microwave inlet is located (as it does in this case). This follows from the tendency of plasmas to propagate toward the source of energy - which has been observed both in laser and microwave plasmas. However, in a resonant cavity, the microwaves bounce back and forth and a preferential propagation direction may not exist. The dynamics of the location of the plasma within the dicharge tube is therefore not obvious. Parametric studies of plasmas in a straight duct, which are presented in a later section, reveal that the plasma may select any one of the three nodes depending on the initial condition. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 174 6.2.2 Plasmas in Sphere-Cylinder Configurations The geometry shown in Fig. 4.1 was used to investigate the plasma in the ex perimental sphere-cylinder configuration. Figure 6.4 shows a representative solution for this geometry. The diameter of the sphere is 0.1 m while the diameter of the inlet pipe is 0.025 m. The overall length of the sphere-cylinder combination is 0.2 m. These calculations were for a pressure of 1 atm and an inlet gas velocity of 0.34 m /s. The corresponding mass flow rate is 8 x 10-8 kg/s. The input power was 3 kW and the absorbed power was 2.2 kW giving a coupling efficiency of 73%. Figure 6.4 shows temperature, velocity and streamline contours. The peak temperature is 9200 K and the plasma is seen to form at the center of the sphere which corresponds to the location of the middle node of the T M 0 1 2 solution. The predicted plasma loca tion agrees with the experimentally observed location. Parametric studies revealed that the plasma forms at the middle node over a wide variety of conditions, which is in agreement with experimental findings. These results are discussed in detail in a later section. Figure 6.5 shows the electric fields corresponding to this solution. The general features of the E-field are similar to those in Fig. 6.3 except that the plasma is located at the middle node. For this case, the inlet coaxial waveguide was located at the bottom end of the cavity near the peripheral wall. For the sphere-cylider, the plasma shows a preference for the middle node while for the straight duct all three nodes may serve as a site for the plasma. Clearly, the flow geometry plays an important role in determining plasma location. For the sphere-cylider, the flow tube is constricted at the location of the two end nodes of the cavity; this causes the flow in these regions to be faster than at the center of the bulb. The plasma is thus more stable at the center than at the two end nodes. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. nax 92G0K Temperature Contours V e l o c i t y Contours Figure 6.4 Representative solution of the sphere-cylinder configuration. P — 1 atm , m = 1 x 10-5 kg/s, Re=10, P,nc = 3 kW. Ol CD ission of the copyright owner. 1.2H Further reproduction prohibited without permission. A x ia l E l e c t r i c F i e l d Radial E l e c t r i c F i e l d Figure 6.5 Distorted electric field solutions for the TM 012 mode. The electric field lines are normalized by E ref = 25000 V /m . 177 6.2.3 Bluff Body Stabilized Microwave Plasmas The second set of the resonant cavity experiments at Penn State have employed a bluff body in the flow as a plasma stabilization device. Further, experiments with waveguide-heated plasmas also employ a bluff body in order to keep the plasma stationary. The numerical investigation of bluff body plasmas would be useful to understand the physics and dynamics of both these devices. In the present analysis, only resonant cavity plasmas are studied but the extension to waveguide plasmas is straightforward. In an earlier chapter, we looked at some cold flow results in order to characterize the recirculating region in the wake of a bluff body. In this section, coupled Navier- Stokes and Maxwell solutions for the bluff body stabilized plasma are presented. The same grid specified in Fig. 4.7 was used in these calculations as well. The length of the flow tube is 0.2 m and its diameter is 0.1 m. The diameter of the bluff body is 0.04 m. In keeping with the experiments, the microwave cavity was set up for the TM ou mode instead of the TM012 mode. Accordingly, the resonant length of the cavity is now 0.074 m or half the old value. The convergence history for a representative case of power 3 kW, inlet velocity 3 m /s and pressure of 3 atm is shown in Fig. 6.6. The corresponding inlet Reynolds number is 1000. The convergence rate of the four primary dependent variables (that is, the vector Q defined in Chapter 4) is linear on the semi-logarithmic plot. The error magnitude is decreased by about 6 orders of magnitude over 600 iteration steps. Considering the difficult nature of the flowfield and the strong tem perature gradients in the plasma, the convergence rate is very good. Just as in the cold flow cases, convergence improved with increase in Reynolds number and deteriorated with decrease in Reynolds number. The detailed flowfield and E-field solutions for this case is shown in Fig. 6.7. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 178 - 3 ---- -10 0 100 200 300 400 500 600 ITERATIONS Figure 6.6 Convergence of coupled Navier-Stokes and Maxwell equations for the bluff body geometry using the TM qu mode. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TEMPERATURE STREAMLINES 179 -6 .E -9 0. 0 P“ 3 a t m Re«1000 AXIAL ELECTRIC FIELD 3.8k J i J i l l RADIAL ELECTRIC FIELD T J 0. Figure 6.7 Representative flowfield and electric field solutions for the bluff body geometry. P0 = 3 atm, P = 3 kW, Re = 1000. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 180 Temperature and streamline contours are shown in Fig. 6.7a while the axial and radial electric fields are shown in Fig. 6.7b. The temperature contours show that the plasma is almost attached to the bluff body. This feature has been noticed in experiments [101] as well. The peak tem perature is 9800K and occurs slightly offcenter leading to a slightly toroidal shaped plasma. The temperature distribution within the plasma is quite uniform. This is again due to the skin effect of the plasma. The streamlines show a recirculating region that is bigger than the corresponding cold flow result (Fig. 4.9d). However, in terms of the recirculating mass flow rate, the recirculation is weaker with the hot plasma. This result may be explained by the significant reduction in density of the hot gases and the increase in local kinematic viscosity, both of which cause a decrease in the local flow Reynolds number. This weakening of the recirculation region by heat addition have also been obtained in the calculations of a bluff body combustion flame stabilizer [102]. The electric field components are shown in Fig. 6.7b. Apart from the fact, that the fields now show the T M qu mode shape, the field lines are similar to the results in the last section. The fields are distorted in the region of the plasma. There are steep gradients at the plasma surface and the electric fields are nearly zero within the plasma (a consequence of the skin effect). The general shape of the TM 0n mode is evident only in regions remote from the plasma. The coupling efficiency of this case is 78% which means that 2.22 kW of the 3 kW that is incident on the plasma is absorbed by the plasma. In order to reduce the reflected component and maximize the power coupled to the plasma, the cavity may be tuned by adjusting the cavity length. The optimization of cavity length results in coupling efficiencies as high as 99% as shown in the following section. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 181 0.3 M odel Validation W ith Experiments In our analysis, the computations have been tailored to match experimental conditions of the companion program in order to enable validation of the numerical model. It would be desirable to make detailed comparisons of temperature and ve locity profiles within the plasma; however, such detailed measurements are not yet available. Recently, radial temperature profiles were measured [101]; however, spa tial resolution is still inadequate and the Abel inversion technique used to generate these radial profiles [101] is known to be inaccurate at the edge of the plasma. At the center of the plasma, where the measurements are reliable, the temperature is uniform. Thus, model validation is limited to comparing the peak temperatures of the plasma. The other global parameter that is available for comparison is the cou pling efficiency. In addition to these, numerical predictions agree qualitatively with experimental observations such as plasma location, shape and size, the influence of discharge pressure, incident power and gas velocity as well as of cavity tuning. The one result th at is not predicted well by the numerical model is the stability boundary or the threshold power—which is the minimum power required to sustain a stable plasma at different pressures. Numerical predictions are consistently higher than observed values; possible reasons for this discrepancy are discussed later. In the present section, experimental data for peak temperature and coupling efficiency are compared with computational predictions; other parametric studies are considered in the following section. 6.3.1 Coupling Efficiency Coupling or absorption efficiency is the ratio of power absorbed by the plasma to the power incident on the plasma. Experiments and numerical analyses show Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 182 that coupling efficiencies up to 99% are possible for resonant cavity plasmas. This is because the resonant cavity configuration allows multiple passes of the microwaves through the plasma thus enabling almost complete absorption of the energy. Figure 6.8 shows experimental results along with computational results as a function of discharge pressure. Other parameters such as mass flow rate or input power do not have any effect on the coupling efficiency. The first of the two sets of experimental data shown on Fig. 6.8 is for the sphere-cylinder configuration (from Ref. 28) using the TM012 mode. Corresponding to this, the maximum coupling efficiency is only about 80%. As mentioned earlier, this lower value is because of cavity wall losses and difficulties in tuning the cavity. The plot also shows a definite trend in the magnitude of coupling efficiency. The maximum coupling efficiencies occur over a range of moderate pressures (1 to 3 atm). At higher and lower pressures, the efficiency starts to tail off. These trends may be explained by considering the collision processes that the electrons undergo. At low pressures, the collision frequency of electrons is low and much of the energy absorbed by the electrons is simply returned to the electric field during the second half phase of the oscillation. At higher pressures, the electron collision frequency is so high that the electrons do not get enough time between collisions to pick up energy from the electromagnetic field. Eventually, at very low and very high pressures, the plasma will be extinguished altogether because of increasingly poor absorption efficiencies. The second set of experimental data (from Ref. 101) is for the bluff body con figuration using the TM on mode. In this case, wall losses are not accounted for in the experimental measurements and the coupling efficiency is determined simply by using the measured values of the incident and reflected powers. Further, sta bilization of the plasma by the bluff body enables fine tuning of the cavity length Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. COUPLING e - -e—-O Expt. s p h e r e - c y l i n d e r &-- 9 - -a Expt. b l u f f body $ - - 3 - -0 Comp. s p h e r e -c y l In d e r * 2 0 Comp. b l u f f body * 1— i— i— i— j— i— i— i— i— j— i— i— i— i— |— i— i— i— i— j— i— i— i— i— | 0 1 2 3 4 5 X I05 PRESSURE, N /m 2 Figure 6.8 Variation of coupling efficiency with discharge pressure. Comparison of computations with experiments (from Balaam [101]). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 184 to obtain maximum power absorption. As a result, these experimental coupling efficiencies are as high as 100%. For this case, measurements for sub atmospher ic pressures are not available. Over the pressure regime measured, the coupling efficiency is nearly constant but it should drop off at higher and lower pressures. The computational results in Fig. 6.8 are for both the sphere-cylinder and the bluff-body configurations. The sphere-cylinder results were obtained by carefully adjusting the cavity length for maximum efficiency as discussed earlier. For the bluff body, a coarser tuning procedure was adopted in order to save computation al expense. In both cases, as a simplification, the wall was taken to be perfectly conducting (and lossless). The situation is, therefore, similar to the bluff-body ex periments [101]. Indeed, the finely-tuned predictions of coupling efficiency match the bluff body experimental results very closely. The coupling efficiency reach es a maximum of about 100% at intermediate pressures (about 1 to 3 atm) and decreases slightly at lower and higher pressures. Figure 6.8 also shows coupling efficiencies from the bluff body computations. These coupling efficiencies are seen to be consistently around 90 to 95% which is somewhat lower than the spherecylinder predictions. Finer tuning of the cavity length for this case should bring these coupling efficiencies up to 99% as well. It should be noted that once the cavity length was adjusted to maximize ab sorbed power for a given case, further adjustments of the length were typically not necessary. In the experiments, the change in the length for different cases amounted to only about 1-5 mm. The computational results showed a similar insensitivity to the cavity length for different cases once the length was optimized. This is an im portant result from a practical viewpoint since the adjustment of cavity length on board the spacecraft would increase the complexity and weight of the system. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 185 6.3.2 Peak Temperatures Experimental measurements and computational predictions of plasma peak tem perature are considered next. Figure 6.9 shows the peak temperatures plot ted against the discharge pressure. Neither experiments nor computations show a significant dependence of the peak temperature on the pressure. There are several sets of experimental and computational results shown. Overall agreement between the numerical results and measurements is reasonably good although the numerical peak temperatures are somewhat underpredicted. The experimental data are from spectroscopic measurements of the electron tem perature of the plasma. As mentioned earlier, radial temperature measurements were attem pted but the spectroscopic data are reliable only for the plasma core where the temperature is quite uniform. Accordingly, only the peak tem perature is used in the comparison study. Two sets of experimental data are plotted on Fig. 6.9. The first set of data is for the sphere-cylinder configuration [28] while the second set is for the bluff-body set-up [101]. There is a marked difference between the two sets of data. For the sphere-cylinder, the peak plasma temperature is constant at about 10,000K, while for the bluff-body, the peak temperature is constant at about 12,000K. This difference is thought to be due to the use of better spectroscopic equipment for the bluff-body measurements [101]. For both sets of data, the electron temperature is quite independent of the gas pressure. Measurements at different flow rates and power levels revealed th at the electron temperatures were insensitive to these parameters as well. Similar results have been reported by other researchers as well. Batal et al. studied argon dis charges [103]. They observed a Te of 10,000K which was insensitive to changes in Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 186 PEAK TEMPERATURE ,K “•000 • □ -e- G— 0- 10000 - B B -S -e — ★ <$■ if. 5 000 ■ 0 B -e 0 —0 —0 Expt. s p h e r e - c y l i n d e r S—B -B Expt. b l u f f body 0 Comp, s p h e r e - c y l i n d e r 0 X X Comp, s t r a i g h t duct $ Comp, b l u f f body A 0 --------- 1— i— i— i— |— i— i— i— i— |— i— i— i— i— |— i— i— i— i— |— i— i— i— i | 0.0 0.5 1.0 B- 1. 5 2.0 2.5 i i i i j 3.0 X105 PRESSURE, N /m 2 Figure 6.9 Variation of peak plasma temperature with discharge pressure. Comparison of computations with experiments (from Balaam [101]). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 187 input power, mass flow rate and pressure. Dymshits and Koretskii [104] reinterpret ed Kapitza’s data [20] for hydrogen plasmas and determined Te to be approximately 10,OOOK. Computational results for microwave plasmas sustained in three different flow configurations—the straight duct, the experimental sphere-cylinder and the bluffbody—are also shown in Fig. 6.9. Flow conditions were matched as closely as possible to the experiments. Since the analysis assumes local thermal equilibrium (LTE), the electron temperature is assumed to be the same as the gas tempera ture. The three sets of data points lie very close to each other indicating little dependence on the flow geometry. The peak gas temperatures range from about 9000 to 10,000K. There is a slight increase in the temperature with gas pressure which was not apparent in the experiments. In general, these computational predic tions are in reasonable agreement with experimental data. The comparison to the sphere-cylinder measurements is excellent; however, the bluff-body measurements are about 20% higher than the numerical predictions. This is somewhat discourag ing since the bluff-body measurements are thought to be more accurate. We note, however, that the electron temperature measurements reported elsewhere [103,104] are about 10,000K. Further, this discrepancy may partly be attributed to the 10% error bounds of the spectroscopic measurements. The underprediction of the peak gas temperature by the computations may also have physical significance. It is possible that the LTE assumption may not be accurate; as a result, the electron temperature may be slightly higher than the gas temperature. Experimental verification of the equilibrium assumption was attem pt ed by measuring the gas temperature by an independent spectroscopic procedure; however, the results of this study proved inconclusive [101]. A more detailed study is necessary to investigate non-equilibrium effects in the plasma. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 188 6.4 Effect o f Parameter Variations on Microwave Plasmas The overall characteristics of microwave plasmas such as plasma location, size and shape as well as temperature profiles and gas flowfields within the plasma are influenced by many parameters. Some of the parameters of interest are microwave power, gas velocity and gas pressure. Experimental investiagtions have concentrated on the effect of discharge pressure and power on the plasma characteristics. How ever, experiments have been limited to low powers and low-to-moderate pressures. Practiced thrusters range from low power (1-5 kW) devices for space station-keeping and moderate power devices (100-500 kW) for orbit trasfer. Numerical computa tions may be used to study plasmas at higher powers and pressures and thus enable the analysis of both kinds of devices. In this section, as a first step to this inves tigation, we present detailed flowfield solutions for a variety of flow conditions and power levels. The sphere-cylinder configuration as well as the bluff body configu ration are used in these calculations. Experimental measurements of such detailed tem perature profiles and flowfields are not yet available and, hence, direct compari son of these data with experiment is not possible. Most of the predicted trends are, however, in qualitative agreement with experimental observations. 6.4.1 Effect o f Variations in Input Power The effect of increasing the incident microwave power on the microwave dis charges is shown in Figs. 6.10 - 6.13. Plasmas sustained in the straight duct, the sphere-cylinder as well as the bluff body set-up are examined. The results for the straight duct are shown in Fig. 6.10 for a series of powers ranging from 1.5 kW to 4 kW. In all cases, the discharge pressure was maintained at 1 atm and the inlet gas velocity 0.1 m /s corresponding to a mass flow rate of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 189 CM 6 "3 u 3 60 *3 P uo d 3 (O cd § to ■ 44 Pi w PL, S w H cti a va> (0 o & 60 t? 2 c •3 Vi 0 1 ° ~X <N s -* ♦> 1 " w •£ «o V (H 3 60 O IA Reproduced with permission o f the copyright owner. Further reproduction prohibited without permission. 190 TEMPERATURE S 't STREAMLINES Figure 6.11 Effect of input power on plasmas in the sphere-cylindrial tube. P0 = 1 atm, m = 2.6 x 10-6 kg/s. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TEMPERATURE (b) 2 kW 3 kW STREAMLINES Figure 6.12 Effect of input power on bluff body plasmas. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 192 PEAK TEMPERATURE, IC 15000 B -e -B s t r a i g h t duct 0 —0 —0 sphere-cylinder A A ^ b l u f f body 10000 5000 0 1 2 3 4 INCIDENT POW ER, kW Figure 6.13 Effect of input power on peak plasma temperature for straight duct, sphere-cylinder and bluff body geometries. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 193 4.2 x 10-5 kg/s. The 1.5 kW plasma in Fig. 6.10a represents the threshold condition for the plasma. It was not possible to maintain the plasma at powers lower than this threshold power for this pressure. As mentioned earlier, experimental findings also indicate such a power threshold. At this lowest power, the plasma is seen to form at the lowermost node of the TMm 2 cavity similar to the result presented in Section 6.2.1. Increasing the power to 2 kW and 3 kW (Figs. 6.10b and 6.10c respectively) basically results in a larger plasma size. The peak temperature of the plasma does not change very much—in fact, it decreases somewhat. Consequently, the temperature in the central core of the plasma becomes more or less uniform. This is a result of the plasma skin effect discussed earlier. For a power of 4 kW (Fig. 6.10d), there is a dramatic change in the solution. The plasma now jumps to the middle node. This shift in location may be explained as follows. W ith increasing power, the plasma grows in size—not just radially but also axially. At some point in this expansion, the gas in the region of the middle node becomes hot enough th at it becomes electrically conducting. This results in absorption of energy at the middle node. Further increase in power, results in the activation of the uppermost node as well. This is the situation in Fig. 6.10d. The inlet gas breaks down in the region of the upper node and starts absorbing the microwave energy. Further energy addition occurs at the middle node. Thus, the plasma is in the shape of an elongated club. All the incident microwave energy is absorbed by the gas at these two nodes and there is no energy remaining for absorption by the time the gas reaches the bottom node. At higher powers, a triple plasma involving heat addition at all three nodes may be anticipated. Figures 6.10e and 6.10f show the results when the incident power was reduced starting from 4 kW. The corresponding powers are now 3 kW and 2 kW respec tively; all conditions are therefore identical to the cases in Figs. 6.10c and 6.10b Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 194 respectively. However, the solutions are seen to be quite different. Cooling the plas ma shows the plasma resting stably at the middle node instead of at the bottom node as in the earlier results. Thus, multiple solutions are possible for the same set of operating conditions. The plasma may be sustained stably at either the middle or the lower node and the choice of node is dependent only on the initial conditions. This is an important result but experimental verification is not available since the straight duct geometry has not been employed in the experiments. For the spherecylinder, experiments showed the plasma consistently forming at the center node. At the other two locations, the plasma was found to be unstable. This result is probably due to the fluid dynamics of the sphere-cylinder tube. The effect of incident power on the computed sphere-cylinder results are shown in Fig. 6.11. The input powers were varied form 1.5 kW to 4 kW as before. The gas pressure was 1 atm and the inlet velocity was 0.1 m /s. The corresponding mass flow rate was 2.6 x 10-6 kg/s. Fig. 6.11a shows the temperature and streamline contours for a power of 1.5 kW. This corresponds to the threshold power for this case. The plasma forms at the middle node. The streamlines show evidence of a large toroidal recirculating bubble in the region of the plasma. This recirculation is due buoyancy forces in the hot gas (i.e., natural convection). For this case, the buoyancy forces dominate over the convective, diffusive and electromagnetic (Lorentz) momentum fluxes. Increasing the incident power to 2 kW (Fig. 6.11b) results in a larger plasma size. As with the straight duct the plasma temperature does not change very much and the temperature becomes uniform at the central core. The streamlines show that for this higher power case, the recirculating region becomes weaker. Since the gas pressure and velocity are the maintained constant, the relative importance of convection and diffusion does not change. However, at the higher power, the electric field at the surface of the plasma increases and, consequently, so do the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 195 Lorentz forces. This increased importance of Lorentz forces at the higher power is the reason for the weakening of the natural convection bubble. Figures 6.11c and 6.11d show the computed solutions for powers of 3 and 4 kW respectively. The same trend continues with these higher powers—the plasma grows in size and the recirculating region becomes weaker. In fact, for the 4 kW case, evidence of a dual plasma similar to the straight duct calculations is noticed. Further, the recirculating region is completely suppressed at this higher power. The streamlines do turn outward from the plasma due to the blockage effect but flow reversal no longer takes place. In the sphere-cylinder computations, the plasma consistently formed at the middle node—a result that is in agreement with experiments. In fact, in the series of calculations described above, multiple solutions similar to the straight duct case, did not result. Efforts to maintain a plasma at the upper or lower nodes were not successful. This is probably because the location of these nodes is in the narrow cylindrical neck of the set-up, where the gas velocities are considerably higher. These higher gas velocities tended to push the plasma downstream and the plasma tended to be unstable in these sections. This situation is quite likely similar to the conditions in the experimental set-up. Further, the tendency to form dual plasmas at higher powers was also observed in experiments. Similar results for the bluff-body are shown in Fig. 6.12. Results are shown for four different powers—1.5, 2, 3 and 4 kW. Once again, the plasma is seen to grow in size at higher powers. The effect is less marked here than for the free-floating plasma. This is because the size of the bluff body also influences the size of the plasma. The effect of heat addition on the size of the recirculating region is also similar. In fact, the recirculating region in the wake of the bluff body is almost completely supressed at the highest power (4 kW, Fig. 6.12d) shown. It should be Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 196 noted that for the bluff-body case, the recirculation is due to buoyancy as well as the result of the bluff body. In particular, the reciculating region exists even in cold flow (Fig. 4.9) and increases in size with the Reynolds number (or bulk convective flux). This result demonstrates that at higher powers the Lorentz forces outweigh the convective fluxes and buoyancy and completely suppress the flow separation. This effect on the nature of the flow in the wake of the bluff body will play an important role in the diffusive heat transfer from the plasma to the bluff body. Since high heat transfer rates to the bluff body are undesirable, optimal conditions for maximum overall thermal efficiency may be determined. The effect of input power on the peak plasma temperature is shown in Fig. 6.13 for all three geometries under consideration. The temperatures for the three geometries are more or less the same. W ith increasing power, the predicted tem peratures decrease slightly. Note that the decrease in temperature corresponds to an increase in the overall size of the plasma. 6.4.2 Effects of Gas Velocity Gas velocity is an im portant parameter for characterizing microwave plasmas. In a practical thruster, we are interested in determining the behavior of the thruster over a range of flow rates as well as in establishing optimum flow rates for maximum thermal efficiencies. This investigation is intended to see how gas velocity affects the nature of the plasma discharge. The effect of change in gas velocity for a straight duct is shown in Fig. 6.14. Four different inlet gas velocities of 0.017, 0.1, 0.4 and 1 m /s are shown. The pres sure was maintained at 1 atm and the microwave power at 3 kW. The mass flow rates corresponding to the three cases are 7 x 10-6 , 4.2 x 10“ 5, 1.7 x 10~4 and 4.2 x 10-4 kg/s respectively. The corresponding Reynolds numbers are approximately 2, 12, Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TEMPERATURE STREAMLINES Figure 6.14 Effect of gas velocity variation on straight duct plasmas. P0 = 1 atm, P = 3 kW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 198 50 and 125. At the lowest gas velocity (Fig. 6.14a), diffusion, buoyancy and Lorentz forces dominate over convection. The plasma forms at the middle node by virtue of the initial profile assumed for the plasma. The temperature contours are more or less symmetrical above and below the heat addition zone. There is however evidence of heat addition occuring at the upper node as well resulting in a clubshaped plasma. The streamlines show two large recirculation regions - one above and one below the heat addition zone. These are a consequence of the buoyancy effects which are significant at the low flow velocity. At the higher gas velocities, convection becomes increasingly important and dominates over diffusion and buoy ancy. Thermal diffusion of energy upstream of the plasma becomes significantly reduced and the plasma is pushed downstream. Heat addition is then limited to the middle node. Radial thermal diffusion is also reduced resulting in a slimmer and more compact plasma profile. As a further consequence of the convection effect, the temperature contours extend further downstream of the plasma. The streamline contours show that buoyant recirculation no longer takes place and the gas flows through the plasma instead of around it. The blockage effect of the plasma is still evident in the streamline contours; these contours also show the downstream dis placement of the plasma. At higher gas velocities, it was difficult to maintain the plasma. Typically, the plasma was blown off downstream of the heat addition re gion resulting in cool down. Such plasma blow-offs represent an important physical threshold for practical operation of microwave devices. The effect of increasing gas velocity on the peak plasma temperature is shown in Fig. 6.15. Over the gas velocity range considered here, the peak temperature is seen to increase slightly with increasing gas velocity. Note th at the increasing gas temperatures correspond to smaller plasma sizes, which is in agreement with the results in Fig. 6.13. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 199 PEAK TEMPERATURE, K 12500 1 0 0 0 0 ----- 7500 5000 2500 0.00 0.25 0.50 0 .75 1 .00 1 .25 1 .50 INLET VELOCITY, m /s Figure 6.15 Effect of gas velocity on peak plasma temperature for the straight duct geometry. P0 = 1 atm, P = 3 kW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 200 Similar results for the sphere-cylinder geometry is shown in Fig. 6.16. The inlet gas velocities are 0.1, 0.3 and 1 m /s and the corresponding Reynolds numbers are 3, 9 and 30 respectively. The pressure is 1 atm and the power is 3 kW. The results are qualitatively similar to the straight duct results. Increasing the gas velocity pushes the temperature contours further downstream of the plasma and reduces thermal diffusion upstream of the plasma. The buoyant recirculation region exists only at lower gas velocities. At higher gas velocities, the gas flows through the plasma instead of around it. The peak temperatures for the three cases shown are about 9200 to 9300 K and there is no significant change with flow rate. Experimental results likewise do no show any significant dependence of peak tem perature on the gas velocity. 6.4.3 Effect o f Gas Pressure The electrical properties of the plasma are a function of the discharge pressure. For this reason, discharge pressure is likely to have an im portant effect on the characteristics of the plasma. In this section, the effect of varying pressure on the tem perature contours and the flowfield is examined. In the next section, the effect of pressure on the stability boundary is considered. As mentioned earlier, local thermal equilibrium is assumed (LTE) and non-LTE effects are neglected. The results for the straight duct geometry are presented in Fig. 6.17. The inlet gas velocity is 0.1 m /s and the power is 3 kW. In Fig. 6.17a, the gas pressure is 0.5 atm. At this pressure a fairly large plasma results. The streamlines show no evidence of recirculation. At the higher pressure of 1 atm (Fig. 6.17b), the plasma appears to have shrunk in size somewhat and the blockage effect of the plasma is more pronounced in the streamline contours. At 2 atm pressure (Fig. 6.17c), these trends continue and at 3 atm (Fig. 6.17d), the plasma appears thinner. This Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.1 m/s 0 . 3 m/s 1 m/s TEMPERATURE STREAMLINES Figure 6.16 Effect of gas velocity on sphere-cylinder plasmas. P0 = 1 atm , P = 3 kW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 6.17 Effect of discharge pressure on straight duct plasmas. P = 3 kW. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 203 shape of the plasma at high pressures has been obscrved-experimentally as well. In Fig. 6.17d, the streamlines also show flow reversal taking place due to buoyancy. Increasing the pressure increases the Reynolds number which decreases the rate of diffusion; the relative importance of buoyancy over diffusion thus increases with pressure (or in other words, the ratio of Reynolds number to Froude number in creases). At large enough velocities, it is unlikely that this recirculation take place because convection would dominate buoyancy; but, at the inlet velocity of 0.1 m /s, the Froude number is small enough for buoyancy to be im portant. At low pressures, this effect is eclipsed by the dominant diffusion effects, while at higher pressures, dif fusion becomes less im portant and buoyancy dominates. The reduction of diffusion effects also explains the more compact plasma profiles at higher pressures. Similar results for the sphere-cylinder and the bluff body are shown in Figs. 6.18 and 6.19. The results are qualitatively similar to the straight duct case. Increasing the pressure reduces diffusion effects. Consequently, the plasma becomes smaller and more compact and buoyancy effects become more im portant leading to flow recirculation. Figure 6.20 shows the variation of peak plasma temperature with discharge pressure for the three geometries. The plasma temperature consistently increases with the discharge pressure. Consistent with earlier results, the plasma tem perature increases as the size of the plasma decreases. Again, the change in tem perature is not very appreciable and has not been observed in experiments. 6.4.4 Effect o f Gas Pressure on Stability Boundary At a given pressure, there is a certain minimum power below which a stable discharge cannot be maintained. This power is called the threshold power or the stability boundary of the plasma. This threshold exists because the electrons lose their energy to the heavy atoms and ions through collisions. The electrons have to Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. TEMPERATURE STREAMLINES Figure 6.18 Effect of discharge pressure on sphere-cylinder plasmas. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 205 1 atm Figure 6.19 Effect of discharge pressure on bluff body plasmas. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 206 12500 PEAK TEMPERATURE, K 10000 7500 ■000---- 2 5 0 0 ------ 0 1 D B A s t r a i g h t duct 0 —0 —© sphere-cylinder A A A b l u f f body 3 2 4- X105 PRESSURE, N /m 2 Figure 6.20 Effect of discharge pressure on peak plasma temperature for the straight duct, sphere-cylinder and bluff body geometries. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5 207 be furnished with a certain amount of energy to make up for this lost energy. When the energy is not sufficient, the electrons lose energy to the heavy particles and slow down. The recombination rate then exceeds the ionization rate and, consequently, the plasma is extinguished. Clearly, since collision processes are dependent on pressure, the threshold power is also a function of pressure. A series of calculations were performed (for the sphere-cylinder geometry) at various pressures and the corresponding threshold powers were determined. The results of this investigation are shown in Fig. 6.21. At the high pressure end, the threshold power is about 2 kW. Decreasing the pressure first reduces the threshold power and at about 1 atm pressure, the threshold power is down to about 1.2 kW. Further decrease in the pressure again increases the threshold power somewhat. Since collision frequency decreases monotonically with pressure, we might expect the threshold power to follow the same trend, that is, decrease continuously with pressure. However, this is not the case. A possible explanation for this behavior comes from the dependence of electrical conductivity on pressure. Figure 2.6 shows the variation of the ac conductivity of the plasma with pressure. At pressures high er than atmospheric, where the collision frequency is high, a decrease in pressure results in higher electrical conductivity. At sub-atmospheric pressures, a decrease in pressure results in a decrease of the electrical conductivity. Thus, the stability boundary seems to be governed by the electrical conductivity rather than the col lision frequency. A high electrical conductivity lowers the threshold power while a low electrical conductivity increases the threshold power. The question that faces us is whether these trends are in agreement with exper imental observations. Experimental measurements [27] show that with decreasing pressures, the threshold power decreases; however, at sub-atmospheric pressures, the threshold power continues to decrease monotonically with decrease in pressure. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 208 2 .5-r- THRESHOLD POWER, kW 2 . 0 ----- . 0 ----- 0 . 5 ----- 0 1 2 3 4 5 PRESSURE, atm Figure 6.21 Stability boundary for plasmas sustained in the sphere-cylinder. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 209 Thus, the experimental trend does not show the turn around at low pressures. Fur thermore, experiments show that it is possible to maintain a plasma discharge at very low pressures and, in fact, the usual method of initiating the discharge is by lowering the pressure to a few torr. At such low pressures, the electrons undergo very few collisions with heavy particles and thus possess much higher energy than the heavy particles. In other words, the electrons are not in thermal equilibrium with the atoms and ions. Maintenance of a low pressure plasma discharge is conse quently a non-equilibrium phenomemon. Since, the current physical model assumes a single temperature for the plasma, such non-equilibrium effects are not included. The equilibrium model then cannot be expected to give the right answers at low pressures. At higher pressures (above 1 atm), it is likely that the equilibrium as sumption is reasonably accurate. For microwave propulsion, the primary interest is with high pressure discharges and it is encouraging that the model qualitatively predicts the right trend at these pressures. Predicted values of threshold power range from about 1 kW to 2 kW. Exper imental measurements give values that are substantially lower—in the range of a few hundred Watts. The precise reason for this discrepancy is not clear. One pos sible source for the discrepancy is the accuracy of the property values used in the computations. The threshold power is very sensitive to transport property values such as thermal conductivity and electrical conductivity. The data used for these properties, especially the electrical conductivity, have wide error margins even for equilbrium gases at the temperatures of interest. Further, the possibility of non equilibrium prevailing in the plasma makes the accurate determination of these properties even more difficult. A related reason is the possiblity of impurities be ing present in the helium used in the experiments which may alter the transport Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 210 properties of the working fluid. In order to determine the sensitivity of the sta bility boundary to electrical conductivity and thermal conductivity, a parametric study was undertaken by arbitrarily varying these properties. Figure 6.22 shows the stability boundary of a plasma at 1 atm as the thermal conductivity was varied. The threshold power decreases with decrease in the thermal conductivity and for a thermal conductivity that is one-fourth the original value, the threshold power is down to about 300 W. Figure 6.23 shows similar results with increasing electrical conductivity (that is, greater amount of ionization). The threshold power is down to about 500 W when the electrical conductivity is increased by an order of magni tude. From the electrical conductivity data, we know that this property is a very strong exponential function of temperature. An order of magnitude difference then corresponds to a temperature change of few hundred Kelvin. Such a tem perature difference may very well be the result of thermal non-equilibrium in the plasma. Besides the accuracy of the data for the transport properties, the complex physics of the plasma may also be responsible for the discrepancy in the results. In the region of the threshold, the plasma is governed by complex kinetic rate reactions governing recombination and ionization. The plasma may be sustained as long as the ionization rate balances the recombination rate. The detailed physics of these reactions has not been included in the physical model and this might explain the inability of the model to make accurate predictions of the power threshold. To account for all the physics, it would be neccessary to incorporate a two-temperature model with finite rate “chemistry” for the ionization and recombination reactions. It should be borne in mind that the primary interest in modeling these plasmas is for propulsion applications. With regard to the feasibility of microwave propul sion, the important issues are concerning how much power may be coupled to the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 211 THRESHOLD POW ER, 1500 —— 000 500 0.0 0.2 0.4 0.6 0.8 1 .0 THERMAL CONDUCTIVITY, K / K 0 Figure 6.22 Reduction in the threshold power with decrease in the thermal conductivity of helium. k0 is the original thermal conductivity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 212 THRESHOLD POW ER, 1 5 0 0 —— 1 0 0 0 ----- 5 0 0 ------ ELECTRICAL CONDUCTIVITY, o-/ Figure 6.23 Reduction in the threshold power with increase in the electrical conductivity of helium. < t0 is the original electrical conductivity. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 213 thermal energy of the plasma, how much of this energy is transferred to the cold pro pellant gas and what thrust levels and specific impulses are available by expanding the heated gas in a nozzle. The global issues of the problem are gasdynamic issues and the reaction kinetics of the plasma do not significantly influence these issues. In particular, we are interested in the operation of the plasma under high pressure conditions where non-equilibrium is not an issue; further, we are interested in power levels that are likely to be significantly higher than the threshold limitations of the plasma device. It is noteworthy th at at higher power levels, the predictions of the equilibrium model are in excellent agreement with experimental trends and mea surements. Thus, under the conditions of interest, more detailed physical models do not seem necessary and the current equilibrium model seems adequate. However, the problem of more accurate knowledge of transport properties remains and their availability will definitely enhance the quantitative accuracy of the model. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 214 CHAPTER 7 CONCLUSIONS AND FUTURE WORK 7.1 Direct Solar Propulsion This thesis is a theoretical study of two advanced propulsion concepts—direct solar thermal propulsion and microwave propulsion. Research attention has focused on solar and microwave propulsion along with other electrothermal concepts in order to find a viable alternative to chemical rockets. The applications of these advanced concepts are for space-based applications such as space-station keeping and maneuvering, orbital transfer vehicles and interplanetary travel. The thrust requirements of these missions range from about 10-2000 N and the corresponding power levels range from low powers (1-10 kW) to high powers (100-10,000 kW). Both solar and microwave thrusters show potential for application over the range of power levels. The primary advantage of solar propulsion is that solar energy is readily available and there is no need for the development of an artificial energy source. Methods of collecting and focusing the solar energy using light-weight materials have been designed and proven to be efficient. In contrast, microwave propulsion requires an independent energy source; however, the technology for high power microwave generation already exists. Furthermore, space-stations are likely to carry low power microwave sources on board for com m u n ication purposes. Thus, microwave thrusters for station-keeping may be used without significantly adding to the weight of the space-station. However, for high power applications such as orbit transfer, microwave generators are likely to add to the weight of the engine, a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 215 factor th at has to be contended with. W ith regard to solar propulsion, both direct and indirect absorption concepts are being considered and the one that shows the most promise will be developed for practical application. Indirect solar propulsion involves solar energy absorption by a metal surface and the transfer of this energy to the propellant via a standard heat exchanger. The technology for indirect propulsion already exists and preliminary studies have indicated that specific impulses of about 600s are possible. Current efforts are directed toward fabricating a solar thruster for the Integrated Space Flight Test [14]. In contrast, the direct absorption concept, in which solar energy is directly absorbed into the propellant gas, is at a preliminary stage of development. The oretically, higher propellant temperatures and, hence, higher specific impulses are possible with direct absorption thrusters. Whether these promises can be fulfilled in a practical configuration is open to question. The physics of solar plasmas are complicated and are non-linearly coupled to the gasdynamics in the thruster. Fur ther, the technological risks with this concept are high and require a greater invest ment for development than the indirect concept. For these reasons, a fundamental theoretical analysis is required to understand the complex radiation-gasynamic in teractions and to direct the course of future experimental studies. In this thesis, direct absorption solar thrusters are studied by means of an analytical model from the viewpoint of assessing their overall feasibility for future space applications. In our study, we focus on two major issues—the efficiency of solar energy transfer to the propellant gas and the performance capability of the solar thruster. Two direct volumetric absorption concepts are considered here—(1) absorption by a suspension of small particles, and (2) absorption by alkali metal vapours. In Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 216 both cases, hydrogen is used as the carrier gas because its low molecular weight makes higher specific impulses possible. Both methods, in principle, enable higher gas temperatures than can be achieved by indirect methods, which also leads to higher specific impulses. However, there are technological issues such as the storage of particles and alkali metals and their introduction into the absorption chamber. Further, the use of particulate m atter in the thrust chamber can lead to serious erosion problems on the walls of the chamber and the nozzle. Likewise, alkali metal vapors may deposit on the inlet window and cause its degradation. These issues are critical to the successful development of the concept and might prove to be deciding factors in the final selection of the absorption scheme. In general, alkali vapors seem more promising than particles; therefore, we focus on this absorption scheme for detailed analyses of efficiency and performance. The analysis of direct solar propulsion involves the modeling of the thermal radiation transfer within the absorption chamber—that is, the absorption of the in cident solar beam and the reradiation from the hot gases to the surroundings. The physical model used in the analysis is two-dimensional in axisymmetric co-ordinates and accounts for beam convergence, radial temperature variation and radial radia tive fluxes. Previous analyses of the problem have all been one-dimensional which neglect all of these effects. For the radiation absorption, a ray-tracing technique is employed while for the radiative transfer problem, an exact integral method as well as two approximate methods—the P -l method and the thick-thin model—are em ployed. The exact integral method is used to validate the two approximate methods. A simple gray gas assumption is used to approximate particulate absorption; how ever, for alkali vapors, the radiative properties are strong functions of wavelength. Thus, a band model is used to represent the spectral dependence. The results from the analysis indicate that total absorption of the incoming Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 217 solar energy is not feasible without unduly high concentrations of alkali metal vapors and/or very large chambers. A combination of direct and indirect absorption does, however, appear promising. This hybrid absorption procedure enables a portion of the incoming radiation to be absorbed by the gas with the remainder being absorbed on a regeneratively cooled wall. In a similar fashion, the reradiated energy that is incident on the chamber walls may also be recovered regeneratively. This regenerative heat transfer from the walls not only preheats the incoming working fluid but also provides the energy for vaporizing the alkali metals without which direct absorption cannot take place. The resulting hybrid propulsive system thus provides better performance and lower metal temperatures than the purely direct absorption option. The results of three different thruster sizes, corresponding to power levels of 1, 5 and 10 MW show that high overall thermal efficiencies are possible when the conditions are optimized. About 30 to 60% of the incident solar beam can be absorbed directly into the gas but parts of this are lost due to reradiation to the walls and inlet window. Using regenerative pre-heating of the inlet gas mixture, with a proper choice of the mass flow rate, overall efficiencies of up to 90% may be achieved. In general, the overall efficiency increases with gas flow rate and optimum performance occurs for that mass flow rate that just enables the incoming flow to be heated to a temperature of 1000 K (the temperature at which the alkali metals remain in the vapor state). For all the three sizes considered, the predicted peak temperatures are around 2800 and 3200 K while the average temperatures range from 2200 to 2500 K. The solar heating results in a hot central “plasma” surrounded by cool outer gases. The strong radial temperature differences and the difference between peak and average temperatures (which are evidence of these gradients) are indicative of the need Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 218 for a two-dimensional analysis of the problem. For a chamber pressure of 3 atm, the calculations predict specific impulses of about 550 s. At this condition, the molecular weight of the hydrogen-alkali mixture is 5. Increasing the pressure to 10 atm reduces the molecular weight of the propellant to 3 while maintaining the same absorption characteristics; for the 10 atm condition, specific impulses increase to about 650 s. The specific impulse predictions for the 3 atm condition are approximately equal to those for indirect solar propulsion, while those for 10 atm are slightly better. However, it should be noted that the two-dimensional predictions for direct solar propulsion are considerably lower than previous one-dimensional estimates, and it is likely that more in-depth multi-dimensional analysis of indirect solar propulsion' will likewise reduce the estimated performance for that configuration. In general, it appears that direct solar propulsion will be competetive with or superior to indirect engines for the 1 to 10 MW size range if proper optimization of mass flux and chamber pressure is conducted. At sizes smaller than this, the indirect option will probably not only be simpler, but will provide equal or better performance. There are several limitations in the present analysis; future research should be directed toward addressing these in order to improve the reliability of the mod el. The accuracy of the present results are dependent on the radiative properties employed in the analysis. There is considerable uncertainty about the absorption cross-sections of alkali metals over some wavelength bands. Currently, there is a research effort at the Air Force Astronautics Lab [105] aimed at determining the radiative properties of alkali vapors. The results of these experiments should be more dependable than the theoretical cross-sections used here. Once these data are available, a more accurate quantification of the performance of direct solar thrusters is possible. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 219 A second limitation is the treatment of non-gray radiative transfer. The band model and the thick-thin model are both approximate models and, while their results are expected to be reasonably accurate, more detailed non-gray analyses are required to fully validate these results. An obvious choice of method is the Monte Carlo method which is best suited for treating spectral dependence. The solution procedure would then involve a full Monte Carlo solution for the radiation field for each iteration step. In general, this would be prohibitively expensive, but with the advent of massively parallel computers, efficient Monte Carlo solutions may well be within reach. Finally, the present analysis approximates the fluid dynamics in the absorp tion region to be quasi one-dimensional. This simplification enables efficient and economical computer solutions of the thermal and radiation fields in the gas. It is expected that the fluid dynamics will have only a second-order effect on the gas temperature and radiative fluxes. Hence, it is important to accurately represent the radiative properties and the non-gray radiative transfer before a full fluid dynamic analysis is incorporated. In principle, it is straightforward to extend the fluid dy namic model to a full two-dimensional analysis. A further extension is to include species equations for the alkali vapors; this would enable us to introduce the alka li metals in the central region of the chamber (where absorption is desired) and, thereby, protect the walls from the high temperature gases. These refinements of the present model will allow detailed analyses of practical thruster configurations and will aid in the actual design process. 7.2 Microwave Absorption in Gases The second problem investigated in this thesis is microwave propulsion. The objective here is to develop a physical model for microwave discharges in flowing Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 220 gases. Similar to solar absorption, microwave heating is also a possible alternative to chemical rockets for low thrust propulsion in space. Low power systems of a few kW for space station keeping as well as moderate power systems of 100 to 500 kW for orbit transfer are both of interest. Microwave propulsion shows promise for both classes of systems. Current research in microwave propulsion is mostly limited to experimental studies and very little theoretical modelling has been attempted. The present work is the first attem pt at a comprehensive two-dimensional analysis of the problem. The analysis involves the coupled solution of the Navier-Stokes equations for the gasdynamics and the Maxwell equations for the microwaves. An implicit ap proximate factorization technique is used for the Navier-Stokes solutions while a time-accurate explicit scheme is used for the Maxwell equations. In both cases, method of characteristics procedures are applied for specifying the boundary con ditions. The unsteady Maxwell solutions are averaged over a electromagnetic wave period and these average values are used for the momentum (Lorentz force terms) and the energy (ohmic heating term) sources in the Navier-Stokes equations. Lo cal thermal equilibrium is assumed for the plasma which means that the electrons and the heavy particles are at the same temperature. Accordingly, the number density of electrons in the plasma is given by the Saha equation for ionization equilibrium. Equilibrium transport properties—thermal conductivity, viscosity and electrical conductivity—are employed in the analysis. The calculations are tailored to match the conditions and geometries of the companion experiments with a resonant microwave cavity. Both free-floating and bluff-body stabilized plasmas are investigated over a range of flow velocities (0.1 to 3 m /s), pressures (0.5 to 5 atm) and incident microwave powers (1 to 4 kW). In all cases, helium was employed as the working fluid. The predictions of the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 221 model are compared with available experimental data such as the coupling effi ciency of the device and the peak plasma temperatures. There is good qualitative agreement between computations and experiments. W ith fine tuning of the cavity length, coupling efficiencies up to 99-100% are observed in both the experiments and the computations. Spectroscopic electron temperatures range from 10,000 K to 12,000 K. Computational predictions of peak plasma temperature are slightly lower—between 9000 and 10,000K. This discrepancy may in part be due to exper imental error bounds of the spectroscopic data. The underprediction may also be a result of some amount of thermal non-equilibrium between the electrons and the gaseous atoms in the plasma. Furthermore, most of the trends observed experimen tally with regard to plasma location, size and shape are qualitatively predicted in the computations as well. Based on these results, the physical model appears to provide a reasonably correct picture of the microwave discharges. The one area where calculations do not match with experimental findings is the predicted threshold power or stability boundary of the plasma. In general, ex perimental thresholds are significantly lower than the predicted values. The precise reason for this discrepancy is not known but there are several plausible causes. One possibility is the plasma model used in the present analysis. The local thermal equilibrium assumption may not be accurate over the range of pressures of interest here. A more general model would include non-equilibrium between electrons and heavy atoms and ions. Such a two-temperature model would involve the solution of additional differential equations for the electron number density and the electron temperature. Effects such as electron diffusion, electron radiation and electron-atom collisions may be accounted for in this model. An even more sophisticated plasma model would include reaction kinetics to describe ionization and recombination reac tion rates. These effects are likely to play an important role in the threshold region Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 222 of the plasma and may therefore be critical to providing an accurate prediction of the threshold power. A further limitation of the present analysis is the accuracy of the transport properties. It has been established that the stability boundary is very sensitive to the magnitudes of the thermal and electrical conductivities. Finally, it is also likely that impurities present in the experiments may be responsible for the lower observed thresholds. There are several possible avenues for future work. The best approach would probably be to pursue these concurrently. In light of the difficulties experienced with the prediction of stability boundary, a natural step would be to refine the plasma model by including a two-temperature model and, if necessary, a reaction kinetics model. A preliminary investigation of these effects would help determine how much of an impact they make on the general nature of the results. Based on this, a decision can be made on whether or not these effects be accounted for in the model. In another direction, exploratory studies on waveguide-heated plasmas and coaxial plasmatrons may also be performed. The basic model developed in this research can be easily extended to study these plasma devices as well. The best concept may then be identified for purposes of further research and development. In addition, generation of microwave discharges in other gases such as nitrogen and hydrogen are also of interest and may be studied theoretically with minimal additional work. The nature of plasmas sustained in molecular gases are likely to be different from discharges in atomic gases. Experiments of nitrogen have indicated this and, in general, it was found to be more difficult to maintain a steady discharge with nitrogen than with helium [101]. Analytical modelling of these plasmas could shed some light on the reasons for this kind of behaviour. The final objective in the investigation of microwave discharges is to study prac tical microwave propulsion systems. This would involve studying realistic thruster Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 223 geometries, including the effect of a choked nozzle as well as using higher pow er levels. This would enable quantification of thrust and specific impulses of the thruster similar to the results presented for direct solar thrusters in this thesis. The study of size scale-up is particularly important because laboratory scale experiments are often limited in the range of powers and pressures th at can be safely handled. Computational results in these regimes can complement experimental studies very effectively. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 224 REFERENCES [1] Stone, J. R. and Bennett, G. L., “The NASA Low Thrust Propulsion Pro gram,” AIAA Paper 89-2492, 25th AIAA/ASME Joint Propulsion Conference, Monterey, CA, July 1989. [2] Birkan, M., “Air Force Basic Research in Plasma Based Propulsion,” AIAA Paper 89-2495,25th AIAA/ASME Joint Propulsion Conference, Monterey, CA, July 1989. [3] Beattie, J. R. and Penn, J. P., “Electric Propulsion - A National Capabil ity,” AIAA Paper 89-2490, 25th AIAA/ASME Joint Propulsion Conference, Monterey, CA, July 1989. [4] Sovey, J. S., Hardy, T. L. and Englehart, M. A., A Bibliography o f Electrother mal Thruster Technology 1984, NASA Technical Memorandum 86998,1986. [5] Morren, W. E., Whalen, M. V. and Sovey, J. S., “Performance and Endurance Tests of a Laboratory Model Multipropellant Resistojet,” J. Propulsion, Vol. 6, No. 1, Jan-Feb 1990, pp. 18-27. [6] Tacina, R. R., “Space Station Freedom Resistojet Study,” J. Propulsion, Vol. 5, No. 6, Nov-Dee 1989, pp. 694-702. [7] Cassady, R. J., “An Analysis of Orbit Maneuvering Capabilities Using Arcjet Propulsion,” AIAA Paper 88-2832, 24th AIAA/ASME Joint Propulsion Con ference, Boston, MA, July 1989. [8] Butler, G. W., Kashiwa, B. A. and King, D. Q., “Numerical Modeling of Arcjet Performance,” AIAA Paper 90-1474, 21st Fluid Dynamics, Plasma Dynamics and Lasers Conference, Seattle, WA, June 1990. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 225 [9] Welle, R., Keefer, D. R. and Peters, C. E., “Laser Sustained Plasmas in Forced Convective Argon Flow, Part I: Experimental Studies,” AIA A Journal, Vol. 25, Aug 1987, pp. 1093-1099. [10] Jeng, S. -M., Keefer, D. R., Welle, R. and Peters, C. E., “Laser-Sustained Plasmas in Forced Convective Argon Flow, Part II: Comparison of Numerical Model with Experiment,” AIAA Journal, Vol. 25, No. 9, Sept 1987, pp. 12241230. [11] Zerkle, D. K., Schwartz, S., Mertogul, A., Chen, X., Krier, H. and Mazumder, J., “Laser-Sustained Argon Plasmas for Thermal Rocket Propulsion,” J. Propusion, Vol. 6, No. 1, Jan-Feb 1990, pp. 38-45. [12] Keefer, D. R., Sedghinasab, A., Wright, N. and Zhang, Q., “Laser Propulsion using Free Electron Lasers,” AIAA 90-2636, 21st AIAA/DGLR/JSASS Inter national Electric Propulsion Conference, Orlando, FL, 1990. [13] Perry, F. J., Solar Thermal Propulsion: An Investigation of Solar Absorption in a Working Fluid, Final Technical Report AFRPL-TR- 84-032,1984. [14] Laug, K., “The Solar Propulsion Concept is Alive and Well at the Astronautics Laboratory,” JANNAF Propulsion Conference, Cleveland, OH, 1989. [15] Etheridge, F. G., Solar Rocket System Concept Analysis, Final Technical Re port AFRPL-TR-79-79,1979. [16] Shoji, J. M., “Potential of Advanced Solar Thermal Propulsion,” Orbit Raising and Maneuvering Propulsion: Research Status and Needs, Progress in Astro nautics and Aeronautics, Vol. 89, 1989, pp. 30-47. [17] Palmer, A. J., and Dunning, G. J., “Study of Radiatively Sustained Cesium Plasmas for Solar Energy Conversion,” Final Report for NASA Contract No. NAS2-10001, Hughes Research Laboratories, CA, July 1980. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 226 [18] Rault, D. G., Radiation Energy Receiver for High Performance Energy Conver sion Cycles, Ph. D. Dissertation, University of Washington, Seattle, WA, Mar. 1983. [19] Micci, M. M., “Prospects of Microwave Heated Propulsion,” AIAA Paper 841390, 20th AIAA/SAE/ASME Joint Propulsion Conference, Cincinnati, OH, June 1984. [20] Kapitza, P. L., “Free Plasma Filament in a High Frequency Field at High Pressure,” Soviet Physics JETP, Vol. 20, No. 6, June, 1970, pp. 973-1008. [21] Meierovich, B. E., “Contribution to the Theory of an Equilibrium HighFrequency Gas Discharge,” Soviet Phys. JETP, Vol. 14, No. 5, May 1972, pp. 1006-1013. [22] Arata, Y., Miyake, S., Kobayashi, A. and Takeuchi, S., “Research of a Station ary High Power Microwave Plasma at Atmospheric Pressure,” Journal o f the Physical Society o f Japan, Vol. 40, No. 5, May 1976, pp. 1456-1461. [23] Hawley, M. C., Asmussen, J., Filpus, J. W., Whitehair, S., Hoekstra, C., Morin, T. J. and Chapman, R., “Review of Research and Development on the Mi crowave Electrothermal Thruster,” J. Propulsion, Vol. 5, No. 6, Nov-Dee 1989, pp. 703-712. [24] Whitehair, S., Asmussen, J. and Nakanishi, S., “Microwave Electrothermal Thruster Performance in Helium Gas,” J. Propulsion, Vol. 3, No. 2, Mar-Apr 1987, pp. 136-144. [25] Whitehair, S., Frasch, L. L. and Asmussen, J., “Experimental Performance of a Microwave Electrothermal Thruster with High Temperature Nozzle Materi als,” AIAA Paper 87-1016, 19th AIAA/DGLR/JSASS International Electric Propulsion Conference, Colorado Springs, CO, May 1987. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 227 [26] W hitehair, S., Asmussen, J. and Nakanishi, N., “Demonstration of a New Elec trotherm al Thruster Concept,” App. Phys. Lett., Vol. 44, No. 10, May 1984, pp. 1014-1016. [27] Balaam, P., Maul, W. and Micci, M. M., “Characteristics of Free-Floating Nitrogen and Helium Plasmas Generated in a Microwave Resonant Cavity,” IEPC Paper 88-099, 20th DGLR/AIAA/JSASS International Electric Propul sion Conference, Garmisch Partenkirchen, W. Germany, Oct 1988. [28] Balaam, P. and Micci, M. M., “Investigation of Free-Floating Nitrogen and Helium Plasmas Generated in a Microwave Resonant Cavity,” AIAA Paper 89-2380, 25th AIAA/ASME Joint Propulsion Conference, Monterey, CA, July 1989. [29] Mueller, J. and Micci, M. M., “Numerical and Experimental Investigat ions of a Propagating Microwave-Heated Plasma,” IEPC Paper 88-100, 20th DGLR/AIAA/JSASS International Electric Propulsion Conference, Garmisch Parenkirchen, W. Germany, Oct 1988. [30] Mueller, J. and Micci, M. M., “Investigation of Propagation Mechanism and Stabilization of a Microwave Heated Plasma,” AIAA Paper 89-2377, 25th AIAA/ASME Joint Propulsion Conference, Monterey, CA, July 1989. [31] Herlan, W. A. and Jassowski, D. M., “Microwave Thruster Development,” AIAA Paper 87-2123, 23rd AIAA/SAE/ASME Joint Propulsion Conference, San-Diego, CA, June 1987. [32] Kuo, K. K., Principles of Combustion, Wiley-Interscience, 1986. [33] Raizer, Yu. P., Laser Induced Discharge Phenomena, Vases, G. C. and Pietrzyk, Z. A., eds. Consultants Bureau, NY, 1977. [34] Kemp, N. H. and Lewis, P. F., Laser Heated Thruster— Interim Report, NASA CR-161665, Physical Sciences Inc., MA, 1980. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 228 [35] Kemp, N. H. and Krech, R. H., Laser Heated Thruster — Final Report, NASA CR-161666, Physical Sciences Inc., MA, 1980. [36] Gulati, A. and Merkle, C. L., “The Absorption of Electromagnetic Radiation in an Advanced Propulsion System,” Journal o f Spacecraft and Rockets, Vol. 21, Jan-Feb 1984, pp. 101-107. [37] Molvik, G. A., Choi, D. and Merkle, C. L., “A Two-Dimensional Analysis of Laser Heat Addition in a Constant Absorptivity Gas,” A IA A Journal, Vol. 23, July 1985, pp. 1053-1060. [38] Merkle, C. L., Molvik, G. A. and Shaw, J.-H., “Numerical Solution of Strong Radiation Gasdynamic Interactions in a Hydrogen Seedant Mixture,” Journal o f Propulsion and Power, Vol. 2, Sept-Oct (1986), 465-473. [39] Eddy, T. L., “Electron Temperature Determination in LTE and Non-LTE Plas mas,” J. Quant. Spectrosc. Radiat. Transfer, Vol. 33, No. 3, 1985, pp. 197-211. [40] Eddy, T. L. and Sedghinasab, A., “The Type and Extent of Non-LTE in Argon Arcs at 0.1-10 Bar,” IEEE Transactions on Plasma Science, Vol. 16, No. 4, Aug 1988, pp. 444-452. [41] Palmer, A. J., “Radiatively Sustained Cesium Plasmas for Solar Electric Con version,” Third N ASA Conference on Radiation Energy Conversion, Progress in Astronautics and Aeronautics, Vol. 61, 1978, pp. 201-210. [42] Mattick, A. T., Hertzberg, A., Decher, R. and Lau, C. V., “High Temperature Solar Photon Engines,” Journal of Energy, Vol. 3, 1979, pp. 30-39. [43] Mattick, A. T., “Absorption of Solar Radiation by Alkali Vapors,” Radiation Energy Conversion in Space, Progress in Astronautics and Aeronautics, Aca demic Press, NY, Vol. 61, pp. 159-171. [44] Rault, D. F. G., “Radiation Energy Receiver for Solar Propulsion Systems,” Journal o f Spacecraft and Rockets, Vol. 22, 1985, pp. 642-648. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [45] Thynell, S. T. and Merkle, C. L., “Analysis of Volumetric Absorption of Solar Energy and its Interaction with Convection,” Proceedings, 1988 National Heat Transfer Conference, HTD-96, Vol. 1, July 1988, pp. 269-279. [46] Venkateswaran, S., Merkle, C. L. and Thynell, S. T., “An Analysis of Direct Solar Thermal Rocket Propulsion,” AIAA Paper 90-0136, 28th Aerospace Sci ences Meeting, Reno, NV, Jan 1990. [47] Venkateswaran, S., Thynell, S. T. anti Merkle, C. L., “A Study of Thermal Radiation Transfer in a Solar Thruster,” Journal o f Heat Transfer, to be pub lished. [48] Boulos, M. I., “Flow and Temperature Fields in the Fire-Ball of an Inductively Coupled Plasma,” IEEE Transactions on Plasma Science, Vol. PS-4, No. 1, Mar 1976, pp. 28-39. [49] Wei, D., Apelian, D. and Farouk, B., “Effects of Coil Location and Injection Flow Rtae in an Inductively Coupled RF Plasma Torch,” AIAA Paper 85-1634, 18th Fluid Dynamics, Plasmadynamics and Lasers Conference, Cincinnati, OH, July 1985. [50] Rhodes, R. and Keefer, D. R., “Numerical Modelling of a Radio Frequency Plasma in Argon,” A IA A Journal, Vol. 27, No. 12, Dec 1989, pp. 1779-1784. [51] Venkateswaran, S., Merkle, C. L. and Micci, M. M., “Analytical Modelling of Microwave Absorption in a Flowing Gas,” AIAA Paper 90-1611, 21st Fluid Dynamics, Plasmadynamics and Lasers Conference, Seattle, WA, June 1990. [52] Venkateswaran, S. and Merkle, C. L., “Coupled Navier-Stokes Maxwell Solu tions for Microwave Propulsion,” 12th International Conference on Numerical Methods in Fluid Dynamics, Oxford, England, July, 1990. [53] M artin, R. E. and Keefer, D. R., “Stabilty Criterion for the Bluff-Body Stabi lized Electrodeless Arc,” Phys. Fluids, Vol. 15, No. 6, June 1972, pp. 1028-1034. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 230 [54] Balaam, P. and Micci, M. M., “The Stabilization and Spectroscopic Study of Microwave Generated Resonant Cavity Plasmas,” AIAA Paper 90-2635, 21st AIAA DGLR JSASS International Electric Propulsion Conference, Orlando, FA, July 1990. [55] Seely, S., Introduction to Electromagnetic Fields, Me Graw -Hill Book Co., Inc., 1958. [56] Johnk, C. T. A., Engineering Electromagnetic Fields and Waves, John Wiley and Sons, 1975. [57] Gandhi, Om. P., Microwave Engineering and Applications, Pergamon Press, 1981. [58] Agdur, B. and Enander, B., “Resonances of a Microwave Cavity Partially Filled with a Plasma,” J. Appl. Phys., Vol. 33, No. 2, Feb 1962, pp. 575-581. [59] Wexler,A., “Computation of Electromagnetic Fields,” IEEE Transactions on Microwave Theory and Techniques, Vol. MTT-17, No. 8, Aug 1969, pp. 416439. [60] Saad, S. M., “Review of Numerical Methods for the Analysis for the Analysis of Arbitrarily Shaped Microwave and Optical Dielectric Waveguides,” IEEE Transactions on Microwave Theory and Techniques, Vol. MTT-33, No. 10, Oct 1985, pp. 894-899. [61] Caxnbel, A. B., Plasma Physics and MagnetoSuid Mechanics, McGraw-Hill Book Co., Inc., 1963. [62] Parodi, R. and Fernandes, P., “Higher Order Modes Computation in RF Cav ities,” IEEE Transactions on Magnetics, Vol. 24, No. 1, Jan 1988, pp. 154-157. [63] Yee, K. S., “Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media,” IEEE Transactions on Antennas and Propagation, AP-14, May 1966, pp. 302-307. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 231 [64] Taflove, A., “Application of the Finite Difference Time-Domain Method to Sinusoidal Steady-State Electromagnetic Penetration Problems,” IEEE Trans actions on Electromagnetic Compatibility, Vol. EMC-22, No. 3, Aug 1980, pp. 191-202. [65] Borup, D. T., Sullivan, D. M and Gandhi, Om. P., “Comparison of the FFT Conjugate Gradient Method and the Finite Difference Time-Domain Method for the 2D Absorption Problem,” IEEE Transactions on Microwave Theory and Techniques, Vol. MTT-35, No. 4, Apr 1987, pp. 383-395. $ [66] Taflove, A. and Brodwin, M. E., “Computation of the Electromagnetic Fileds and Induced Temperatures Within a Model of the Microwave-Irradiated Human i Eye,” IEEE Transactions on Microwave Theory and Techniques, Vol. MTT-23, No. 11, Nov 1975, pp. 888-896. [67] Merkle, C. L., Computational Fluid Dynamics o f Inviscid and High Reynolds Number Flows, Dept, of Mechanical Engineering, Penn State, 1990. [68] Anderson, D. A., Tannehill, J. C. and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer, Hemispherical Publishing Company, NY, 1984. [69] Shankar, V., Hall, W. and Mohammadian, A. B. H., “A CFD-Based Finite Vol ume Procedure for Computational Electromagnetics - Interdisciplinary Appli cations of CFD Methods,” AIAA Paper 89-1987, 9th CFD Conference, Buffalo, NY, June 1989. [70] Goorjian, P. M., “Algorithm Development for Maxwell’s Equations for Com putational Electromagnetism,” AIAA Paper 90-0251, 28th Aerospace Sciences Meeting, Reno, NV, Jan 1990. [71] Jahn, R., Principles o f Electric Propulsion, McGraw-Hill Book Company, 1968. [72] Vincenti, W. G. and Kruger, C. H., Introduction to Physical Gas Dynamics, Robert E. Kreiger Publishing Company, FA, 1986. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 232 [73] Frost, L. S., and Phelps, A. V., “Transport Coefficients of Monoatomic Gases,” Physics Review, 136, 1964, pp. 1538-1544. [74] Lick, W. J. and Emmons, A. W., Transport Properties o f Helium from 200 to 12,000 K, Harvard University Press, Cambridge, 1965. [75] Siegel, R. and Howell, J. R., Thermal Radiation Transfer, second edition, McGraw Hill Book Company. [76] Niemax, K. and Pichler, G., “New Aspects in the Self-Broadening of Alkali Resonance Lines,” Journal o f Physics B, Vol. 8,1975, pp. 179-184. [77] Woerdman, J. P. and de Groot, J. J., “Emission and Absorption Spectroscopy of High Pressure Sodium Discharges,” ACS Symposium Series, Vol. 79, 1982, pp. 33-41. [78] Chetroprud, V. E., “Optical Properties of Potassium Vapours at High Temper atures and Concentrations,” High Temperature, Vol. 7, 1969, pp. 1009-1013. [79] Chetroprud, V. E., “IR Absorption Spectra of the Potassium Molecules,” High Temperature, Vol. 14, 1976, pp. 195-197. [80] Benedict, Drummond and Schlie, “Absorption Spectra of the C s 2 Molecule,” Journal o f Chemical Physics, Vol. 86, 1977, pp. 4600-4607. [81] JANNAF Thermochemical Tables, National Bureau of Standards, 1971. [82] Marr, G. V. and Creek, D. M., “Photoionization Absorption Continuafor Alkali Metal Vapours,” Proc. of the Royal Society, Vol. 304,1968, pp. 233-244. [83] Moskvin, Yu. V., “Photoionization of Atoms and Recombination of Ions in the Vapours of Alkali Metals,” 1963, pp. 316-318. [84] Zel’dovich, Ya. B. and Raizer, Yu. P., Physics of Shock Waves and High Tem perature Hydrodynamic Phenomena, Vol. 1, 1966. [85] Pomraning, G. C., The Equations o f Radiation Hydrodynamics, Pergamon Press. 1965. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 233 [86] Ozisik, N., Radiative Transfer, Wiley, NY, 1973. [87] Modest, M., Radiation Transfer in Participating Media, Dept, of Mechanical Engineering, Penn State Univeristy, 1989. [88] Thynell, S. T., “The Integral Form of the Equation of Transfer in Finite, TwoDimensional Cylindrical Media,” J. Quant. Spectrosc. Radiat. Transfer, Vol. 42, No. 2, pp. 117-136. [89] Thynell, S. T., “Treatment of Radiation Heat Transfer in Absorbing, Emitting, Scattering, Two-Dimensional Cylindrical Media,” Numerical Heat Transfer, (in press). [90] Venkateswaran, S., An Analysis of Solar Thermal Rocket Propulsion - A Pro posal for Research, Dept, of Mechanical Engineering, Penn State, 1988. [91] Warming, R. F. and Beam, R. M., “On the Construction and Application of Implicit Factored Schemes for Conservation Law,” SIAM-AMS Proceedings, Vol. 11,1978, pp. 85-129. [92] Briley, W. R. and McDonald, H., “On the Structure and Use of Linearized Block Implicit Schemes,” Journal of Computational Physics, Vol. 34,1980, pp. 54-77. [93] Mac Cormack, R. W., “Current Status of Numerical Solutions of the NavierStokes Equations,” AIAA Paper 85-0032, 23rd Aerospace Sciences Meeting, Reno, NV, Jan 1985. [94] Rai„ M. M., “Navier-Stokes Simulations of Blade Vortex Interaction Using High Order Accurate Upwind Schemes,” AIAA Paper 87-0543, 25th Aerospace Sciences Meeting, Reno, NV, Jan 1987. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [95] Briley, W. R., McDonald, H. and Shamroth, S. J., “A Low Mach Number Euler Formulation and Application to Time Iterative LBI Bcheme,” A IA A Journal, Vol. 21, Oct 1983, pp. 1467-1489. [96] Guerra, J. and Gustafsson, B., “Numerical Method of Incompressible and Com pressible Flow Problems with Smooth Solutions,” Journal o f Computational Physics, Vol. 63, Apr 1986, pp. 377-397. [97] Douglas, J. and Gunn, J.E ., “A General Formulation of Alternating Direction Method — Part I — Parabolic and Hyperbolic Problem,” Numerische Mathematik, Vol. 82, 1964, pp. 428-453. [98] Rai, M. M. and Chaussee, D. S., “New Implicit Schemes and Implicit Boundary Conditions,” AIAA Paper 82-0123, 21st Aerospace Sciences Meeting, Reno, NV, 1983. [99] Merkle, C. L. and Choi, Y.-H, “Computation of Low-Speed Flow with Heat Addition,” A IA A Journal, Vol. 25, Jun 1987, pp. 831-838. [100] Choi, Y.-H., Computation o f Low Speed Compressible Flows, PhD. Disserta tion, Dept, of Mechanical Engineering, Penn State, 1988. [101] Balaam, P., Experimental Development of a Microwave Resonant Cavity Elec trothermal Propulsion Device, Ph. D. Dissertation, Dept, of Aerospace Engi neering, Penn State, in progress. [102] Fuji, S., and Eguchi, K., “Comparison of Cold and Reacting Flows Around a Bluff Body Flame Stabilizer,” Journal o f Fluids Engineering, Vol. 103, June 1981, pp. 328-334. [103] Batal, A., Jarosz, J. and Mermet, J. M., “A Spectrometric Study of a 40 MHz Inductively Coupled Plasm a VII. Argon Continuum in the Visible Region of the Spectrum,” Spectrochimica Acta, Vol. 36B, No. 10, 1981, pp. 983-992. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 235 [104] Dymshits, B. M. and Koretskii, Y. P., “Temperature of a Free Plasma Filament in a High Frequency Field at High Pressure,” Optics and Spectroscopy, Vol. 33, No. 1, July 1972, pp. 17-18. [105] Larson, C. W., “Kinetics, Thermodynamics and Spectroscopy of Mixtures of Hydrogen and Alkali Metal Vapors at Temperatures and Pressures Approaching 3000 K and 100 atm ,” presented at 1987 Joint Meeting of the Western States and Japanese Section of the Combustion Institute, Honolulu, HA, 1987. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 236 APPEN D IX A B S O R P T IO N B A N D S O F A LK A LIS B a n d 1. 0.1 to 0.3 fim. Ground-state photoionization. fci = 0.1 x 10“ 22JVjvo + 0.12 x 10~2z N k + 0.18 x 10~22N Ct B a n d 2. 0.3 to 0.4 fim. 6p photoionization of Cs. k2 = 25 x 10~22i\TCaexp(—2.3 x 10~12/fcBT) B a n d 3. 0.4 to 0.55 fim. 5d photoionization of Cs and X —* A dimer band of Na. fc3 = 25 x 10-22JVc* exp(—2.3 x 10~u / k BT) + 183.3 x lO"22^ , B a n d 4. 0.55 to 0.65 fim. Atomic resonance line of Na. Jb4 = 0.8 x 10~46[JVjv«]2 B a n d 5. 0.65 to 0.72 fim. X —►B dimer band of Na. h = 108. x 10~22N Naj B a n d 6. 0.72 to 0.82 ^m. Atomic resonance line of K. fee = 0.25 x 10"45[Ar*-]2 B a n d 7. 0.82 to 0.95 fim. Atomic resonance line of Cs. h = 0.5 x 10-45[JVc«]2 B a n d 8. 0.95 to 1.15 fim. X —*■B dimer bands of K and Cs. k6 = 220. x 10-22Nc*2 + 250. x lO-22# * , B a n d 0. 1.15 to 1.3 fim. X C dimer band of Cs. kg = 50. x 10~22N Ctl Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. VITA Sankaran Venkateswaran was born on July 4, 1963, in Madras, India. He received the degree of Bachelor of Engineering (Hons.) from the Madurai Kamraj University, India, in December 1984. Mr. Venkateswaran has been a graduate student in the Department of Mechanical Engineering, Penn State from August 1985 to December 1990. During this period, he was also a Graduate Research Assistant under Prof. C. L. Merkle. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

1/--страниц