. .. p (grit-.0. .61 1).. ... Z!» , !.l!. C! .(l. \v ~r .1; {Iflrfittp iii-(tr . . a .. L ‘4 1 i J. . .«EA’Etfiéa rt: .. S.\...:.....I; (trip ..., ...,ri. a cré.‘l inv’lr tf.‘ .6! 4.51:. (.274 91:23: 5.. fitting... i 1 . I. .... ..l/I... t. .1. . .... E)...» {tcirvicriil r: {21.) I... r 51...? 11 .7: ......r (I 2......57 .55.: . . .1 r s: ((1! 1. . .5. 7.. 12:. (11:31: (1111.33.53... r)? 1315.111}... . 211.27.... 51.1.1! .52.. u; ‘14 v1.5; \ 5.. ...). '4! ( $317 (5..) (x . 3.. 511:: . tr;§¢-\.Qlt)l.:.x).i (liltrla; 1.11 r 1"? . . If x . .. ...: 2.5 . t. (1.} ‘49... (...: .tv‘tk? . '1‘. (Tire. I; s. v . . .1 \t.\: .1 . .(z ..r I . a an.) 5 . (lit . ... sx 1;: «i . I‘r Fr .... v... \tan iv Ar... .... {1.3.1.1 , ...).i....cntczcit5iztrLisztckzlt. . . .c ‘ . 2;...» L; r i). u\L..(c4 ~- ‘0 \/ r1]: 0.) 3 m m (D catm U C) mmml NOMENCLATURE elliptical major axis length neutral molecule ionized molecule excited molecule mth order constant arbitrary vector transition probability elliptical minor axis length ionized molecule excited molecule magnetic induction (vector) mth order constant arbitrary vector arbitrary constant arbitrary constant average number of molecular collisions mLh order constant generic heat capacity heat capacity of air heat capacity of water Comet Rendezvous Asteroid Flyby separate variable electric induction (vector) mth order constant electron charge electron electrical energy electric field (vector) energy absorbed by air energy absorbed by gas (propellant) xi [11910101 X t (D ‘< N mmmmmm CD EWSK f(t) th m order constant energy radial energy energy energy of electron state component of E lost by radiation from microwave source absorbed by microwave wall x component of E y component of E axial component of E angular component of E East-West StationKeeping time dependant function variable dependant function force flow flow flow flow flow flow flow rate of air rate of gas rate of energy into element rate of energy out of element rate of water of excited species of neutral species gravitational acceleration of earth variable dependant function degeneracy term of transition (statistical weight) Planck’s constant magnetic field radial component of H x component of H y component of H axial component of H angular component of H heat conduction / convection xii 8p INC 0 0 w 77 w cqcu ale. 7? '0 M N m 03 3J:z a I? t” x measured relative emission line intensity specific impulse incremental step size in computer program ,__i - complex number [= V -1 J current density (vector) mth order Bessel function zero order Bessel function Boltzmann constant phasor constant electromagnetic constant 32 in phasor notation Runge Kutta formula microwave probe length microwave short length mass of small object; mth order value mass of large object magnetization density (vector) mass of spacecraft mass of spacecraft plus propellant momentum loss per molecular collision average momentum loss molecular weight electron density molecular density in neutral state molecular density in excited state gas density; normalized electron density th . . . . n order Bessel function; den81ty of exc1ted spec1es National Aeronautics and Space Administration North—South StationKeeping polarization density (vector) conduction power convection power electric field power xiii (D .— exc p. ion 3 'J "I rad rec rok sel @*U"U*U"U"U"U"U'U'U'tl (D (D D FUFUFUH>HHWQ£I elastic collision power excitation power incident microwave power ionization power nth root of JO reflective microwave power radiation power electron recombination power rocket thruster power superelastic collision power total charge of object partition function momentum cross section radius within cylinder / sphere generation rate of electrons generation rate of energy radial vector component radius function of variable radius of cylinder spectral response function real component of function flow of radiation radio frequency time time function of variable air temperature electron temperature electronic temperature initial temperature water temperature transverse magnetic mode arbitrary variable scalar velocity of object xiv col dis Ndfid‘im rad tangential velocity of orbit scalar electron velocity gravitational velocity of orbit scalar velocity of propellant relative velocity of gas weight function of integral x-axis vector component y-axis vector component axial position in cylinder axial vector component axial function of variable solid angle percent energy absorbed by air coolant percent energy absorbed by plasma gas percent energy absorbed by water coolant arbitrary variable nth term variable arbitrary variable nth term variable propagation constant electric susceptibility magnetic susceptibility absolute permittivity ionization coefficient laboratory collision energy permittivity of free space relative error in numerical analysis excitation coefficient excited ionization coefficient electron recombination coefficient collisional e- recombination coefficient dissociative e- recombination coefficient radiative e- recombination coefficient K V ('D '1 m .... D 3 5 CC:~E~C>’>’XX-¢ u-Igr-MQ 05 N QQQQQ® = N VrZ n ossine d9 2.53 The average momentum loss can be found by integrating the momentum loss times the average collision over all of 6. To simplify the expression, a momentum cross section term (QM) is defined as28 n Q = 2 n J (1 - cos 6) sin 9 0 d9 2.54 0 Using this cross section term, the average momentum loss . 29 can be expressed by the following equation : = u t) N V. Q 2.55 m e r m The scattering differential cross section can be expressed considering only coulomb collisions. For this type of collision, the cross section can be calculated from the differential elementzez _ do 2 2 3 d9 8 n e m Ve sin [ The excitation rate (OVe> can be broken down into specific excitation or de-excitation terms. The following . 8 equation should hold true2 23 Experimentally, the QIn cross section term can be found from . . . . . 27 ion beam scattering uSing the follow1ng expreSSion : Q ~ 2 (a - azln en)2 = 2 Q 2.59 m 1 The Qex is the total cross section for charge exchange. The a2 term is calculated from the linearity of meQ versus the log of 6“ where En is the laboratory collision energy. A complete explanation of the theory and experimental procedure of ion beam scattering may be found in Mason and McDonald’s book27. Deionization is needed to complete the conservation of momentum. Electron recombination is the process of deionization. This recombination is highly dependant upon the electron temperature (Te). For simplicity, only approximations will be given for three primary types of recombinations (radiative, dissociative, and collisional)27. These values are approximated from experimental data of gases in various environments. 7... = Te 2.60 7 . = T'°-5 2.61 dis 9 -4.5 , 1°01 - Te 2.62 2.3.5 Conservation of Energy. Another important aspect of modeling the plasma region is developing an accurate model showing the temperature gradient within the plasma. The conservation of energy equation can be used for this model. This energy equation can be written on a differential element of the plasma at the microscopic level. The accumulation of energy (GE/6t) is equal to the flow of energy in (F1) minus the flow of 24 energy out (F0) plus the generation of energy (ren) within the differential elements. 6E , 2 This energy can be expressed in terms of power at the microscopic level. The flow terms and the generation term . . 8 can be sub-diVided by the follow1ng power terms2 F. - F = P i-I) + P 2.64 C C nd nv The radiation, conduction, and convection terms are functions of temperature, heat capacity, and boundary conditions (such as cavity wall temperature and material 9 . . type) . The electron recombination term may be neglected for low pressures. The electric field power can be written 26 as : 2 P =OE 2.66 . 28 with the conductivity coefficient written as : 2 n e e m (€+ J'w) The excitation and ionization terms are functions of the density of the species and of their excitation / ionization rates. These two terms can be written for the electron 28 as 2 25 Finally, the collision terms provide both generation and degeneration of energy. The generation type involves superelastic (or de-excitation) collisions and can be express as28 P = n n 6 <0 UL> 2,70 sei e x x d e The degeneration term involves elastic collisions of . . 28 electrons and neutral speCies which can be expressed as ° P = n n [ 2m ] <0 v > 2.71 e) e n M se [‘3 .4 Discharge Properties. 2.4.1 Energy Distribution The measurement of the energy distribution for this research involved a macrOscopic frame of reference. The macroscopic energy balance only covers energy entering and leaving the entire system. The system is defined as the microwave resonance cavity (see Figure 2.4). An energy . 9 balance written around this system can be written as E = E + E + E + E 2.72 a g a w r E8 is the energy entering the system from the microwave power source and is represented by the following expression: Ea is the energy absorbed by the air cooling and is measured by the air flow rate (F8) and its temperature rise (AT3) using the following equation (assuming the heat . 2 . 9 capaCity, LP a, is constant over the temperature range) : ) 26 Figure 2.5 Calorimetry System 27 E' is the energy absorbed by the water cooling and is measured by the water flow rate (Fw) and its temperature rise (ATV) using an equation similar to 2.74. Er is the energy escaping the system as emission radiation. Because the system is almost entirely enclosed, this radiation is neglected. Rearranging the above equations yields the following: E = P. - P - C F AT - C F AT 2.75 9 1 l' p,a a a p,w w H For this research, the above energy expressions are reported in percentages, as opposed to absolute values. The following equations indicate how these values are computed: Ea %E = x 100% 2.76 a E 8 EH %E = x 100% 2.77 w E S E 26 = 9 x 100% 2.78 g E As determined from reading (measurement) errors and statistical deviations, the following absolute error range estimates are provided: 2% for both air and water energy, 3% for power source energy, and a composite of 7% for gas energy. These percentages are given with respect to the microwave source energy. 2.4.2 Electromagnetic Resonator Modes in Microwave Cavity The experiments for this thesis use a cylindrical cavity resonator. A transverse magnetic (TMfinp) mode is 28 used. To understand the electromagnetic field in the cavity, one must analyze the three cylindrical components of the electric (E) and magnetic (H) fields. For a circular'fliflw mode, the axial (z) component of the magnetic field (Hz) is zero. A solution to Maxwell’s equations will be calculated using phasors to characterize the time dependance of the wave fields. Before solving Maxwell’s equations for the waveform, it may be clearer to define a few constants (terms). = jw (phasor constant) 2.80 Using Maxwell’s equations (number 2.27 and 2.28) and equation 2.31, one could calculate the electric and magnetic fields in the cavity. First, one must calculate the waveguide equations before considering the resonator. Using the boundary condition that Hz is zero, the radial (r) and angular (6) components of the electric and magnetic fields are calculated as follows (in phasor and cylindrical . 33 coordinates) : 6E2 r 86 - k2 E6 = -' J. (I) ‘1 Hr 2083 BE: k2 Er - 3r = - j 0.) )1 H9 2.84 1 37‘ .11 II c... 1E: m [*1 '1 N 00 01 These four equations can be simplified by separating the unknown field components and rewriting them as functions of one variable (E2). The following are the simplified equations33 kz 8E2 Er = 2 2.87 k ar C k2 6E2 E = 2088 9 k2 r 66 C j w 6 6E2 Hr = 2 2.89 k r 56 C _ j w e 8Ez H = 2090 9 k2 ar C The solution to the above equations could easily be calculated if we knew E2. Using the Laplacian operation (V2) and separation of variables (radial and axial functions) on E2, we can derive a general solution. This derivation is similar to that done in section 2.3.3. The . . . . 2 33 follow1ng is the general solution to V E2 = 0 . E = [C J (k r) + D N (k r)][A sin(m6) + B cos(m9)]exp(k z) 2.91 2 mile mmc m m 2 For finite values at the origin (r=0), one obtains a reduced solution for E233. E = B C J (k r) cos(m9) exp(k z) 2.92 In m m c 2 ['1‘] II + . Em Jm(kcr) cos(m6) exp(-j sz) 2.93 Therefore, 30 the remaining five components of this waveform are: 2 33 R E m P P P r .2 p 0 1 um run an . = - - J 6 . .’ Er p 2 [er[ R0 ] R0 “H1( R0 ]]cos(m ) exp(ijz) 2 94 mn 2 Pmn _j BPRO m E1 Jm[ R0 E = sin(m6) exp(jB z) 2.95 6 2 p R mn O ' w e R2 E J [Pmnr) -J 0 m 1 m R0 . . . H = Sin(m6) exp(jB z) 2.96 r 2 p P R run 0 . 2 -Jw EROEI In, Pmn pmn Pmn H6 = 2 P m[ R r] — R m*1[ R r] cos(m9) exp(ijz) 2.97 P 0 0 In“ H = O 2.98 2 A table showing the mth zero of Jm(Pmn) follows. waveform is characterized by this zero and is written as TM . Inn th . . Table 2.2 m Zero of Jn(Pmn) n m g 0 1 2 3 1 2.405 3.832 5.136 6.380 2 5.520 7.016 8.417 9.761 3 8.654 10.173 11.620 13.015 4 11.792 13.323 14.796 16.200 31 For a resonator, one must consider boundary conditions for E2 at the base and height of the cavity. For z=0 and z-L8 (L8 equals the height of the cavity - otherwise known as the microwave short length), the following is true: For the radius at the edge of the cavity (r=RO), the electric field reaches a maximum value. 6E Z . 8r — 0 2.100 The exponential term of E2 (see equation 2.93) may be transformed to a sine function of 2 if we take the following value for the propagation constant: The E2 term then reduces to: E = E J (k r) cos(n6) sin [w] 2.102 2 m m c L 8 Thus, the solution to this equation defines a Tanp1mode. Additionally, the frequency of the microwave must be known to calculate kc (see equations 2.79 through 2.82) . The frequency (w) is 2.45 GHz. This frequency also characterizes the electric and magnetic field components. For this experiment, TM and TM modes were used (see 011 012 Figure 2.5 for wave patterns of those modes). 2.4.3 Spectroscopic Analysis. . 39 Free electrons are an integral part of a plasma . Thus, ionization and deionization processes within a plasma occur frequently and may be thought of as the plasma’s 32 tal TM Modes Figure 2.6 Experimen 33 sustainment (or driving force). As such, the electron temperature can be seen as one of the few properties characterizing the state of a plasma. One way to measure the electron temperature is through spectroscopy. There are two distinct methods for conducting spectroscopic measurements, absorption and emission. For this research, emission spectroscopy is used. The spectrometer detects the spectral lines emitted by the gas as the atom decays from the excited state to the ground state. In the plasma, the atoms are excited from . . . . . . 38 electromagnetic radiation and from colliSions of speCies . The intensity of the spectral line is determined by two factors, the Boltzmann distribution and the transmission moment (with Einstein coefficients). The Boltzmann distribution law characterizes the population of various excited levels (Nn) using the following equation (with the degeneracy term, g“, added)38: —E . n Nn — gn N exp[ k T ] 2.103 The transition moments are determined through the interactions of the species with the emitted radiation. Using the relative emission line intensities for spectroscopic measurements, an emission line intensity going 8 through the spectrometer can be given by : 2.104 N h 8n Anm RA um. (:19 exp[ -En J 111 Q 47! RT For ground state transitions when Q becomes independent of temperatures, this equation can be reduced to provide an easier equation to calculate the electronic temperature. 34 f I x. l, E m nm n .. lnl R1 g A ‘J — constant - k T .105 n ma [‘3 For this thesis, local thermal equilibrium is assumed. Although this assumption may not provide accurate results, it does allow for a quick "ball-park" calculation of the electron temperature. This assumption allows us to assume that the electron temperature is approximately equal to the electronic (excitation) temperature. T H T 2.106 Errors are involved in this techinique to measure the electronic temperature. Although we can easily measure the relative population of an excited species, the absolute value is not accurately obtained. This comes from the procedure of measuring relative emision line. An absolute emision line can not be measured in this spectroscopy experiment. Additionally, emission spectroscopy does not provide information about the ground state species, which happens to be the most populated state. According to an analysis by Chapman, we can expect errors in excess of 50% 8 because of these limitations . A better way of measuring the electron temperature is measuring line width species such as hydrogen which undergoes a linear Stark broadening and which does not assume local thermal equilibrium. CHAPTER I I I ELECTRIC PROPULSION 3.1 Characteristics. Electrical propulsion systems use charged particles accelerated by an electric field as a working fluid. These systems are capable of creating greater specific impulses than chemical or nuclear propulsion systems. Likewise, this capability further supports their suitability for steering rockets in space. There are three basic types of electric rocket thrusters: electrothermal, electrostatic, and . 4 electromagnetic . Electrothermal thrusters use electric energy to power an arc or resistance heater to heat a conventional working fluid. Ions or colloidal particles make up the working fluid in an electrostatic system and are accelerated by an electrostatic field. Finally, a travelling magnetic and electric field system accelerates a plasma in an electromagnetic system. 3.2 Propellants. As mentioned in chapter two, specific impulse is an important parameter for propulsion systems. A higher specific impulse means lower propellant flow rate to produce a given thrust (see equations 2.1 and 2.4). In comparison, chemical systems have specific impulses in the range of 250 - 450 seconds, while electric systems have the range of 300 - 5000+ secondszo. 36 On the negative side, electric propulsion systems are characterized by low thrustzo. In other words, these systems do not possess the force needed to quickly overcome a strong gravitational counteraction. This requires longer operations to achieve the desired velocity change than chemical propulsion. For example, a mission may require several hours of operation per day for electric propulsion. An equivalent mission for chemical propulsion would require several minutes per week. Why is propellant amount so important? The propellant is a major factor in the cost of a spacecraft mission. In one example, the propellant accounted for 43% of the mass in the Galileo mission. In another example, 76% of the mass in the Comet Rendezvous Asteroid Flyby (CRAF) mission was propellant. 3.3 Applications. Currently, electrothermal propulsion performs North-South StationKeeping (NSSK) on many geosynchronous satellites’s’z’. American Telephone & Telegraph Company (AT&T) include electric propulsion options on their newest satellites - Telstar IV. This type of propulsion will be used on the Space Station Freedom as well as other man . 16,21 occupied platforms . Figure 3.1 schematically illustrates a basic NSSK satellite. The solar panels (power source) define the North-South axis. For a NSSK, the rocket thrusters act in the direction of the solar panels. This thruster configuration is more efficient than East-West StationKeeping (EWSK). Typically, two pairs of rocket thrusters are used. These pairs are located at the north and south side of the satellite. Although only one 38 thruster is needed for each side, two are used with one as a backup in case the other malfunctions. Propellant contamination is also a concern for application purposes. Problems arise when communication must be transmitted through some portion of the thruster plume. The radio frequency (RF) signal may interact with charged species of that plume. This interaction is serious as the following may occur: reflection of the signal, attenuation and phase shift of the signal, or generation of noise in the signal. 3.4 Proposed System. One type of electrothermal propulsion system uses a microwave induced plasmas. Although this system uses an electromagnetic environment, it is classified as an electrothermal because it uses a nozzle (not the electromagnetic field) to accelerate the propellant. Schematically shown in Figure 3.2 is a version of this system. The power is beamed to the spacecraft from an outside source (such as a space station or planetary base) as microwave or millimeter wave power. This power is focused onto a resonant cavity to sustain a plasma in the working fluid via heating it. The hot gas would expand through a nozzle to produce thrust. Alternately in a self-contained situation, power from solar panels or from nuclear reaction could be used to run a microwave frequency oscillator to sustain the plasma. 39 Deemed Ilium" er Millimeter 'eve Power Neale Energy “It“. lei E... ..ij Abeorpuon Chute: Propellant Storage m - Velocity Gaseou- Propellant Figure 3.2 Electrothermal Propulsion System CHAPTER IV EXPERIMENTAL SYSTEM 4.1 Introduction. The system used in this experiment was designed to conduct diagnostic measurements of three elements of plasma characteristics. At the macroscopic level, the power distribution and plasma dimensions were determined using thermocouples and visual photography respectively. And at the microscopic level, the electron temperature were measured using an optical emission spectrometer. See Figure 4.1 for the overall set-up. 4.2 Microwave Cavity. An electromagnetic system was needed to generate a plasma. The microwave cavity body was made from a 17.8 cm inner diameter brass tube. As seen in Figure 4.2, the cavity contained a sliding short and a coupling probe (the two major mechanical moving parts of the cavity). The movement of this short allowed the cavity to have a length varying from 6 to 16 cm. The coupling probe acted as an antenna which transmitted the microwave power to the cavity. The sliding short and coupling probe were adjusted (or moved) to obtain the desired resonant mode. A resonant mode represents an eigenvalue of the solution to Maxwell’s equations. Two separate resonance modes were used in these . _ 8 (Ls - 7.2 cm) and TM01 (Ls — 14.4 cm) . experiments. TM01 2 1 Additional features of this cavity included: two copper screen windows located at 90 degree angles from the coupling probe (which allowed photographic and spectral 40 41 ‘- Gee Exhauat Preeeure Gage _\ Gas Pump ' , .. - - - ) Thermocouple m I I A '-----Ia- ' ) Airlxhauet ~ m - I’ll]; 'IIIII‘ .. ...). .. aE— _ _ _ Spectroecopy a " >' i r Higowave ; 5 Ca 1’. I I y i I - i , Ineulaticn L. m - - l " " " ’ Thermocouple Water Outlet Figure 4.1 Experimental Set-up 42 LEGEND b- “fl 1. Cavity Wall 7. Sliding Short 8. Base Plate 4. Plasma Discharge Pg. 5. Viewing Window Lp. 6. Discharge Chamber Ls. Figure 4.2 Microwave Cavity Microwave Power Coupling Probe Air Cooling Chamber Gravity Force Probe Length Short Length 43 measurements), and two circular holes (in both the base and top plates) to allow propellant and cooling air flows through the cavity. 4.3 Plasma Containment. The plasma was generated in quartz tubes placed within the cavity (see Figure 4.3). The inner tube is 33 mm outer diameter and was used for the propellant flow. The outer tube was 50 mm outer diameter and was used for air cooling of the inner tube. Both tubes were about 2 1/2 feet long and was epoxied to aluminum collars. These collars fed the gas and air to and from the cavity. For additional protection, water cooling was done on the collar downstream of the cavity. 4.4 Flow System. Flow of 99.99% pure nitrogen and helium was controlled using a back pressure regulator and a 3/4 inch valve in front of the vacuum pump. A Heise gauge with a range from 1-1600 torr was used to measure the pressure of the plasma chamber. Four sets of flow meters were used to measure the gas, water, and air flows. Thermocouples were used to measure the temperature of the air and water both entering and exiting the cavity. 4.5 Microwave Power. A Micro-Now 420B1 (0—500 watt) microwave power oscillator was used to send up to 400 watts of power at a fixed frequency of 2.45 GHz to the cavity (see Figure 4.4). Although rated for 500 watts, energy was lost from the microwave cable, circulator, and bidirectional coaxial coupler. 44 Gas Outlet Water Inlet/Outlet Air Outlets Aluminum Input Collar Air Cooling Passage ) i ) _—.L;--_.l__ Thermocouple .________ Quartz Tube (50 mm 0.1).) Quartz Tube (33 mm 0.1).) Plasma Gas _)_:_ Passage 1 -————- Aluminum Output Collar .. Gas Inlet can... Air Inlet Figure 4.3 Plasma Containment Tubes 45 I—P4flla leneeeed Pains Keane hmudeat Penn PM"! I—Plulll Power lessees Adtanuetere lldl hard-.lnereunm seas IMMIedflennl Cessna Coup“. lure-Nmrlltll (e-seo ') OICHJQGOf Phanne seas Cunnnater Figure 4.4 Microwave Power Source 46 Connected to the microwave oscillator was a Ferrite 2620 circulator. This circulator provided at least 20 dB of isolation to each the incident and reflected power sensors. The circulator protected the magnetron in the oscillator from reflected signals and increased the accuracy of the power measurements. The reflected power was absorbed by the Termaline 8201 coaxial resistor. The incident and reflected powers were measured using Hewlett-Packard 8481A power sensors and 435A power meters. 4.6 Temperature Probes. Type T thermocouples (copper constantan) with braided glass insulation were placed at the inlet and outlet for the water and air cooling (see Figure 4.1). An Omega 400B Digicator was used to measure the temperature at these four locations. 4.7 Spectroscopy. The radiation emitted by the plasma was measured using a McPherson Model 216.5 Half Meter Scanning Monochromator and photomultiplier detector. A high voltage of 900 volts was provided to the photomultiplier tube (PMT) using a Harrison (Hewlett—Packard) Model 6110A (DC) power supply. The output from the PMT was processed through a Keithly Model 616 digital electrometer. The processed output is sent to a Metrabyte data acquisition & control system and recorded on a Zenith 80286 personal computer (see Figure 4.5). The monochromator was positioned about 100 cm from the plasma. The emission radiation was focussed on the monochromator using two 25 cm focal length glass lenses. This lens system concentrated the emission radiation on the entrance slit opening of the monochromator. To optimize the 47 O O Zenith — Metrabyte xeithly cl 0 80200 —- Ins-13 Digital Ilectrcmeter liicrewave Cavity . Run' (Ell McPherson no.8 ‘ .-. ..-"...- n liensehrometer Focusing Lenses B-P 01104 (DC) Power Supply Figure 4.5 spectroscopy System 48 intensity of the spectroscopic emissions, the slit widths for this experiment were set at 100 microns for the entrance slit and 50 microns for the exit slit. The atomic spectra was taken using the 1200 grooves per mm grating (plate) with a range of 1050 - 10000 A. This groove setting allowed for a large range of wavelengths to be observed. The reciprocal linear dispersion was 16.6 A per mm. The focal length of the spectrometer was one half meter. CHAPTER V ENERGY DISTRIBUTION 5.1 Pressure Dependence. The energy distribution within the microwave cavity was analyzed over various pressures, ranging from 200 to 800 torr. Using the same microwave system, Hoekstra had conducted this experiment using the two separate modes (TM011 and TM012’ and the two pure gases (nitrogen and helium) . In each of these experiments, the power used was about 250 watts with air cooling flow of 2 SCFM and water cooling flow of 5.75 ml/sec. Hoekstra’s experiments involved gas flow in the direction of the gravitational force (F9). An additional experiment was done reversing the flow of the gas. Using this new data and Hoekstra’s results for the same experiment, an estimate is done to determine the gravitational effects on the other three Hoekstra pressure experiments. Reverse of the gas flow showed a slight change for energy absorbed by the gas (see Figure 5.1). Nevertheless, a significant difference in the power distribution was observed. About 5% of the total energy, which was absorbed by the cavity wall, was redistributed to the air cooling. Thorough cleaning of the cavity before the new experiment may be the main cause of this difference. Dirt and oil deposits on the cavity walls would absorb a large amount of input energy. Thus, removing this dirt and oil would account for the large energy distribution change from the wall. This new energy distribution agrees with the work done by Chapmane. Figures 5.2 - 5.4 shows estimated energy 49 50 0:03 Snowmen. CON? COOP com com 00¢ CON 0 _ _ _ c t F t L 9:32.35 28.3w I t . o Aobmxoczv .2 III - 35.9.35 «00 xlx Acpczncnoxv 2363 ole 10F Accesses—31v a? mlm Acctanctoxv mow *WVR 1 ION 4 Ion . n n ” Toe m\m\m\1)m_ tom $88 NEE .. 626$ mczmmotd mnmco> 220545ng ..oBoa paqiosqv Jemod % Calorimetry Graph I Figure 5.1 5‘) Ath Rammed OONF 000? com 000 00¢ CON 0 _ _ u L e t L - b 1—( n - Acbmxoozv .863 o Acnymxoozv .__< 4 639.005 moo AouoEmmmv .303 1 o P @6823 .2 . AouoEmmmv wow 0) Q 0 )C 1. ON .\\I.I|III r I on T 0+. T T on 1 00 r 2. A6668 :os: .. 826$ mtzmmotd msmto> cossngma toe/0d peqiosqv JGMOd % Figure 5.2 Calorimetry Graph II 52 AtBV otammota 0mm com Om; Accumxooxv 2303 639.005 .__< 659.85 «00 AcaoEfimmv .303 5305883 :4. AcuoEwmuv moo CON _ LP IIIIII, R V/Ise I 388 NEE .. ceasezv mtzmmoca msmco> cocbctpmé toe/om o no? Tom pquosqv JGMOd % Figure 5.3 Calorimetry Graph III 55 €65 0czmm0tn. com com T t u . 6.39.005 .303 ele 6.59.005 ...< Ts 6.0.9.005 moo xlx 66:50.5 800.03 ole 66E$mmv ._.< 810 \M 66:505 000 I 00¢ CON _ u h 0 0 T 4 / I A6668 :05 .... eoootezv oesmmotd mam..0> c0::£.:.m.n_ ..0>>on. 10F WON peqiosqv JGMOd % Figure 5.4 Calorimetry Graph IV 54 distribution values for the reversed gas flows compared to Hoekstra. 5.2 Flow Dependence. The power absorption by the helium gas was calculated over a range of various flow rates at a constant pressure of 200 torr using an air cooling rate of 2 SCFM and water cooling rate of 5.75 ml/sec. This experiment was only conducted in the TM012 mode. Once again using Hoekstra’s results, an estimate is made for each mode (see Figure 5.5). Hoekstra conducted a similar experiment in the TM011 mode. With the exception of the power distribution (which may be caused by dirt on the cavity wall), the gravitational force had an insignificant effect on the plasma. The difference between the data for flows with and against gravity were within the range of experimental error. A difference caused by gravity may be observed by significantly reducing this error. Unfortunately, that difference can not be observed using the current experimental system. 0.3 Plasma Power Absorption. The power initially absorbed by the plasma is the combination of the power absorbed by the gas and by the cooling air. For helium, power absorbed by the plasma was about 67% in the TM012 mode and about 80% in the TMO11 mode. The remaining power was that power which is absorbed by the cavity walls through radiation and convection. This plasma power absorption was the initial anticipated power we can expect to use for propulsion. As an example, Figure 5.6 shows the total power absorbed by the 55 300$ 23. so: _ comm com F 00.0 F :o E 6363 86.8: I :o 3 $88.82 62326: To NE E 668308 85.00: I «.0 E. 6303 283on In com. o _ . o )b 805 Kq peqiosqv JeMOd % A26; omm .. 83:65 06W. 30.... m:m..0> comatomnd. tosod 000 Figure 5.5 Calorimetry Graph V 56 A2008 2.6m 2,2... com P 000 F com com 004 r p p L _ a _ b b 0mm h VA 6668 N55 .. 82.65 b o om Tom Two row owsold Kq peqiosqv JQMOd % 08.0w. so... mzmt0> cochtomn< t0>>0d Figure 5.6 Calorimetry Graph VI 07 plasma for various flow rates, which was dramatically different than the power absorbed by the gas as shown in Figure 5.5. The power absorbed by the wall can be recovered by such methods as using heat exchangers. A heat exchanger can transfer energy to the propellant before it enters the cavity. CHAPTER VI PLASMA DIMENSIONS 6.1 Introduction. An important element of microwave generated plasmas is its size and shape. Its size varies about pressure, flow, and power changes. These measurements can be used to model heat transfer and chemical processes within a plasma. Plasma dimensions were calculated from 35 mm color slide pictures. The projected images were measured and compared with a known grid pattern. This grid pattern is photographed and projected at the same distance as the plasma. Water cooling flow rate of 5.75 ml/sec and air cooling flow rates of 2 SCFM for helium gas and 3 SCFM for nitrogen gas were used. A 35 mm camera, on a tripod, was position about 2 cm from the microwave cavity window (see Figure 4.2). Pictures were taken using a 50 mm lens, using Ektachrome 200 ISO color slide film. Several aperture and shutter speed settings were analyzed to determine the optimal settings. The following data were reported for aperture and shutter speed settings of f2.8 and 1/2503, respectively. Using only the TM012 mode, four sperate measurements were made for each plasma: height and width for strong and weak regions. Plasmas for both helium and nitrogen were observed to have two distinct regions - a more intense lighter color inner region surrounded by a diffuse darker outer region. Thus, it was appropriate to record two separate measurements, identifying the inner and outer regions as strong and weak ionization regions respectively. 58 59 6.2 Pressure and Flow Dependency. Measurements were done on helium plasmas in the pressure range of 400 - 1000 torr using a power source of about 250 watts. Three gas flow rates (0, 572, and 1144 SCCM) were used for each pressure point (see Figures 6.1 - 6.3). Nitrogen was observed in the pressure region of 400 - 500 torr using the same power. Three gas flow rates (0, 102, 204 SCCM) were used for each pressure observation (see Figures 6.45 - 6.6). It was noticed that plasma size decreases with increases in pressure and increases slightly with increases in gas flow rate. 6.3 Power Dependency. Measurements were done on helium plasmas at a constant pressure of 400 torr with no gas flow for both nitrogen and helium. The power was varied between 100 - 250 watts for helium (see Figures 6.7 - 6.8) and between 210 - 260 watts for nitrogen (see Figure 6.9). The primary observation was that small changes in power at these pressure ranges affected the plasma dimensions. Increases in power resulted in a larger plasma, as expected. 6.4 Plasma Color and Shape. For both helium and nitrogen plasmas, the strong ionization regions displayed an intense white color. The color difference is observed in the weak ionization region. Helium displayed a purple color, while nitrogen displayed an orange color. Additionally, these two gases displayed a different shaped plasma. Both shapes can simply be represented by an oblate ellipsoid. The helium plasma differed slightly in the middle with a small indentation, resulting in a "dumbbell shaped" figure. Figure 6.10 shows 60 0:05 0.300021 00NP 000v 000 P IF lb) P L 20.0.. c060. 0:030 ele .30.; 5.00.. 0:030 III 50.0.. c060. 0.003 Ole .323 c0608 v.00; film 000 00¢ 00m h h m))m./.w/m [of .0606 NE E .. so... 060 2000 8 ..00 To... 3 -04 (mm) qibue'] $1.3mean m3m..®> mCO_mC®E_D UEmUF. E3:®I Plasma Dimensions Graph I Figure 6.1 61 CON F FIJHW1 «:90; .362 953% £23 .368 9.8% £90: :29: x33 £23 :29: x02: Ath 0.5881 000— com a . _ 3111 L omm $11.71/- Tim/m/m / Om: h com 385 NS 5: ... 20m 80 208 3.8 mazmmmi m3m$> chmcmEE oEmoE E26: 0 T9. fl .3 (ww) qibuej Figure 6.2 Plasma Dimensions Graph II 62 com _‘ W 0:80 whammoli 000— com 000 00¢ 00m 0 r H? . _ . F . _ . a . £20: :29! 0:95» I 0 £23 c208 ocean I a. 290: 539. x003 ole Fm £23 569. x003 E s T0? a 1.: I Jul JI/l .. Tm/m/m 10m ImN Ton rm», 10¢ 10¢ AoUoE NE 5:. l 26.... moo 200m 9315 oSmmoLa mamao> mcoficmcEQ oEmoE E320: (ww) Lnfiue‘] Figure 6.5 Plasma Dimensions Graph III 63 0:05 ounmmoi 0mm 6m 2.: 00m Emma... :20? 9.9% I 0 £23 c069. macho I r 2905 :20? x83 Ole rm £23 .339. xooz mlm - T0— r 1m: 1. f as U 70m pm. I/I I U. ) m. m ImN w 1 w ( 100 100 .ll/ f 10+ T 0 Lo 7. Se $88 NS 5: ... 25: 80 zoom 8 oSmmoi m3m$> mcoficoch oEmoE comobaz Figure 6.4 Plasma Dimensions Graph IV 64‘ foot oSmmoi 0mm JIM 0mm 0m; 00m. 29oz c209. 9.8% I 0 £23 c038 mcobm Ill T 290; .362 x00; olo In £23 c802 3003 Elm 0. nor a Tm? f 10m .ll/I ,. rmm mTIIIlml/lm f In a QYI\IIIIA'III////-. IAwmw r0¢ Q\\\1c o a row 385 NE E ... 22... 80 zoom 80 (mm) qibuej oSmmoi m3m$> chmcmEE oEmoE comobLZ Figure 6.5 Plasma Dimensions Graph V 65 0mm C83 Samoan. r: 0mm 0a.; 00m. «:30... c209. 9.5.3» I 0 £23 .368 9.93m I f 290: :29: x83 @Io Tn £23 :20»: x00? film. 0. mo? 10 F f TON T Wmm E} U 1m 1 10m. a 100 0}]. 1 10+ I\® a rm¢ 388 NS 5: ... 26; 80 zoom :5 (mm) qibue'l ousmmoi m3m$> chmcmEE oEmoE comobLZ Figure 6.6 Plasma Dimensions Graph VI 66 9.93 Loaoa 00m .1 00m com on? 00? on b ,. Li + u b h b b _ k 3%.: 0 £25 4 onoE NS 3:. I COZONEB 95me 0 0 1 am 1 10F (ww) qibue‘] Egon. msmuo> mco_mcoE5 oEmoE 83:0: Figure 6.7 Plasma Dimensions Graph VII 67 9520 season 00m om F 00 P on P h b b - com omm r r L L p u r 22»: o £25 < GEE NS 5: I 8:322 x85 Loioa msmuo> chmcoEE oEmoE E26: rom (ww) qifiuej Figure 6.8 Plasma Dimensions Graph VIII 68 90.95 Egon. 0mm mwm 0mm 5505.3. 93.3% I «:90: ole cozofico. 93.5w I 523 «In cozofico. x8? I «:20: I cozofico. x33 I 523 I 0mm 4\4\I\L 1\«\4\\4 % ESE NE s: .. so; 80 020 com r00 (ww) qibue'l 630a mzmuo> mco_mcoE_0 oEmoE 808:2 Plasma Dimensions Graph IX Figure 6.9 69 Figure 6.10 Photograph of Nitrogen Plasma in TM012 Mode 3. 70 a photograph of a nitrogen plasma in a TM012 mode (courtesy of Maris Mantenieks of NASA - Lewis Research Center). 0.5 Plasma Volume. The plasma volume was easily obtained assuming axial symmetry and neglecting the indentation of the dumbbell shape. The volume was calculated as though the plasma region behaved as an oblate ellipsoid. The simple equation for this calculation is: . .. 2 Volume = g“ X {_ES%EEE_] x [-El%£2-J 6.1 The data from Figures 6.1 - 6.6 were calculated for this volume determination. This calculated volume data is provided in Table 6.1. As seen from this data, the volume was dependant upon both pressure and flow rate. When pressure decreased or flow rate increased, the volume of the plasma increased. 6.6 Mechanical Observations. A small mechanical measuring device was constructed by Hoekstra and modified by Haraburda to measure objects from a distance. As shown in Figures 6.11 and 6.12, an 8 x 7 inch object constructed of two steel plates with a fixed hole on one side and an adjustable iris on the other was a device designed to provide immediate feedback information concerning plasma dimensions. Using this device, one could obtain a reasonable measurement within minutes. On the other hand, a more accurate measurement could be obtained through photography within hours. Thus, this device could be used during an experiment to determine any irregularities in the dimension. This was helpful because the experimental 71 Table 6.1 Helium and Nitrogen Plasma Volumes Cavity Strong Weak EXPERIMENT presses s... s Resign J Itorr) (Cm ) (Cm I 47+ 400 2.50 11.68 HELIUM 600 4.83 9.75 a . - 800 4.63 8.24 ‘0 SLCM FLOW) 1000 4.70 8.10 400 6.35 13.15 HELIUM 600 5.10 10.29 _ a“ l. 800 4.85 8.53 (.512 SLLM FLuw) 1000 4.77 8.01 400 6.76 13.18 HELIUM 600 5.55 10.69 , -- 800 5.05 8.45 (1144 SCCM FLUW) 1000 4.77 8.37 NITROGEN 400 10.13 15.23 - a - 450 9.10 15.54 ‘” SLCM FLUW’ 500 8.84 14.70 NITROGEN 400 11.74 16.39 450 9.91 15.83 (102 500” FLOW) 500 10.14 16.34 NITROGEN 400 11.87 16.51 A q - 450 10.79 16.74 (204 SLCM FLO“) 500 10.52 16.38 72 Figure 6.11 Mechanical Measuring Device - Side View 73 ~~-—-— Iris Adjuster - Iris Meter Figure 6.12 Mechanical Measuring Device - Rear View 74 run was normally terminated before the photographic measurement. To use this device, the iris side is placed against the microwave viewing window. The person looks through the other hole at the plasma and adjusts the iris to contain the plasma, both its height and width. Using the scale provided, a plasma dimension can be obtained by multiplying the reading value by a calibrated number. This device has a measurement error of 10%. CHAPTER VII SPECTROSCOPY 7.1 Introduction. The electron temperature was an important parameter characterizing the plasma. The spectroscopy experiment was used to calculate (or approximate a calculation of) that temperature. Chapter two described the theory behind this calculation and chapter four described the experimental system. Before doing spectroscopic experiments, two things were done. First, a computer program was written (in Basic language) to use the data acquisition system. The adjustment coefficients in the program and the frequency (timing) of data acquisition were Optimized such that the measurement errors were less than 0.5% from the digital electrometer. Second, the spectrometer was calibrated using a Model 2450 tungsten lamp (from Optronics Laboratories, Inc.) with known intensity readings. The known intensity readings were obtained from the 17 AUG 82 calibration of the lamp. A current of 6.500 t 0.001 amps was supplied to the lamp. Instrumental readings were taken for the spectrometer over a wavelength range of 3000 3 through 9000 A. The results of this calibration are provided in Table 7.1. Linear interpolation was used to calculate the spectral response function (Rx) for specific wavelengths. The Rx’s in Table 75 76 7.1 were calculated using the following non-dimensionalized equations: Measured Intensity - Known Intensity _ R - . X . 1.1 1 Known Inten51ty x Measured Inten31ty =6300 _ Known Intensity - RX - 120 [ Measured Intensity ] ('2 To use equation 2.105, several coefficients must be provided. Because the electron temperature was calculated from the slope of the equation calculations, several wavelengths must be measured. From several scans of the spectrum, four strong and recognizeable transition regions (wavelengths) were observed. The data for those wavelengths are provided in table 7.21. Although equation 2.105 requires those transitions to be ground state originating, not enough strong transition regions were observed to calculate an electron temperature. Therefore, equation 2.105 was used for non-ground state transitions. Table 7.2 Electronic Transition Values for Helium xnm RX gn Anm En 0.0320 3.851E-18 0.0948 3.687E-18 0.2460 3.803E-18 2945.11 0.087 3888.65 0.086 4471.48 0.220 1 8361.69 25.51 (DUICOCO 7 7 Table 7.1 Spectrophotometer Calibration (using 6.5 amp Tungsten Lamp) Known Measured Known Measured 1 Intensity Intensity X X Intensity Intensity X 3000 0.00358 3.06 0.141 6100 0.4692 63.33 0.890 3100 0.00547 2.77 0.238 6200 0.4914 62.64 0.942 3200 0.00735 2.84 0.311 6300 0.5136 61.66 1.000 3300 0.0109 4.14 0.316 6400 0.5358 49.67 1.295 3400 0.0145 9.86 0.177 6500 0.5580 34.26 1.956 3500 0.0180 21.04 0.103 6600 0.5796 21.78 3.195 .3600 0.0234 35.97 0.078 6700 0.6012 13.70 5.268 3700 0.0288 52.37 0.066 6800 0.6228 9.40 7.958 3800 0.0374 61.66 0.073 6900 0.6444 6.48 11.93 3900 0.0461 62.64 0.088 7000 0.666 4.89 16.34 4000 0.0547 63.13 0.104 7100 0.683 3.89 21.08 4100 0.0680 63.62 0.128 7200 0.700 3.55 23.69 4200 0.0812 63.87 0.153 7300 0.717 3.35 25.68 4300 0.0945 63.87 0.178 7400 0.734 3.45 25.54 4400 0.1077 63.87 0.202 7500 0.751 3.47 25.95 4500 0.1210 64.11 0.227 7600 0.763 3.45 26.55 4600 0.1348 64.36 0.251 7700 0.775 3.74 24.85 4700 0.1586 64.11 0.297 7800 0.788 4.06 23.29 4800 0.1774 64.36 0.331 7900 0.7998 4.21 22.82 4900 0.1962 64.36 0.366 8000 0.812 4.04 24.15 5000 0.2150 64.36 0.401 8100 0.819 3.82 25.76 5100 0.2372 64.36 0.443 8200 0.825 3.87 25.62 5200 0.2594 64.36 0.484 8300 0.832 3.82 26.17 5300 0.2816 64.36 0.525 8400 0.839 4.01 25.10 5400 0.3038 64.11 0.569 8500 0.846 4.04 25.16 5500 0.3260 64.11 0.611 8600 0.852 4.43 23.10 5600 0.3502 64.11 0.656 8700 0.859 4.70 21.95 5700 0.3744 63.87 0.704 8800 0.866 5.14 20.23 5800 0.3986 63.87 0.749 8900 0.872 5.80 18.05 5900 0.4228 63.62 0.798 9000 0.879 6.41 16.46 6000 0.4470 63.38 0.847 78 7.2 Experimental Results. An experiment was done using helium in the TM012 mode with no gas flow at a power of 220 watts. The pressure range was 400 - 800 torr. For this experiment, the electron temperature was assumed to be that of the electronic temperature under the assumption that local thermal equilibrium occurs (see equation 2.106). The results are provided in Figure 7.1. Because there was a shape and magnitude difference between my data and that of Hoekstra, I plotted the same data for a similar experiment from Mueller and Micciaz. Although not shown, Chapman’s results were important in that his results showed the electron temperature dependence upon pressures. Additionally, it shows results using Hoekstra’s technique of data acquisition and interpretation. The volume electron temperature is expected to decrease as pressure increases. It takes less energy to maintain a smaller volume of plasma (which decreases with increasing pressure). However, the peak electron temperature is expected to increase with pressure increases. Chapman’s results for very low pressure (0.5 - 10 torr) showed a . . . . . 8 decrea51ng temperature With an increa51ng pressure . The magnitude difference can be explained because a different technique in data acquisition occurred. Hoekstra used the same technique Chapman did by measuring the peak heights on a chart recorder. Chapman’s results had temperatures around 5,0000K. My data were calculated measuring the area under the peak using a data acquisition system with a computer. 79 000 P owm Trot Samoan. ooo 00¢ r _ CON 1? 1- o 039.00: min 85 a. 3:3: 7.. 09:50.61 olo 1(\m\m E] ? l/o 4:] 36; 0mm .. 26: oz .. 82.2.: maze/EMQEE zoEomd (>1 000 L) eJnioJedwei Electron Temperature of Helium Figure 7.1 80 Therefore, the electron temperature is expected to be near 13,0000K for pressures near atmospheric. Additionally, the electron temperature is expected to decrease with increasing pressure. CHAPTER VIII COMPUTER MODEL OF ELECTRON DIFFUSION 8.1 Introduction. An explanation of the diffusion of the electron in the plasma region was briefly covered in chapter two. An analytical solution was obtained by neglecting the generation term in the continuity equation. As previously mentioned, this is typically assumed for very low pressures. Unfortunately, this assumption is not valid for high pressures, such as atmospheric. As pressure increases, the density of the gas increases. This increased density results in an increase in the collision rates, such as recombination and ionization. The purpose of this chapter is to develop a numerical method to see the effects that electron recombination and electron ionization have upon the distribution gradient of the electron density within the plasma. 8.2 Development of Mathematical Model. An assumption is made that the ion and electron densities are equal. As a result, equation 2.34 (continuity equation for electrons) can be re-written as (with N defined as the normalized electron density): 81 82 This equation includes ambipolar diffusion. For no generation, this equation can be simplified by separating N into separate variables as mentioned in chapter. However, this method will not work for generation because of the non-linear recombination term on the right hand side. This thesis only considered the radial part of this equation. To calculate the electron density along the radial axis about the center plane (z=0) can be done by saying that: This is valid because the change of N is symmetric about the origin. The resulting equation is a second order differential equation: 2 - —d—1:— : —_ dN + (I N2 - VIN 8'3 dr r dr This equation was re-written as a set of two first order differential equations by setting: dN dr Therefore, the following set of first order differential equations were used to determine the electron density distribution in the radial direction along the center plane: dN _ - dD = - D + a N2 - U N 8.6 dr r I These equations were solved using the normalized density. The normalized density was the non-dimensional quantity with 83 N equal to one when the electron density is at its maximum value (or the center of the plasma). For this simulation, the maximum density will be assumed to be 1x1012 electrons / cma. This density uas used because Chapman observed electron densities around this values. Using this assumption, the recombination coefficient must be normalized. This coefficient is a function of temperature. For the assumed maximum electron density and for a temperature range of 250 - 64,0000K, the recombination coefficient has a range of 1x10.13 - 1x10"6 . . 3 40 recombinations cm / sec . The two boundary conditions required to solve this problem are: 8.3 Numerical Analysis. A variable-step Runge Kutta method will be used to solve these two differential equations with the given boundary conditions. A fourth and fifth order evaluation with six evaluations per step will be done. These two evaluations will be done to estimate the relative error of each step. The following formulas for . . - . . 34 solV1ng equation 8.0 Will be used in the computer program dN. - _ l f(ri,Ni) - dr 8.9 84 K2 = f(r}+- ‘* , N +--fr— K1) 8.11 K3 = f(ri+ % INC, N + INC[§% K1+‘§% K2] ) 8 12 K4 = f(ri+ %§ INC, N + INC[:?3$ K1- gig? K2+ :13: K3) ) 8.13 K5 = f(rl+ INC, N + INC[%%% K - 8 K2+ 3g§g K3 4:3: K4) ) 8.14 N1+1 = N + INC(§%§ K1 + :22: K3 + :13: K4 - % K5) 8.16 x:::"= N + INC[T%§ K1+ 13:3? K3+ :2335 K4- 3% K5+ 3% K6] 8.17 The above formulas are the save for solving equation 8.6. Because those two equations (8.5 and 8.6) are coupled, they must be solved simultaneously. The relative error per unit step can be written as: 15:::t) ' Ni+1| S €rel TEL§QET 'Nil 8'18 The variable step in this program is the recalculation of the step size if the relative error is exceeded for that step. The gamma value (Ire!) is the fraction of the old increment (INC) step size predicted to satisfy the error limitation. Because this method is a fourth order estimation, the fraction required will be the fourth root of the error equation (8.18) and will be defined as: INC (Nil 0.25 e 1 7 7 (2;) 8'19 IN - N. I (b - a) J 1+1 1+1 85 The new step size will be defined as: INC("9"’ = 0.8 7 INC(°ld) rel As needed, the program will recalculate the step size until the size is small enough to satisfy the error limitation. The program for this numerical analysis is located in appendix B. The boundary conditions for this differential equation are not both initial value conditions. Because they contain an initial and final value, a "shooting method" approach was used for the numerical method. An estimation for the initial derivative of N at R=1 is predicted. This prediction is continually changed until the final boundary condition is obtained. 8.4 Computer Simulation. Simulations were done using the above method (see appendix B). One simulation looked at the effects of recombination while the other looked at the effects of ionization. The step sizes for the simulations were reduced until the step size required for the no generation simulation resulted in a solution matching the analytical solution of equation 2.49. In the first simulation, the ionization frequency was set to zero. Four runs of the program were done varying the dimensionless recombination coefficient from 0 to 1000 (or from 0 to 1x10_9 recombinations cm3 / sec in dimensional form). For low recombination coefficients, the plasma electrons extend to the wall. As seen in Figure 8.1, higher values of recombination cause the plasma to contract. As 86 Samoa 8831 04 md one .50 Nd 0.0 3.600 cozocBEooom mEbo> £25 ommzoctozv ZOEmOd ._<_Q :52le Zomeomjm Kigsueg 110110913 Radial Electron Density Gradient - Figure 8.1 Recombination Effects 87 observed in chapter six, the plasma width becomes smaller as pressure increases. This simulation illustrates why that happens. Additionally, the difference between strong and weak regions can be implied by the shape of those curves (in Figure 8.1). The curve clearly shows that the center of the plasma has the highest density of electrons. The outer region becomes less dense. Depending upon the cutoff density for each observed region, an electron density profile may be generated for plasma conditions observed in chapter six. The second simulation held the dimensionless recombination at 1000 and varied the ionization frequency from 0 to 30 ionizations / sec. Because of the large amount of energy required to ionize a neutral atom, multiple collisions are required by the electrons for ionization to occur. Therefore, the ionization is expected to be small and almost negligible. As seen in Figure 8.2, the electron density near the wall becomes slightly larger when ionization occurs. These simulations illustrate that an accurate model describing the electron distribution must include the generation term. These simulations do not provide a very accurate model as they used constant values for recombination and ionization. These values are a function of temperature and pressure, which in turn depends upon position. Nevertheless, the simulations do provide some insight into what the electron gradient is and why it is for various situations. 88 Cosmo; 888m 0. F 0.0 . owo .30 N0 0.0 IN.0 4.0 10.0 0 4.0.0 0.. ON 0n 1 0; 6883021.. cozoNEB mEbo> ct; ooN__oE._ozv ZOEmOd 4,304.1 m3mmm> :5sz ZOwFOmjm Kilsueg uonoelg Figure 8.2 Radial Electron Density Gradient - Ionization Effects CHAPTER IX CONCLUSIONS Several different conclusions can be drawn upon the work of this thesis. Potential use of this type of thruster is important. Also, the need for conducting a microscopic theoretical analysis for plasma phenomena is addressed. The three diagnostic techniques discussed provide needed information to characterize the macroscopic energy transfer within a microwave-induced plasma. Finally, computer numerical models provide useful information concerning plasma transport phenomena. Because the microwave-induced plasma electrothermal rocket system experimentally displays similar characteristics to other electrothermal rocket systems, it has shown its potential application for spacecraft propulsion. Electrothermal thrusters lack high thrust, but have high specific impulse which is very useful for platform station keeping. Investigating the microscopic theory of plasma particles is very important. From the conservation of particles and momentum, the shape of the plasma can be predicted. From the conservation of energy, the gradients of important parameters can be predicted, such as temperature and pressure gradients. Therefore, modeling transport properties must include a microscopic theoretical analysis. 89 90 The energy distribution experiment provided us with the importance behind the condition of the containment wall. A dirty wall can absorb a large amount of power. Reversing the gas flow showed a negligible difference in the energy distribution. The power distribution appears to stabilize for flow rate beyond 500 SCCM. Finally, initial power absorption to the plasma is 67% in the TM012 mode and 80% in the TM011 mode. The plasma dimensions experiment provided us with the effects that pressure, flow rate, and power have upon the volume and shape of the plasma. The volume is decreased as pressure increases, flow rate decreases, or power decreases. Additionally, high flow rates tend to elastically elongate the shape of the plasma. Unlike nitrogen plasmas exhibiting an ellipsoidal shape, helium plasmas exhibit a "dumbbell" shape. The spectroscopy experiment provided us with an estimate for the electron temperature and the effects pressure had upon that temperature. The electron temperature is approximately 13,0000K in the atmospheric pressure region. Finally, the electron temperature decreases as pressure increases. The computer model provided us with the importance of the recombination and ionization effects. For low pressures when recombination and ionization coefficients are negligible, the wall recombination plays a very important role in the shape of the plasma. The plasma reaches the wall. For high pressures when those coefficients are important, wall recombination becomes less important. The plasma forms an ellipsoidal shape and does not reach the wall. CHAPTER X RECOWIENDATIONS Improvements can be made on future research in the areas of equipment modification, experiments, and theoretical modeling. 10.1 Equipment Modifications. An accurate model for the plasma region must account for the flow pattern (i.e. the velocity profile). A major area which needs to be addressed is flow in and out of the plasma. For example, does the flow by-pass or pass through the plasma discharge? If both occur, how much of the flow by-passes the discharge? To answer these questions, modification to the collars of the plasma tubes should be done. The collars should be modified such that the flow velocity contains an angular component (swirling flow). Furthermore, this modification should include an adjuster such that the axial helix interval can be varied. Additionally, a modification to the collars should include a way for an additional gas to enter. This additional gas can be used to conduct streamline and pulse change experiments. Another important problem with modeling plasmas in a microwave resonance cavity deals with power absorption to the walls. For an efficient thruster system, power absorbed by the walls must be minimized. As seen in chapter 5, the condition of the cavity wall influences power absorption to the wall. Two main types of energy are 91 92 absorbed by the wall — radiation and convection. To minimize radiation energy to the wall, the inside wall of the cavity should act as a mirror. A mirror coating, such as silver, should reflect most of the radiation energy to the plasma. When the investigations of analyzing the plasma region and the thruster region result in an accurate model, one should consider combining the two experiments. For example, a nozzle should be connected to the cavity system to investigate thrust and efficiency measurements. This should result in optimizing the location of the plasma relative to the nozzle. The quartz tube provides a large unknown temperature gradient. Presently, we can estimate the outer wall temperature of the inner tube as being the temperature of the exiting air coolant. An infrared temperature probe can be used to measure an intermediate level temperature of the plasma. This temperature can be approximated to be that of the inner wall temperature. The wall effects on a plasma are important. One such effect is wall recombination. This wall recombination affects the size and shape of the plasma. To investigate this effect, one should conduct experiments using various sizes and shapes for the confinement tube. Finally, with part of the experimental system being automated, one should consider automating the entire system. The next step in automating the system includes connecting both the temperature probes from the calorimetry experiments and the power meters (incident and reflective) from the microwave power source to the data acquistion 93 system. Once this is completed, a computer program should be written to numerically analyze the data as it is being collected. 10.2 Proposed Experiments. Several new areas of experiments should be considered. Some of them involve the equipment modifications as mentioned above. Dependent upon those modifications, the following experiments can be conducted. Flow pattern experiments can be done to determine the flow to and from the plasma as a function of the environmental variables (i.e. pressure and flow conditions). Power absorption to the walls of the mirror coated cavity should be compared to the mirrorless cavity. Simultaneous thruster and plasma experiments should be conducted. Estimates of the inner wall temperatures should be compared to changes in the environment. Finally, simultaneous experiments should be done using a completely automated system to provide more accurate data for constant environmental conditions. Experiments which can presently be done involve power increases and gas mixtures. Increasing the power to above one kilowatt should be done. This would provide a more realistic model characterizing the dependence upon power. Finally, gas mixtures should be done starting with non-reactive binary gases. This work should continue into investigating tertiary reactive gases. This evolution of experiments may someday branch off into a new type of rocket thruster mention by Pollardls. This new type of thruster is a hybrid chemical / electric thruster. The reactive gases could simulate a chemical thruster with high thrust, while the non-reactive gases could simulate an electric thruster 94 with high impulse. This is a new concept which should seriously be considered. 10.3 Theoretical Modeling. Given the large amount of data collected, work should be done in correlating the analyzed variables. For example, we would need a model to estimate the power distribution and plasma dimensions for a given environmental condition (i.e. pressure and flow pattern). An initial step towards modeling the plasma region was taken in chapter eight. The weakness of that model stems form its unrealistic assumptions (constant temperature and pressure throughout the plasma discharge). A more accurate model should be developed which accounts for the non-zero temperature gradient. Additionally, because the electron is not the only important species in the plasma, the gradient profile of the other species (i.e. ions and neutrals) should be determined. Another important model involves energy and species transfer to and from the plasma region. A model should be developed to estimate the energy and mass transport properties between the plasma and its adjoining fluid. 10.4 Recommendation List. The following summarizes and lists the above recommendations: 1. To modify the plasma tube collars to allow experiments to be conducted for investigating flow patterns near the plasma discharge. 95 2. To mirror coat the resonance cavity for minimizing power absorption to the wall. 3. To combine the cavity and thruster experimental systems for analyzing the optimal location of the plasma with respect to the nozzle. 4. To use an infrared temperature probe for estimating the inner wall temperature. 5. To use different size containment tubes for investigating the wall effects upon the plasma. 6. To completely automate the experimental system so that one can conduct simultaneous experiments. 7. To conduct increased power experiments for providing data with power above one kilowatt. 8. To conduct mixture experiments for providing an expanded data base beyond pure gas experiments. 9. To correlate the given and future data for predicting power distribution and plasma dimensions. 10. To develop a theoretical model for determining the species and temperature gradients within the plasma discharge. 11. To develop a theoretical model for predicting the transport properties between the plasma and the adjacient fluid. [\3 REFERENCES CRC Press, Inc., Handbook 9; Che_istry and Physics, 68th ed. [1988). National Research Council Panel on the Physics of Plasmas and Fluids, Plasmas and Fluids, Washington D.C., National Academy Press [1986]. Hellund, E.J., The Plasma State , New York, Reinhold Publishing Corp. [1961]. Dryden, H.L., "Power and Propulsion for the Exploration of Space," Advances in Space Research , New York Permagon Press [1964]. Langton, N.H., ed., Rocket Propulsion, New York, American Elsevier Publishing Co. [1970]. Hawley, M.C., Asmussen, J., Filpus, J.W., Whitehair, S., Hoekstra, C., Morin, T.J., Chapman, R., "A Review of Research and Development on the Microwave-Plasma Electrothermal Rocket", Journal 9f Propulsion and Power, Vol 5, No 6 [1989]. Moisson, M. and Zakrzewski, Z., "Plasmas Sustained by Surface Waves at Microwave and RF Frequencies: Experimental Investigation and Applications", Radiative Processes in Discharge Plasmas, Pletnum Press [1987]. Chapman, R., "Energy Distribution and Transfer in Flowing Hydrogen Microwave Plasmas.‘ Ph.D. Dissertation, Michigan State University [1986]. 96 10. 11. 12. 13. 14. 16. 17. 97 Bennett, C., and Myers, J., Momentum, Hgaty and Mass Transfer, 2nd ed., McGraw-Hill, Inc., New York [1974]. Hoekstra, C.F., "Investigations of Energy Transport Properties in High Pressure Microwave Plasmas." M.S. Thesis, Michigan State University [1988]. Q Eddy, T.L., "Low Pressure Plasma Diagnostic Methods,' AIAA/ASME/SAE/ASEE 25th Joint Propulsion Conference [1989]. Eddy, T.L., and Sedghinasab, A., "The Type and Extent of Non-LTE in Argon Arcs at 0.1 - 10 Bar," IEEE Transactions 93 Plasma Science, Vol 16, No 4 [1988]. Cho, K.Y., and Eddy, T.L., "Collisional—Radiative Modeling with Multi-Temperature Thermodynamic Models," Journal Quant. Spectrosc. Radiat. T an fer, Vol 41, No 4 [1989]. Eddy, T.L., "Electron Temperature Determination in LTE and Non-LTE Plasmas," Journal Quant. Spectgosc. Radiat. Trangfer , Vol 33, No 3 [1985]. Pollard, J.E., and Cohen, R.B., Hybrid Electric Chemical Pro ulsion, Report SD-TR-89-24, Air Force Systems Command [1989]. Stone, J.R., Recent Advanceg in Low Thrust Propulsion Technology, NASA Technical Memorandum 100959 [1988]. Sovey, J.S., Zana, L.M., and Knowles, S.C., Electromagnetic Emiggion Experienceg Using Electric Propulgion Systems ; A Surve , NASA Technical Memorandum 100120 [1987]. 20. 21. 22. 24 26 98 Hawkins, C.E., and Nakanishi, 5., Free Radical Propulsion Concept, NASA Technical Memorandum 81770 [1981]. Aston, G., and Brophy, J.R., ”A Detailed Model of Electrothermal Propulsion Systems," AIAA/ASME/SAE/ASEE 25th Joint Propulsion Conference [1989]. Beattie, J.R., and Penn, J.P., "Electric Propulsion — National Capability," AIAA/ASME/SAE/ASEE 25th Joint Propulsion Conference [1989]. Stone, J.R., and Bennet, G.L., The NASA Low Thrust Propulsion Program, NASA Technical Memorandum 102065 [1989]. Morin, T.J., "Collision Induced Heating of a Weakly H Ionized Dilute Gas in Steady Flow, Ph.D. Dissertation, Michigan State University [1985]. r, M.B., "Life Support Systems," Military Posture ; 1985, Joint Chiefs of Staff [1985]. 0 r11 «1:» w Whitehair, S.J., ”Experimental Development of a 7 Microwave Electrothermal Thruster,’ Ph.D. Dissertation, Michigan State University [1986]. Halliday, D., and Resnick, R., Physics, New York, John Wiley and Sons [1978]. Samaras, D.G., Theory 9f Ion Flow D namics, Englewood Cliffs, N.J., Prentice—Hall, Inc. [1962]. Mason, E.A., and McDonald, E.W., Transport Properties 9f Ions lg Gases, New York, John Wiley and Sons [1988]. 28. 29. 30. 31. 32. 33. 34 36. 99 Cherrington, B.E., Gaseoug Electronics and Gas ngggg, New York, Pergamon Press [1979]. Cambel, A.B., Plasma Physics and Magnetofluid Mechanics, New York, McGraw Hill Book Company, Inc. [1963]. Panofsky, W.K.H., and Phillips, M., Classical and Magnetism, 2nd ed., Reading, Massachusetts, Addison-Wesly Publishing Co. [1962]. Electricity Johnson, L.W., and Ross,R.D., Numerical A al sis, 2nd ed., Phillipines, Addison-Wesley Publishing Co. [1982]. Mueller, J., and Micci, M., "Investigation of Propagation Mechanism and Stabilization of a Microwave Heated Plasma," AIAA/ASME/SAE/ASEE 25‘h Joint Propulsion Conference [1989]. Atwater, H.A., Introduction t9 Microwave Theory, New York, McGraw Hill Book Company [1962]. Goodger, E.M., Principleg 9f Spaceflight Propulgion, Oxford, Pergamon Press [1970]. Jahn, R.G., Physics g: Electric Pr0pulsion, New York, McGraw Hill Book Company [1968]. Davis, H.F., and Snider, A.D., Vector Analysis, Dubuque, Iowa, Wm. C. Brown Publishers [1988]. Kreyszig, E., Advanced Engineering Mathematics, New York, John Wiley & Sons [1988]. 38. 39. 40. 41. 100 Bromberg, J.P., Physical Chemistr , Boston, Allyn and Bacon [1980]. Nicholson, D.R., Introduction Lg Plasma Theor , New York, John Wiley & Sons [1983]. McDaniel, E.W., Collision Phenomena lg Ionized Gases, New York, John Wiley & Sons [1964]. Haraburda, S., and Hawley, M., "Investigations of Microwave Plasmas (Applications in Electrothermal . t. h . Thruster Systems)," AIAA/ASME/SAE/ASEE 25 J01nt Propulsion Conference [1989]. APPENDIX A MATHEMATICS A.1 Introduction. The following appendix highlights the vector operations, phasor transformations (from time domain), and series solutions / orthoganalility used in this thesis. Only cartesian and cylindrical coordinates will be discussed. It is noted that spherical coordinates can be used, but are not used in this paper. . 36 A.2 Vector Operations The following expressions and equations are used throughout the presentation of this thesis: Cartesian Coordinate Vector: — A A A E = E x + Eyy + Ezz A.1 X Cylindrical Coordinate Vector: E = E f + E60 + E Q A.2 I' 2 Vector Addition: E + H = (E + H )Q + (E + H )9 + (E + H )2 A.3. x x y y 2 z 101 102 Scalar Vector Multiplication: -' A CE = 0Ex+aEy+oEzz Vector Dot Product: E - H = E H + E H + E H xx yy 22 Vector Cross Product: — A E x H = (E H -E H )2 + (E H -E H )9 + (E H -E H )z yzzy zxxz xyyx Gradient of a function: Vf(X)YyZ) = ax X 1' 8y y + 7Z— _ af A 6f 6f A ”We,“ - W1" + W8 + oz 2 Divergence of a Vector: _ 6Ex 6Ey 6E V'E(XIY)Z) = a—x 7' a—y + 62 _ 6(rEr) 8E9 aE V'E(r,0,z) = -;—§;f—— + r 89 + 62 Curl of a Vector: _ atz aEy)A aEx 6122]A aEy aEX A vxE(x,y,z) = 5;— - 5;_Jx + — Jy + 5§_ - z VXE e - BEE _ 359 A , IE1 _ 6E2 @ + l filgfii _ SEE A (r’ ’z) ' rae 82 r 82 or r ar 88 Z 103 Laplacian of a function: 2 2 2 a ‘ a V2f(x,y,z) = V-(Vf) = —:+ :4. g 4.13 0x 6y 62 2 2 6 6 a vzf‘r'e’z’ = r6r[r 8:) 285+ : 4-14 r 89 az A.3 Phasor Transformation. Phasors are defined as transformation of a differential or integral of a time dependant function into an algebraic expression. To accomplish this, one must rewrite the time domain function as the real component of a complex phase domain function. Let: f(t) = cos(wt) A.15 Rewritten as real component of complex phase domain function: f(t) = Re[cos(wt) + j sin(wt)] A.16 Differential Differential equation: d f(t) _ d cos(wt) d [j sin(wt)] A ——dt—__ - Re[ dt + dt 4.16 —QE%L£1 = Re[-w sin(wt) + jw cos(wt)] A.17 d f” = -o sinwt) 4.18 dt 104 Define Differential Phasor: _§t_ = 30 A.19 Substitution: d—d31 = 3(1) f(t) A.20 —QE%L£1- = jw Re[cos(wt) + j sin(wt)] A.21 -ga§i£l- = Re[jw cos(wt) - w sin(wt)] A.22 %(t) - = J1 f(x) g(x) dx A.36 Inner Product of Orthoganal Functions: = 0 for f(x) ¢ g(x) 4.3T Rewriting y(x): y(x) = C1f(x) + 023m) + . o . A.38 To solve for C1, one must multiply each side by f(x) and a 10’ weight function, w(x), and integrate over the domain of the function. fw(x)f(x)y(x)dx = Iclw(x)f(x)f(x)dx + [C2w(x)f(x)g(x)dx + °°° A.39 By orthoganality of the series: Jw(x) f(x) y(x) = fC1w(x) f(x) f(x) dx + 0 A.4O Rearranging the terms: I w(x) f(x) y(x) dx I W(x) f(x) f(x) dx - . . 37 Some useful integrals of Bessel functions Ix“ J (x) dx = x“ J (x) + C A.42 n-l n 1 [J (x) dx = [J (x) dx - 2J (x) A.43 . n+1 n-1 n -n —n Ix J (x) dx = -x J (x) + C A.44 n+1 n 1 R 2 R2 2 [x J (x) dx = -—- J (R) A.45 n 2 n+1 o . 39 A.5 Useful Vector Properties The following are several useful vector relationships which become very helpful when solving vector related equations such as Maxwell’s equations for electromagnetics. ) = (X x E) - E = B - ( x X) A.46 OI OI I - (8 x TX(EXC>=(X-E)B-(X B)C 447 (AXE) (EXB)=(7§ E) (E E)-(X 5) (E C) 448 V(fg)=ng+gi A.49 V-(fX)=IV-I+X-Vf 1.50 VX(fX)=vaX+foX A.51 VZX=V (V-K) -Vx (WK) A.52 Vfo=v-(VxX)=o A.53 APPENDIX B COMPUTER PROGRAM The following is the computer program used for calculating the normalized electron density for chapter eight. PROGRAM ARMYl ************************************************************** * x * PROGRAM: Runge Kutta Solution for Electron Diffusion * * NAME: Scott S. Haraburda * * DATE: 1 January 1990 * *DESCRIPTION: This program is designed to estimate electron * * concentration gradient within the plasma region * using the fifth order Runge Kutta numerical method* using variable step size control. This control is* ******************************* VARIABLES:HSTART HMIN HMAX RTOL RERR Q R RQ RR NQ NR DT ROLD RNEW NOLD NNEW GAMMA IFLAG RKl-6 NK1-6 EST A1-5 R1-6 I , IR RECOMB REST accomplished using an absolute error test. The The The The The The The The The The The The The The The The The The The The The The The The The The The Initial Step Size for the Program. Minimum Step Size for the Program. Maximum Step Size for the Program. Relative Tolerance Level. Relative Error for a Given Step. Initial N Value in the Interval. Latter N Value in the Interval. Value of R at Point Q. Value of R at Point R. Value of N at Point Q. Value of N at Point R. Given Step Size for the Program Old Value of R. New Value of R. Old Value of N. New Value of N. Calculated Fractional Step Size Termination and Write Flag. Runge Kutta Functions. Runge Kutta Functions. Estimated Error. Runge Kutta Coefficients for New Runge Kutta Coefficients for EST. Loop Integer for Main Program. Collision Rate Frequency Value. Recombination Rate Frequency Value Estimated Value for R. 109 ***************************** 110 * NEST — The Estimated Value for N. * * NT - The intermediate New Value of N. * * RT — The Intermediate New Value of R. * * * ************************************************************** INTEGER I, IFLAG REAL HSTART,HMIN,HMAX,RTOL,RERR,Q,R,DT,RECOMB,IR,REST, + NERR,NEST DOUBLE PRECISION RQ, RR, NQ, NR OPEN (9, FILE = ’ARMYl', STATUS = ’NEW’) 1 PRINT *, ’Input the Recombination Coefficient [RECOMB]:’ READ *, RECOMB PRINT *, ’Input the Collision Frequency [IR]:’ READ *, IR PRINT *, ’Input the Initial Derivative Value [Nle’ READ *, NQ PRINT *, ’Do You Want the Results Saved [1-yes, O-no]?’ READ *, IFLAG HSTART = 0.01 HMIN = lE-6 RTOL = 1E-6 HMAX 0.01 R = 1. RQ = 0. DT = 0.01 PRINT 75 PRINT 100, R, NQ, RQ, R, R, HSTART Q = 1. IF (IFLAG .EQ. 1) WRITE (UNIT = 9, FMT = *) Q, RQ DO 5 I = 1,99 R = l. - (I - 1) * DT Q = 1. - I * DT CAL DESOLV (Q,R,NQ,NR,RQ,RR,RTOL,HSTART,HMIN,HMAX, + IFLAG,NERR,RERR,RECOMB,IR) IF (IFLAG .EQ. 3) GOTO 10 PRINT 100, Q, NR, RR, NERR, RERR, HSTART IF (IFLAG .EQ. 1) WRITE (UNIT = 9, EMT = *) Q, RR NQ = NR RQ = RR 5 CONTINUE IF (IFLAG .EQ. 0) GOTO 1 GOTO 50 10 PRINT *,’Algorithm failed. Last point was [R,NOLD,ROLD]:’ PRINT *, Q, NR, RR CONTINUE FORMAT (//,T5,’r’,T14,’N(r)’,T28,’R(r)’,T39,’N Rel Err’, + T53,’R Rel Err’,T66,’HSTART’,/,73(1H-)) 100 FORMAT (T2,F5.2,T9,E12.5,T23,E12.5,T37,E12.5,T51,E12.5, + T65,F7.4) CLOSE (9) END ‘1 0| Cl 0 111 +SUBROUTINE RKF (NOLD, NNEW, ROLD, RNEW, NEST, REST, INC, RECOMB, IR, R) ************************x************************************* * * * The Runge Kutta fourth / fifth order method to estimate * * the next value of N(r) & R(r) and their estimated error. * * * ************************************************************** INTEGER I, IFLAG REAL NEST,REST,INC,NK1,NK2,NK3,NK4,NK5,NK6,Q,R,DT,IR, + RK1,RKZ,RK3,RK4,RK5,RK6,RECOMB,RTOL,RERR,NERR,RN DOUBLE PRECISION ROLD,RNEW,NOLD,NNEW,NT,RT *******************x****************************************** s * * Define the Runge Kutta coefficients. * * * ************************************************************** B21 = 1. / 4. B22 = 3. / 8. B23 = 12. / 13. B31 = 3. / 32. B32 = 9. / 32. B41 = 1932. / 2197. B42 = -7200. / 2197. B43 = 7296. / 2197. 851 = 439. / 216. B52 = —8. B53 = 3680. / 513. BS4 = -845. / 4104. B61 = —8. / 27. B62 = 2. 863 = -3544 / 2565. B64 = 1859. / 4104. B65 = -11. / 40. A1 = 25. / 216. A3 = 1408. / 2565. A4 = 2197. / 4104. A5 = —1. / 5. R1 = 1. / 360. R3 = -128. / 4275. R4 = -2197. / 75240. R5 = 1. / 50. R6 = 2. / 55 NT = NOLD RT = ROLD 112 ************************************************************* * * * Begin the Runge Kutta numerical analysis. * * *********************#****x********************************** RN = R NKl = -NT/RN + RECOMB*RT*RT - IR*RT RKl = NT NT = NOLD + INC * B21 * NKl RT = ROLD + INC * 821 * RKl RN = R + INC * B21 NK2 = -NT/RN + RECOMB*RT*RT - IR*RT RK2 = NT NT = NOLD + INC * (B31*NK1 + B32*NK2) RT = ROLD + INC * (B31*RK1 + B32*RK2) RN = R + INC * 822 NK3 = -NT/RN + RECOMB*RT*RT — IR*RT RK3 = NT NT = NOLD + INC * (B41*NK1 + B42*NK2 + B43*NK3) RT = ROLD + INC * (B41*RK1 + B42*RK2 + B43*RK3) RN = R + INC * 823 NK4 = -NT/RN + RECOMB*RT*RT — IR*RT RK4 = NT NT = NOLD + INC * (B51*NK1+B52*NK2+B53*NK3+B54*NK4) RT = ROLD + INC * (B51*RK1+B52*RK2+B53*RK3+B54*RK4) RN = R + INC NKS = -NT/RN + RECOMB*RT*RT — IR*RT RK5 = NT NT = NOLD+INC*(BB1*NK1+B62*NK2+B63*NK3+B64*NK4+B65*NK5) RT = ROLD+INC*(861*RK1+B62*RK2+BS3*RK3+864*RK4+B65*RK5) NNEW = NOLD + INC * (A1*NK1 + A3*NK3 + A4*NK4 + A5*NK5) RNEW ROLD + INC * (A1*RK1 + A3*RK3 + A4*RK4 + A5*RK5) RN = R + INC / 2. NK6 = -NT/RN + RECOMB*RT*RT — IR*RT RK6 = NT NEST = R1*NK1 + R3*NK3 + R4*NK4 + R5*NK5 + R6*NK6 REST = R1*RK1 + R3*RK3 + R4*RK4 + R5*RK5 + R6*RK6 RETURN END 113 SUBROUTINE GAMCAL (Q,R,NOLD,ROLD,NEST,REST,RTOL,GAMMA) ************************************************************* * * * * * * * * The method to evaluate the necessary step size for the * given tolerance levels and precalculated estimated error. * The mixed control is evaluated here using both the relat. * and absolute tolerances. * * * *********************************************************** DOUBLE PRECISION NOLD,ROLD REAL T1,T2,ABSEST,NEST,REST,Q,R,RTOL IF (NEST .NE. 0.) THEN T1 = ABS (NOLD / NEST) ELSE T1 = 1E30 ENDIF IF (REST .NE. 0.) THEN T2 = ABS (ROLD / REST) ELSE T2 = 1E3O ENDIF ABSEST = MIN (T1,T2) ************************************************************* * * Calculate the gamma value for the variable step algorithm.* * * 1: ************************************************************* IF (ABSEST .EQ. O.) GOTO 1 GAMMA = (RTOL * ABSEST / (R - Q)) ** 0.25 RETURN T1 = MAX (NEST,REST,RTOL/10) GAMMA = (RTOL * RTOL / (T1 * (R — Q))) ** 0.25 RETURN END 114 SUBROUTINE DESOLV (Q.R,NQ,NR:RQ:RR,RT0L1HSTART.HMIN. + HMAX,IFLAG,NERR,RERR,RECOMB,IR) ************************************************************ * * * The main subroutine combining the Runge Kutta algorithm * * along with the variable step size calculation for a given* * interval [Q,R]. Minimum and Maximum step sizes are neede* * This subroutine is terminated if the calculated step size* * is less than the minimum step size provided. * * * ************************************************************ INTEGER IFLAG REAL Q,R,HSTART,HMIN,HMAX,RTOL,RERR,NERR,NEXT,REST, + HOLD,DT,T,NEST,RECOMB,IR DOUBLE PRECISION NQ,NR,RQ,RR,ROLD,RNEW,NNEW,NOLD HOLD = HSTART T = R NOLD = NQ ROLD = RQ ************************************************************ * * * Using given parameters, the new value and required step * * size is calculated. * * * ************************************************************ 1 CALL RKF (NOLD,NNEW,ROLD,RNEW,NEST,REST,HOLD,RECOMB,IR,T) CALL GAMCAL (Q,T,NOLD,ROLD,NEST,REST,RTOL,GAMMA) HNEW = (0.8) * GAMMA * HOLD IF (GAMMA .GE. 1.) GOTO 2 IF (HNEW .LT. HOLD/10.) HNEW = HOLD / 10. IF (HNEW .LT.HMIN) GOTO 3 HOLD = HNEW GOTO 1 *********************************************************** a x * Calculated new value for R is acceptable, and next * * parameters are established. * * * *********************************************************** 2 IF (HNEW .LT. HMIN) GOTO 3 IF (HNEW .GT. 5.*HOLD) HNEW = 5. * HOLD IF (HNEW .GT. HMAX) HNEW = HMAX IF (T-HOLD .LE. Q) GOTO 4 T = T - HOLD HOLD = HNEW NOLD = NNEW ROLD = RNEW l GOTO 115 *********************************************************** * * * Algorithm failure flag and last values are set. * * * *********************************************************** 3 IFLAC = 3 R = T NR = NOLD RR = ROLD RETURN *********************************************************** * * * Computation of final step in interval, along with the * * relative and absolute errors (of that final step). * * * ****************************************************#****** 4 CONTINUE HSTART = HNEW HOLD = T-Q CALL RKF (NOLD,NNEW,ROLD,RNEW,NEST,REST,HOLD,RECOMB,IR,T) NR = NNEW RR = RNEW IF (RR .88. 0.) THEN RERR = 1. ELSE RERR = ABS(REST / RR) ENDIF IF (NR .EQ. 0.) THEN NERR = 1. ELSE NERR = ABS(NEST / NR) ENDIF - RETURN END TINIll/()TTTIII