close

Вход

Забыли?

вход по аккаунту

?

Analysis of solar and microwave thermal propulsion systems

код для вставкиСкачать
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.
Документ
Категория
Без категории
Просмотров
0
Размер файла
9 265 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа