3 1293 01 93 0547 LIBRARY Michigan State University ‘mmimwmnmnmmm This is to certify that the dissertation entitled THE ENERGETICS OF HYDROGEN ATOM RECOMBINATION : presented by John Waldemar Filpus has been accepted towards fulfillment of the requirements for Ph.D. degree in Chemical Engineering “777% 6% Major professor } Date &! ”-1425"?- ucrn...n:a__..- - - m‘ ”7, , . . .1 0.127" MSU LIBRARIES “ RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES will be charged if 500E is returned after the date stamped below. THE ENERGETICS OF HYDROGEN ATOM RECOMBINATION by John Waldemar Filpus A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1986 ABSTRACT THE ENERGETICS OF HYDROGEN ATOM RECOMBINATION by John Waldemar Filpus The dissociation-recombination cycle of hydrogen is likely to be important in the operation of a proposed microwave-plasma electrothermal rocket for deep-space propulsion and maneuvering. A series of prelimi- nary experiments on a small-scale microwave-plasma electrothermal rocket system and two modelling studies, one involving vibrational and rota- tional molecular excitation and one on more macroscopic phenomena downstream from the plasma, are described. Plasmas at pressures less than 100 Torr were studied in hydrogen, nitrogen, and argon, at flow rates up to 2 mg-mol per second, and inci- dent microwave power between 300 and 600 W. Thrusts generated by the system were measured, by a vane-type thrust stand, and temperature profiles below the plasma were estimated, by a gas-dynamic technique, involving pressure ratios in flow through a nozzle. The measured thrusts were proportional to the square root of the gas-dynamic temperatures. The temperatures dropped rapidly from 1200 K in the plasma to 300 to 400 K just downstream, depending on the direction of flow of the coolant. Cooling of the plasma tube dominates the temperature profiles obtained. Directions for further improvement of the experimental arrangement are suggested, in the realm of improved matching of the John Waldemar Filpus input power and the capacity of the gas, moving the nozzle closer to the plasma, and increased plasma pressure. The influence of vibrational and rotational energy of the recom- bined molecules on the rate of recombination and release of the energy of recombination is studied. The relaxation times for vibrational ex- citation are on the order of hundreds of times shorter than the resi- dence time for the recombination reaction over a range of conditions that cover the range of practical designs. Molecular internal energy is found to be negligible for a practical electrothermal rocket system. The development of computer models of the recombination reaction zone, downstream of the plasma, is described, and the models are compared with the experimental results. The models confirm the suggested improvements for the experimental arrangement. The results indicate that a two-dimen- sional, axially-dispersed model, with a catalytic surface, is necessary to model all the characteristics of the experimental results. The axial dispersion model is evaluated by a search technique on the initial derivatives. Initial work provides qualitative confirmation, though quantitative agreement is lacking, and directions for improvement of the model are indicated, by combining the two-dimensional model and the dispersion model, and by fitting model parameters and initial conditions to the results. To my Parents iv ACKNOWLEDGMENTS This work would not have been possible without the assistance of the people at NASA-Lewis Research Center, Cleveland, Ohio, especially Shigeo Nakanishi, Aerospace Engineer, and John Mallinkey, Technician. The text was written using the facilities provided by the Engineering Computer Facility, College of Engineering, MSU, and most of the figures were created using the PLOTIT program on ECF equipment. Thanks must go to Martin C. Hawley, for his steady encouragement and critiquing of the final product. TABLE OF CONTENTS List of Tables List of Figures Nomenclature 1. Introduction 2. Technical Background 3. Equilibrium Considerations 4. Experiments Apparatus Procedure Gas-Dynamic Temperature Theory Experimental Program Specific Impulse - Measured and Theory Thrust Measurement Study Discussion 5. Vibrational Modelling 6. Computer Modelling Adiabatic, Gas-Phase, Plug Flow Model Cooled Model Surface Reaction Two-Dimensional Model Axial Dispersion Model 7. Conclusions vi Page viii ix xii 14 17 24 24 35 37 40 57 62 68 82 83 87 93 114 121 133 vii Appendices . Computer Program Descriptions Variable Dictionary . Computer Program Listings Program 1 - Vibrational Model Program 2 - Reactive Surface Program 3 - Two-Dimensional Model Program 4 - Axial Dispersion Model Program 5 - Runge-Kutta Integrator . Physical Properties Bibliography 140 158 169 170 175 182 189 206 209 215 LIST OF TABLES Table . Specific Impulse for Various Propellants. . Calculated Specific Impulses - Pressure Dependence. . Calculated Specific Impulses - Dissociation Dependence (Medium Pressure) . Calculated Specific Impulses - Dissociation Dependence (High Pressure) . Experimental Results. . Theoretical Thrust Efficiencies . Temperature Offsets for Collisional Relaxation - Recombinational Excitation Vibrational Steady State . Parameters for Modelled Run viii Page 20 21 22 58 63 78 85 10. 11. 12. l3. 14. 15. 16. 17. LIST OF FIGURES Figure . Schematic of Proposed Microwave-Plasma Electrothermal Propulsion System . Temperature and pressure dependence of equilibrium dissociation for H2 . Experimental Apparatus - Schematic - Plasma Generator . Schematic of Resonant Cavity . Cross-Sections of Plasma Tubes . Experimental Apparatus - Schematic - Thrust Stand . Gas-Dynamic Temperature vs. Flow Rate in Nitrogen - 6.4 cm - Counter-Current Air . Gas-Dynamic Temperature vs. Distance From Plasma In Hydrogen - A11 Flow Rates . Gas-Dynamic Temperature vs. Mass Flow Rate In Hydrogen - 6.4 cm Below Plasma Gas-Dynamic Temperature vs. Distance Below Cavity in Argon - 22 mg/s Flow Rate Percentage of Incident Power Absorbed by Cavity vs. Plasma Pressure in Hydrogen Specific Impulse vs. Mass Flow Rate In Hydrogen - Effect of Plasma and Tank Pressure Modelled Translational and Vibrational Temperature Profiles - Low Temperature Modelled Translational and Vibrational Temperature Profiles - X0 - 0.99 - Y - 1. (Low Pressure) Modelled Translational and Vibrational Temperature Profiles - X0 - 0.99 - Y - 1. (High Pressure) Modelled Temperature Profile - Adiabatic Gas-Phase Reaction Modelled Temperature Profile - Non-Reactive ix Page 18 25 28 29 32 46 51 53 55 60 73 74 76 86 88 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. Modelled Temperature Profile - Non-Reactive - Realistic Heat Transfer Modelled Temperature Profiles for Modelled Run Without Surface Recombination Modelled Temperature Profile - Reactive Lossy Surface Modelled Dissociation Profiles - Lossy Reactive vs. Inert Surface Modelled Temperature Profiles - Effect of Surface Reaction Energy Distribution Modelled Temperature Profiles - Hot Wall Assumption Modelled Temperature Profile - Hot Wall Assumption Modelled Temperature Profiles - Hot Wall Assumption, Diffusion Controlled Modelled Temperature Profiles for Modelled Run Plug-Flow Model - One-Dimensional - Fully Developed Modelled Dissociation Profiles - Effect of Surface Reaction Fully Developed Model Modelled Temperature Profile - Counter-Current Coolant Flow Modelled Temperature Profiles - Flow Rate Dependence Modelled Temperature Profiles - Initial Dissociation Dependence Modelled Temperature Profiles - Initial Temperature Dependence Modelled Temperature Profiles Initial Dissociation Effect Conductivity 0.1 Hydrogen - Modelled Temperature Profile - Non-Reactive, Conductivity - .15 Hydrogen vs. Nitrogen Experiments Modelled Temperature Profiles for Modelled Run - Two- Dimensional Model - Five Stations Modelled Average Temperature Profile - Low Flow - Two-Dimensional Model - Five Stations - Cool Sheath 91 94 96 97 98 100 101 103 104 106 107 108 110 111 112 115 119 120 36. 37. 38. xi Modelled Temperature Profiles - Axial Dispersion vs. Plug-Flow Modelled Dissociation Profiles - Axial Dispersion vs. Plug-Flow Modelled Temperature Profiles - Dispersion Model - Flow Rate Effect 129 130 131 NOMENCLATURE Cross-sectional Area of Flow Reactor Outer Surface Area of Cooling Jacket Surface Area of Reactor per Length Area of Nozzle Throat Heat Capacity of Reacting Mixture Molar Heat Capacity of Coolant Vibrational Heat Capacity Diffusion Constant, Atomic Hydrogen Through Molecular Hydrogen Heat of Recombination Internal Diameter of Reactor Change in Dissociation From Reaction Change in Dissociation From Volume Reaction Change in Dissociation From Surface Reaction Vibrational Energy Density Molar Flow Rate of Hydrogen into Plasma Fanning Friction Factor Molar Flow Rate of Coolant Standard Gravity Gravitational Constant Local Mass Flux Planck's Constant xii xiii Heat-Transfer Coefficient from Surface to Coolant Heat-Transfer Coefficient from Gas to Surface Mass-Transfer Coefficient Measured Specific Impulse Theoretical Specific Impulse Equilibrium Constant Boltzmann's Constant Rate constant for Vibrational-Translational Collisional Relaxation Surface Reaction Rate Constant Reaction Rate Constant, Hydrogen Atom as Third Body Reaction Rate Constant, Hydrogen Molecule as Third Body Molecular Weight Molecular Weight of Molecular Hydrogen Avogadro's Number Pressure Stagnation Pressure Heat Flux to Coolant Gas Constant Radius Volumetric Rate of Generation of Molecular Hydrogen by Recombination Temperature, Translational Temperature xiv TC Coolant Temperature Te Temperature Based on Total Relaxation of Vibrational Excitation Tg(e) Estimated Gas-Dynamic Temperature To Temperature of Environment TV Vibrational Temperature TVO Initial Vibrational Temperature Tw Surface Temperature TO Stagnation Temperature ub Bulk Velocity of Flow Uo Heat-Transfer Coefficient to Environment UO Overall Heat-Transfer Coefficient X Fraction of Hydrogen Dissociated Xw Dissociation at Surface Y Fraction of Recombination Energy Stored as Vibrational Excitation 2 Distance Below Plasma Subscripts c Coolant V Vibrational w Surface Conditions XV Greek Letters Specific Heat Ratio Thermal Conductivity Fundamental Frequency of Vibration Others Concentration of Species Chapter 1 - INTRODUCTION Most currently used or envisioned concepts for systems for space- craft propulsion and maneuvering rely on what may be called, in general, reaction propulsion systems. These consist of a means for accelerating a mass to some velocity relative to the spacecraft and releasing it. By Newton's Third Law, the spacecraft recoils, gaining momentum equal to that imparted to the released mass. If the released mass, also called reaction mass, is released in a continuous flow, or in frequent, short bursts, the spacecraft experiences a force, or thrust, equal to the product of the reaction mass flow rate and the relative velocity imparted to the reaction mass. The reaction propulsion system most commonly used in space flight, to date, is, of course, the chemical rocket, which uses a stream of gas as the reaction mass, heated under pressure by chemical reactions within the gas, then accelerated by thermodynamic expansion through a nozzle. Alternatively, the reaction mass, sometimes called the working fluid, can be heated externally, for example, as in the NERVA nuclear rocket, by passage through a nuclear fission reactor. A totally different con- cept, under development to the point of prototype flight tests, ionizes the working fluid and uses electrostatic forces to accelerate the ions - the ion engine. A third, under study, would use electromagnetic forces to accelerate solid bodies, then release them in rapid succession - the mass driver (O'Neill [1977]). These are only a few of the reaction propulsion systems that have received serious scientific and engineering study, while nearly every system that can perform the basic function has probably been considered. For example, a toy balloon, when blown up and released uses the elastic forces of its rubber skin and the pressure of the air held inside to force the air out and the reaction causes the balloon to fly around. A similar technique is exemplified by the toy rockets that are filled with water, then pressurized by a hand air pump. When the rocket is released, the air forces the water out, and the rocket flies off. One particularly memorable science fiction novel (Anderson [1962]) describes the use of an aqueous solution of fermented malted barley, pressurized by dissolved carbon dioxide, and expanded through a nozzle to provide emergency spacecraft propulsion. Much recent study has been directed at electric propulsion systems, ones in which the energy needed to accelerate the reaction mass is, at some point, in the form of electricity. These systems have certain advantages for certain applications. The performance of a chemical rocket is limited by the nature of the reactions involved, by the amount of energy available per unit mass of propellant and by the fact that the propellants have to include the products of the reaction. A nuclear thermal rocket, such as NERVA, has its own problems. An electrical propulsion system can be limited, in the amount of energy available per unit of reaction mass, only by the amount of energy and reaction mass flux available. In near-Earth orbits, the electricity could come from solar panels, making the potential energy store all but unlimited. For deep-space missions, nuclear reactors may be required to provide the electricity, but an electric propulsion system may be much more flexible than a nuclear thermal system, such as NERVA. Unlike the chemical rock- ets, non-chemical systems, such as NERVA and electric propulsion sys- tems, can be somewhat less selective in their choice of working fluid. NERVA, for instance, was designed to use hydrogen as a propellant, for its advantages in specific impulse (xlgg infra), to provide a large amount of thrust while using a small amount of propellant. No conventional chemical rocket can use hydrogen alone for reaction mass. Electric propulsion systems that have been under study group them- selves into three categories, depending on how the electrical energy is transferred to the working fluid. These groups are electrostatic sys- tems, such as the ion engine, electromagnetic systems, such as the mass driver and rail gun, and electrothermal systems, in which the electri- city is used to heat the working fluid and, as in a chemical rocket, thermodynamic expansion through a nozzle provides the thrust. An elec- trothermal propulsion system is, in fact, the first electric propulsion system to be used in practical application in space. This is the resis- tojet, in which the working fluid is heated by passing over and around a coil heated resistively, much as in an electric hair drier, before it flows out the nozzle. Certain of the latest series of INTELSAT commu- nication satellites use augmented hydrazine thrusters for station- keeping and attitude control. A standard hydrazine thruster uses the catalyzed decomposition of hydrazine to heat the decomposition products and force them through a nozzle to produce the thrust. The augmented thruster places a resistance heating coil before the nozzle, to further heat the working fluid electrically. The resistojet has been studied for use with a wide range of pro- pellants, besides the ammonia-nitrogen-hydrogen-hydrazine mixture used in the augmented hydrazine thruster. Early studies looked at using the waste products, water, carbon dioxide, ammonia, and methane, from a manned space station as propellants for the station-keeping and attitude control thrusters (Greco and Byke [1969], Halbach and Yoshida [1971]). Deposition and other reactions between the working fluids and the sur- face materials led to the concept's abandonment. The resistojet has other technical problems, involving the transfer of heat from the resis- tance coil into the gas. Since the coil must be hotter than the working fluid, the materials selection problem is most critical, to provide the mechanical strength and electrical performance required at the tempera- tures called for. To allow the coil surface temperature to be as little above the fluid temperature as possible, the flow path is often tortu- ous, to provide as much contact area as possible. This aggravates the problem of deposition. An alternative electrothermal propulsion concept, the microwave- plasma electrothermal rocket, is the subject of this work. In this system, shown schematically in Fig. 1, the electricity, from a solar panel, would power a microwave-frequency oscillator, and the microwaves generated would sustain a plasma in the working fluid, inside a resonant cavity, heating the working fluid, which would expand through a nozzle to produce thrust. This concept would have an important advantage over the resistojet, in that the working fluid would be heated from within, as in a chemical rocket, and the main unusual requirement for the plasma chamber wall is that it be transparent to microwaves. The heat and chemical resistance required are no different than for chemical rockets, and might be handled by techniques already standard for chemical rocket design and construction. An important mechanism by which energy might be coupled from the electromagnetic fields sustained by the microwaves to the working fluid is the dissociation-recombination cycle of a diatomic gas, such as hydrogen. In the plasma, a significant population of free atoms, Shaw SOLAR PANEL MICROWAVE F j OSCILLATOR —-> H2 TANK —>| ‘§~ NOZZLE CAVITY RECOMBINATION CHAMBER Figure 1. Schematic of Proposed Microwave-Plasma Electrothermal Propulsion System [1959] reported 90% dissociation in a hydrogen plasma, might be main- tained. On leaving the plasma, the atoms recombine, re-forming the molecules, and releasing the energy absorbed in their dissociating, and heating the gas still further. This mechanism could play such an impor- tant role that the system is sometimes referred to as the "free radical" [Hawkins and Nakanishi, 1981] or "atomic hydrogen" rocket. However, it should be clear that the system is, in its fundamental principles, electrothermal in nature. A commonly-used figure of merit in evaluating and comparing the performance of reaction propulsion systems is the specific impulse, which has units of seconds. Experimentally, the specific impulse of a given experimental rocket is determined by measuring the thrust of the rocket in pounds (force), and dividing by the flow rate of propellant in pounds (mass) per second. In the 81 system of units, the force in Newtons divided by the flow rate in kilograms per second must be divided by the standard acceleration of gravity, 9.8 m/sz, to give the correct units, and allow for the difference between the pound (force) and pound (mass) in the English system. An engine with a high specific impulse can produce a larger total impulse, or change in velocity, from a given amount of reaction mass than one with a lower specific impulse. Hence, all else being equal, as high a specific impulse as is consistent with other constraints is desirable in a spacecraft propulsion system, where total spacecraft mass is likely to be a critical constraint. The specific impulse can also be computed by dividing the exhaust velocity, in meters per second, by the standard acceleration of gravity. This permits estimation of the specific impulse for theoretical systems, by assuming that all the energy added to the propellant is converted to linear kinetic energy, and computing the velocity from E - % m V”. Table 1 lists the specific impulse for various common chemical rocket propel- lant mixtures, specifically those used in the first stage of the Saturn V and the Atlas boosters, the Titan and Ariane rockets, and the Centaur, Saturn V upper stages, and the Space Shuttle Main Engines, the design value for the NERVA nuclear rocket, and a range of values for various ion engine designs. Although the chemical rockets have lower specific impulse than the nuclear or ion engines, they will remain the propulsion system of choice for the foreseeable future, at least for surface to orbit operations. Neither the ion engine nor the NERVA nuclear rocket can generate enough thrust to lift themselves against Earth's gravity, though other nuclear rocket designs might get around the problems with NERVA, and be usable for direct launch (Kingsbury [1975]). Ion engines, and indirectly the NERVA atomic rocket design, are both limited in the amount of reaction mass flow they can handle for a given amount of system dry mass, thus limiting the total thrust they can generate. The NERVA design is limited by heat transfer from the reactor core to the working fluid/propellant. I Although the microwave-plasma electrothermal rocket is not purely a chemical rocket system, the importance of the dissociation-recombination cycle in the performance of the system, particularly when hydrogen is the working fluid, allows it to be treated as one, with the main reaction the recombination reaction of hydrogen: H + H - H2 (1) .The heat of reaction of this reaction is very high, 436 MJ/kg-mol of product at 298 K, and, combined with the low molecular weight of the Table 1. Specific Impulse for various Propellants. Propellant O2 - Kerosene °2 ' H2 N204 - UDMH - Hydrazine Nuclear (NERVA) Ion Atomic Hydrogen H:H2 Ratio 1:10 1:5 1:4 1:1.6 1:0.54 1:0 Specific Impulse (sec) 300 391 278 825 1500-8000 472 660 704 1,046 1,470 2,103 Reference Sutton Sutton Sutton Sutton Sutton Palmer Palmer Stewart Stewart Stewart Stewart [1976] [1976] [1976] [1976] [1976] [1959] [1959] [1962] [1962] [1962] [1962] product, 2.016, if all the energy released in the reaction is trans- formed to the kinetic energy of the exhaust, this suggests that the exhaust velocity from a rocket fueled by atomic hydrogen alone would be as high as 21 km/s, with a corresponding specific impulse of 2100 5. By comparison, burning hydrogen with oxygen releases 242 MJ/kg-mol, produc- ing an ultimate exhaust velocity of 3.7 km/s, due to water's molecular weight of 18, and specific impulse of 370 5. Table 1 also shows this sort of theoretical specific impulse for a rocket with a propellant of various proportions of molecular and atomic hydrogen, assuming total recombination of the free atoms. The table shows the promise of atomic hydrogen as a propellant: higher specific impulse than any commonly used chemical rocket, into the range of nuclear rockets and to match some ion engines.‘ It seems likely that a rocket based on the recombination reac- tion of hydrogen can avoid the political and technical problems that killed NASA's NERVA nuclear rocket project. In addition, due to the different concepts behind the two systems, it seems likely that the microwave-plasma electrothermal rocket can generate higher total thrusts than an ion engine can, at similar system masses. The electrostatic heart of an ion engine is so limited as to how much reaction mass it can accelerate per unit area of the grids used that it becomes very expen- sive, in money and in engine mass, to scale the engine up to generate more than very low thrusts. These low thrust levels mean that all maneu- vering must be leisurely, though the high specific impulse makes the ion engine ideal for long-duration, unmanned flights. NASA's proposed, though cancelled, probe to Comet Halley was intended to use an ion engine, for example. The microwave-plasma electrothermal rocket may be nearly so economical in reaction mass as as ion engine, with the added 10 benefit of easy scaling to thrust levels currently the regime of chemical rockets. Atomic hydrogen, because of the large amount of energy released in its recombination reaction, and the high specific impulse potentially derivable from this energy, has long been desired as a rocket propellant (Stewart [1962]). However, the normal state of hydrogen, under reas- onably accessible conditions, is as molecules, the equilibrium of the reaction in equation (1) is shifted far to the right. Hence, the free atoms must be generated from the molecular state. Atomic hydrogen is, hence, not truly a source of energy for propulsion, since the dissocia- tion of the molecules consumes as much energy as the recombination would release. This isn't without precedent, since hydrogen itself is not found free in large quantities on Earth, so the hydrogen used in the liquid-hydrogen fueled rockets used today, such as the Space Shuttle Main Engines, must be liberated first from some hydrogen-containing compound. The most easily accessible source of hydrogen is, of course, the water of the oceans. The electrolytic dissociation of water to generate hydrogen is precisely the reverse of the combustion reaction of hydrogen with oxygen, used in the rockets. Given that the hydrogen atoms must be generated before they can be used for propulsion or for other purposes, there are two obvious options as to how and where this can be done, in the case of a spacecraft with an atomic hydrogen-fueled engine. The atoms may be either generated on board, at need, from a supply of molecular hydrogen or generated at a base, either on-planet or in orbit, and stored on board until they were needed. Since most methods for generating atoms were both inefficient and heavy, early studies (Ibid., and references therein) concentrated on 11 the second option, of storing the free atoms. This proved difficult, as free hydrogen atoms are extremely reactive, under anything like normal conditions. A typical storage scheme involved freezing the atoms in a matrix of molecular hydrogen, at 20 K, at about one atom per ten mole- cules (Golden [1958]). A later storage method involves using magnetic fields to sort the spins of the elementary particles making up the atoms so that recombination is quantum-mechanically forbidden (Silvera and Walraven [1982]). This also involves very low temperatures, a major component of the system being super-fluid liquid helium, at less than 3 R. All storage methods have proven to be unwieldy and impractical for spaceflight, so the concept has been abandoned, and the concept of an atomic-hydrogen fueled rocket was neglected as impractical. The concept has been revived, however, in the wake of the results of Shaw [1959], Mearns and Ekinci [1977], inter alia, showing that the microwave plasma could be an efficient source of free hydrogen atoms. In addition, recent advances in electronics and solar-cell technology, driven by the space program, have produced such a reduction in size and mass of such equipment and improvements in efficiency that a microwave- plasma electrothermal propulsion system, as described above, now appears feasible for interplanetary propulsion. Although hydrogen's low molecular weight would produce the highest specific impulse, and hence might be preferred as a propellant, any diatomic gas could be used. It seems likely, from cosmochemical consi- derations, that oxygen might be the propellant of choice in the inner solar system. Although hydrogen is the most abundant element in the universe, as stars and giant planets like Jupiter are largely made of 12 hydrogen, and the ice moons of the giants are made largely of hydrogen- containing compounds, such as ammonia, methane, and water, the smaller, warmer, and hence rocky inner bodies of the solar system could not retain as much of the primordial hydrogen or hydrogenous compounds. Oxygen, the third most abundant element in the universe, after the inert helium, is the predominant element in the small, rocky bodies of the solar system, in its silicate, alumina, and other compounds. It would, no doubt, be a major by-product of any lunar or asteroidal mining and smelting operation. Hydrogen, on the other hand, is comparatively rare in the inner solar system, and would be far too valuable to throw away as reaction mass, being mainly brought from Earth, or the ice moons of Jupiter or Saturn, unless an asteroidal source of water is located. Research work has also been done using nitrogen (Asmussen, g§_§l. [1984], Whitehair, gt al. [1984]) and even helium (Whitehair, g§_al. [1986]) as working fluids in a microwave-plasma electrothermal rocket. As part of a study of the microwave-plasma electrothermal propul- sion system, and of other potential uses for microwave—frequency plasmas in propulsion, NASA-Lewis Research Center contracted with the plasma chemistry study group at Michigan State University to conduct experimen- tal and theoretical studies of the microwave plasma (Grant # NSC-3299). This paper will describe three aspects of this joint study: First, a series of experiments, performed at NASA-Lewis Research Center, Cleve- land, Ohio, to test the total concept and to measure the thrust that can be generated by a microwave-plasma electrothermal rocket, second, a theoretical examination of the last stage in the dissociation-recombina- tion cycle, the energetics of hydrogen atom recombination, and last, a 13 series of computer models of the recombination reaction to analyze and explain the experimental results. Chapter 2 - TECHNICAL BACKGROUND A number of technical, engineering, and scientific problems remain to be solved before the microwave-plasma electrothermal rocket takes any sort of role as a spacecraft propulsion option. The problems, in the order in which they're met along the energy pathway, include the effi- ciency of generating the microwaves, the transfer of the energy to the plasma, maintaining the plasma at operational flow rates and pressures, thermalizing the energy that may have been absorbed through non-thermal mechanisms, and controlling the heat of the exhausting working fluid. The first of these is a matter of electrical engineering, based on much prior work on microwave systems, and outside the scope of this work. The second may create interesting problems, for the plasma chamber required to pass a large flow rate of working fluid may be so much larger than our experimental systems that its primary resonant modes will be excited by radiation outside the microwave frequency band. How sensitive the plasma's behavior is to the frequency of the exciting radiation is a question for other work. It is clear that much of the detail design of the system is intimately related to the frequency that is used, and these considerations are also outside the scope of this work. The performance of the plasma under operational conditions will require investigation once those conditions are defined, though experi- ment and modeling, such as that reported in this work, will supply some preliminary guidelines. The physics and chemistry within the plasma are outside the scope of this work, though the experimental system used can contribute some data for these studies, and this modelling and experi- mental work, and extensions of them, can help define the direction in 14 15 which final design should lead. The final design of the thermalization section and nozzle of an operational microwave-plasma electrothermal rocket will be based on a detailed understanding of the processes at work in these sections, which will come through an interactive combina- tion of experiment and numerical modelling. For example, full dissocia- tion and complete recombination would heat the working fluid more than 15,000 kelvin. Just what temperatures would be attained and where, and how to handle them is a topic for a long design process, and this work is intended to provide a beginning for this process. This work will address some aspects of this phase of the overall problem, specifically, the energetics of the working fluid after it leaves the plasma to flow through the nozzle. This paper completes work published by Chapman, g; a1. [1982'], Filpus and Hawley [1984], and Morin, 341. [1982]. This work first reports a series of experimental work, at the NASA-Lewis Research Center, Cleveland, Ohio, to investigate the performance of a small-scale model of a microwave-plasma electrothermal rocket, attempt- ing to measure, or at least estimate, the temperature profile below the plasma and the thrust generated by the experimental system. Next, I consider the question of thermalization of non-thermal, vibrational and rotational excitation, energy transfer in the plasma. In the dissociation-recombination cycle, it has been suggested (Villermaux [1964a,b]) that a large fraction of the recombination energy may be manifest as vibrational energy of the recombined molecules for long times. This work employs a computer model to estimate the amount and the relaxation time for such excitation under expected operational condi- tions. Last, this work describes an attempt to construct a numerical 16 model of the recombination system, in an effort to predict the perfor- mance of proposed systems, and to verify the model by fitting its predictions to the experimental results. Chapter 3 - EQUILIBRIUM CONSIDERATIONS The recombination of free hydrogen atoms releases 436 MJ per kg-mol of molecular hydrogen produced. For this energy to be useful in a propulsion application, it must be turned to directed kinetic energy of flow of the gas. The object of this study is to determine the conditions under which the production of such kinetic energy, either directed, or random on the molecular level, which can become directed by thermody- namic effects in flow through a nozzle, can be maximized. The first, and simplest, approximation is to assume that a stream of hydrogen gas is dissociated completely, then allowed to recombine fully, with all of the energy released in the recombination serving to heat the molecular gas, which is then allowed to expand through a nozzle. This produces the figures in Table 1, an exhaust velocity of 20.9 km/s, and a specific impulse of 2100 seconds, as described above. However, this is equivalent to a temperature rise of some 15,000 K, from the kinetic theory heat capacity for molecular hydrogen of 7/2 R, 29,100 J/kg-mol K. Besides the staggering, though not totally insurmountable problems involved in handling gases at such temperatures, about two and a half times the temperature of the surface of the sun, thermal disso- ciation will occur, at any reasonable temperature, limiting the extent of the recombination. At 1 atmosphere pressure, the equilibrium disso- ciation at 5000 K is over 90 per cent (Fig. 2, from Snellenberger [1980]), for instance. The next step, then, is to consider the effects of chemical equi- librium on the extent of the recombination reaction. As a simplified model of the system, the inlet gas was assumed to be at 300 K, and 17 l8 80? u 000m OOON .m: toe cofiuovuommwu Earta_paaao co mucmccmamv meanness ecu ma:aocoa5¢p--.~ ocamau x6 .oLBmcanma. 89. -9 10m 10m .00 low 10% Jom 10v. uogsueAuoo we: Jed 19 allowed to react adiabatically to equilibrium, then no further reaction occurred in the nozzle, frozen flow. Again, from the kinetic theory of gases, the heat capacity of atomic hydrogen was taken to be 5/2 R, or 20,800 J/kg-mol K. The equilibrium constant was computed from the exact heat capacities (see appendix C). Table 2 shows the equilibrium tempera- ture, dissociation, and specific impulse produced by a fully dissociated hydrogen stream over a range of pressures. The specific impulses were computed based on the theoretical adiabatic expansion into a perfect vacuum, assuming no further reaction, or frozen flow (Equations 3a and 3b, p. 39). Tables 3 and 4 show the same values for a range of initial dissociations at pressures of 100 kPa, 1 atm, and 2500 kPa, 25 atm, respectively. These tables show that thermal equilibrium does limit the perfor- mance of an atomic hydrogen rocket short of that which a simpler model would predict, particularly sharply at high atom concentrations and low pressures. However, these limits are still high enough, even at low pressures and relatively low concentrations at higher pressures, to generate performance, in terms of specific impulse, above that of chemi- cal rockets and comparable to nuclear rockets. It should be clear from Table 1 that a specific impulse of 1000 seconds is a useful benchmark to shoot for, well above chemical rocket performance and comparable to nuclear rockets. A fully dissociated stream would produce this impulse at a recombination pressure of 10 kPa (1/10 atm), while at 100 kPa (1 atm), 30% dissociation suffices to generate the same level of specific impulse. 20 Table 2. Calculated Specific Impulses - Pressure Dependence. Adiabatic equilibrium recombination from fully dissociated 300 K Hydrogen stream. Pressure Equilibrium Equilibrium Specific Temperature Dissociation Impulse = 105 Pa) (K) (sec) .1 3559 .706 1200 .2 3707 .692 1223 .5 3920 .673 1258 1.0 .4096 .658 1285 2.0 4287 .641 1314 5.0 4564 .616 1356 10. 4794 .595 1389 25. 5131 .565 1437 21 Table 3. Calculated Specific Impulses - Dissociation Dependence (Medium Pressure) Adiabatic equilibrium recombination from 1 atm 300 K Hydrogen stream. Initial Equilibrium Equilibrium Specific Dissociation Temperature Dissociation Impulse (K) (see) .1 1707 7.5:1r10'5 731 .2 2656 .023 940 .3 3049 .086 1026 .4 3282 .161 1080 .5 3424 .241 1122 .6 3598 .324 1158 .7 3726 .407 1191 .8 3847 .491 1223 .9 3968 .575 1254 1.0 4096 .658 1285 22 Table 4. Calculated Specific Impulses - Dissociation Dependence (High Pressure) Adiabatic equilibrium recombination from 25 atm 300 K Hydrogen stream. Initial Equilibrium Equilibrium Specific Dissociation Temperature Dissociation Impulse (K) (sec) .1 1703 1.5*10'5 731 .2 2830 .009 973 .3 3465 .051 1099 .4 3850 .113 1176 .5 4133 .182 1233 .6 4365 .256 1281 .7 4570 .332 1324 .8 4760 .410 1363 .9 4944 .487 1400 1.0 5131 .565 1437 23 To improve further on this approximation, the finite reaction rate of the recombination must be considered. This is done by making a mathe- matical model of the reaction system, and solving the equations involved. A numerical solution, using a computer, is called for. Detailed description of such modelling is included in Chapter 6 of this paper, below. Chapter 4 — EXPERIMENTS As part of the NASA-supported study of the microwave-plasma "free radical" rocket concept, an experimental program was initiated at NASA- Lewis Research Center, Cleveland, Ohio. The portion of this program that pertains to this work had three primary goals. First, the measurement of thrust from an experimental model of a microwave-plasma rocket, and investigation of some of the parameters that might affect the thrust as measured. Increases in the measured thrust with the plasma over that meausred in cold flow, the magnitude of the increase and what parameters affected the increase were primary targets of the investigation. Second, determination of some of the conditions downstream of the plasma, for verification of numerical models. Specifically, temperature estimations 21; a gas-dynamic method were sought, especially if a profile of this temperature below the plasma could be established. Third, investigation of different configurations of the plasma tube, to decide which might overcome some of the limitations of the experimental arrangement, or would be better suited for practical applications. Apparatus The experimental apparatus consisted of two separate sections, the microwave-plasma generator, shown in schematic form in Figure 3, and the thrust stand, shown in schematic in Figure 6. In the microwave-plasma generation section, the gases used were commercial grade bottled gases, fed through micrometer variable leak control valves, with the flow rates 24 25 uououosou mammam - ouuofiosom - maueuunm< Huucoaauonxm .m madman mmDh m0m.P .>._<> ¥ mw02_._>o “3420 304k. museum”... .. .. . mozmmzmo \\ w><>>0m0=2 Ill ~ ll... -7-» I H ,1_ HM: ‘1 55.58 _ ES... / mwAdDOU 4420:0920 >._..._ .9: EDDO<> 26 measured by Hastings-Raydist flow meters, calibrated in air with conver- sion factors for other gases. The gases were fed to a flow manifold, where the gases were selected and mixed, and thence to the discharge tube. In transpiration flow experiments, the gas was fed to the jacket of the discharge tube, and flowed through the discharge tube walls into the discharge tube. Then, the gas flowed from the discharge tube through its nozzle into the tank of the vacuum facility. The thrust stand was located inside the vacuum tank, and the vane was positioned to intercept most of the discharge jet. The plasma pressure was assumed to be the pressure measured at the upstream end of the discharge tube by a Barocel capacitance manometer, connected to the discharge tube by a further eighteen inches -of rubber hose. A Bourdon tube manometer served to measure the jacket pressure in transpiration flow experiments, and to measure the upstream pressure of the venturi meter used to measure the flow rate of the cooling air. A water manometer was used to measure the throat pressure drop of the venturi meter, and the standard equation, with a calibrated efficiency, was used to determine the flow rate. Copper-constantan thermocouples were used to measure the temperature of the cooling air and, on some occasions, the temperature at the upstream end of the discharge tube. The alumina tube used for transpiration flow experiments was opaque, so a glass window was installed in the upstream end of the discharge tube to permit visual monitoring of the plasma. The microwave oscillator used had a frequency of 2.45 GHz, filtered to a narrow band, and the output power could be adjusted between 300 and 700 watts. The microwaves were passed through a water-cooled, ferrite core isolator, and a 40 dB directional coupler, where the power incident 27 on and reflected from the plasma cavity were measured, with thermistor- type power meters. About one meter of flexible wave guide allowed for the adjustments required in tuning the cavity, and led from the direc- tional coupler to the right-angle wave-guide-to-coaxial transition to the coaxial probe that led into the microwave resonant cavity. The cavity (Figure 4), of a Michigan State University-developed design described in Mertz, gt_§l., [1974, 1976] and Mallavarpu, g£_§l., [1978], was cylindrical, diameter 8 inches, with the discharge tube running coaxially through it. It was water cooled, though the water flow rate and temperature were not monitored. The length of the resonant section of the cavity and the depth to which the coaxial probe could be inserted into the cavity could be adjusted independently to tune the cavity to attain resonance, so the minimum amount of the microwaves incident on the cavity were reflected. A number of different configurations of the plasma tube were used in this study. The first one used was called a venturi tube (Figure 5(a)). Constructed of quartz, the design included an plasma chamber of about 22 mm. inside diameter, then it necked down sharply, to a throat of some 3 mm. diameter, then the inner wall expanded as a straight cone over some 23 cm. in length back to about 22 mm. diameter. There was a second quartz tube coaxial with the main plasma tube, over the plasma chamber and nozzle, with an overall outside diameter of 34 mm. This outer jacket was ported so that air could be blown over the inner tube, to cool it, with the air inlet above the microwave cavity, and outlets on the end of the nozzle, so that the air and the plasma gas flowed co- currently. 28 huw>mu occcomom mo ouumaonom .e ouswum mmeh 02.4000 mwh<>> hmOd 02_§w_> mwhmeOx‘ zeonatemozmr memo“. 0E 20.5535 may/F .lfil- m0h02m ousumuomaoa oHEm:ho-wmo .w muswwm .Eo. >t> - (2 c, To / MW) / g (3a) or I (ch) - (2 1 R T / {7 - 1) / M )1/2 / (ab) sp 0 w g The efficiencies given in Table 5 were computed by dividing the measured specific impulse by the theoretical specific impulse given by equations (3), assuming the zero-dissociation gas-dynamic temperature, and zero dissociation in the stream. The corrected temperatures and estimated. dissociations were calculated by assuming the nozzle efficiency remained constant from the cold to the hot flows, and solving equations (2) and (3), with the dependence of specific heat on composition, for the con- centration and temperature that provided the measured back pressure and specific impulse ratio. 57 Thrust Measurement Study An extensive study was done, at 600 W incident microwave power, using the 7-cm cavity mode, with hydrogen as working fluid over the range of flow rates 220 to 1400 pmol/s and corresponding plasma pres- sures of 12 to 57 Torr, with 4.7 g/s cooling air flow, the nozzle at the end of the tube, some 36 cm below the cavity, and into both ranges of tank pressure. The results are shown in Table 5, and the specific impulse is plotted vs. flow rate in Figure 12. The theoretical specific impulses are computed based on a non-dissociated hydrogen flow at room temperature, about 300 K, for the cold flow cases and at the gas-dynamic temperature for hot flow. There was some increase in measured specific impulse with flow rate, both with and without plasma, at high tank pressures. At low tank pressures, the measured specific impulse declined very slightly with flow rate. There was a slight kink, in that the very lowest flow rate showed a higher specific impulse than did the next higher ones. The measured specific impulse at low tank pressures was' from two to five times as high as that the corresponding flow rates- produced at the high tank pressures, when there was no such dependence expected. This was all independent of whether there was a plasma or not, and the ratio between the hot and cold flow values remained roughly constant, and comparable to the pressure ratios, consistent with using the gas-dynamic temperature as the gas temperature. At the low tank pressures, as can be seen in Figure 12, the cold flow specific impulse was measured to be higher than that predicted by an ideal, adiabatic, expansion into a perfect vacuum (as after Shapiro [1953]). Similarly, so was the hot flow value, if one assumed zero dissociation and used the 58 Table 5. Experimental Results. [a] High Back Pressure Gas Flow Rate (10 mol/s) Cold Flow: Thrust (MN) Pressure (Torr) Isp [m] (S) Isp [Th] (5) Efficiency Tank Pressure (nflbrr) Hot Flow: AbsorbedPower (W) Cooling Air: Flow (g/s) Outlet Temp. (K) mat Load (W) Pressure (Torr) Gas Dynamic Tauperature (K) Thrust (mN) Isp [m] (s) Isp [Th] (8) Efficiency 220 .30 10.0 72 290 .25 1.9 539 4.7 375 392 11.7 405 107 340 360 .52 15.0 75 290 .26 2.8 600 4.6 385 423 17.6 405 .62 89 340 .26 500 .72 20.0 74 290 .26 3.2 599 4.6 383 414 23.2 405 .91 93 340 .27 600 1.02 25.0 81 290 .28 3.6 600 4.7 382 407 28.9 395 1.17 94 336 .28 1100 2.41 40.0 113 290 .39 4.4 596 4.7 376 392 46.24 394 2.77 130 336 .39 1400 3.58 50.0 132 290 .46 5.2 597 4.7 374 384 57.5 391 4.07 150 334 .46 Table 5 (cont'd.) [b] low Back Pressure Gas_§low'Rate 220 (10 mol/s) Cold Flow: Thrust (mN) 1.44 Pressure (Torr) 9.8 3 Isp [m] (s) 3 7 I Th sp [ ] (s) 290 Efficiency 1.16 Tank Pressure .048 (mTorr) Hot Flow: Absorbed Pmer (W) 532 Cooling Air: Flow (g/s) 4.7 Outlet Temp. 376 (K) Heat load (W) 392 Pressure (Torr) 11.6 Gas MC 413 Taper-attire (K) Thrust (mN) ‘ 1.61 Isp [m] (s) 381 Isp [Th] (s) 344 Efficiency 1.11 Dissociation .183 Corrected 349 Gas whamic Tenperauu'e (K) 59 350 3.29 15.0 333 290 1.15 .074 548 4.7 377 397 17.3 408 2.64 385 342 1.13 .064 384 500 3.21 20.0 331 290 1.14 .14 556 4.7 378 400 23.1 404 3.67 378 340 1.11 370 640 4.14 25.0 333 290 1.15 .16 . 566 4.7 378 400 28.8 407 4.66 374 341 1.10 .175 346 1100 6.85 40.0 324 290 1.12 .24 558 4.7 374 384 46.0 398 7.78 367 337 1.09 364 1400 8.63 50.0 318 290 1.10 ' .3 553 4.7 373 376 57.2 393 9.74 359 335 1.07 .082 363 60 400 O O O o O O o "-0-3"; '0' - - - -;TT+EO'RETTCXL-VALTJE' T's-30.0%? e 300 J — - — - — —TH-E'0R'E-TICTI'. Vida-5505i .55 g 200 - 5' a. E D D U 2 D I 5,- CI I I ' 3 D I a. D l n U I 100 ‘ C] U E] I I I I I I 0 HIGH BACKGROUND PRESSURE, FLOW WITH PLASMA I HIGH BACKGROUND PRESSURE, COLD FLOW 0 Low BACKGROUND PRESSURE, FLOW WITH PLASMA 0 Low BACKGROUND PRESSURE, COLD FLOW I I I T 0 1 2 3 4 5 HYDROGEN FLOW (mo/Si Figure 12. Specific Impulse vs. Mass Flow Rate In Hydrogen Effect of Plasma and Tank Pressure 61 zero-dissociation gas-dynamic temperature to compute the specific impulse. SO, for that matter, was the hot flow value higher than that computed by assuming a constant nozzle efficiency to estimate the atom concentration and a corrected temperature, the last two lines in Table 5, since the cold-flow efficiency was over 100 per cent. The increased specific impulse, at constant gas-dynamic tempera- ture, at the higher flow rates and tank pressures, is most likely due to pluming effects, with higher flow rates producing narrower jets of gas. In addition, at low flow rates and high tank pressures, the flow through the nozzle can separate from the nozzle walls, and the environmental gas flows into and increases the boundary layer, thus decreasing the noz- zle's efficiency of expansion. A third effect may occur when the pres- sure is 1ow enough that free molecular flow occurs, and the recoil from the molecules bouncing Off the target may as much as double the thrust measured with a target-type thrust stand above that measured by a stand that measures the thrust imparted to the engine itself. Even these tests have been reported to show some dependence of the measured specific impulse on the background pressure. Since no effort was made to hold the' tank pressure constant, the vacuum system maintained different tank pressures at different flow rates. This might explain the observed flow rate dependence at low tank pressures and at the lowest flow rates at high tank pressures, as a small fraction of the difference between the two tank pressure ranges. Attempts to operate a transpiration flow tube proved unsuccessful. Thermal shock would crack the porous ceramic (alumina) tubes, changing the uniformity of the porosity. Breakdown in the inlet_jacket also proved problematic. Last, and most important, the plasma had a tendency 62 to migrate into the dead area of the tube above the transpiration flow region, so much of the gas would end up by-passing the plasma entirely. Discussion To gain some idea of the performance of a practical microwave plasma electrothermal rocket, based on the experimental results, it is useful to construct a theoretical engine concept whose performance may be extrapolated from the performance of the experimental system. For this theoretical performance, all of the energy that was lost to the cooling air, in the experimental system, is assumed to be confined to the working fluid, and heating it further, with all other energy losses assumed to be constant. Table 6 compares this theoretical thrust with the thrust measured in the experiment for a selected number of runs from those included in Table 5. For each case, the measured thrust and cool- ing air heat load are listed, with the measured specific impulse. Then, a computed specific impulse and thrust, with the ratio between the measured thrust and that computed, an efficiency, for three conceptual,- idealized, models of the engine. First, for ideal nozzle flow, assuming the zero-dissociation gas-dynamic temperature and zero dissociation in the working fluid. This efficiency is called the nozzle efficiency, though it obviously includes the characteristics of the thrust stand. Then, the ideal nozzle performance of the theoretical system, in which all the heat taken off by the cooling air is retained within the working fluid. Last, the performance that might be measured if the nozzle effi- ciency, the ratio of the the measured performance to the ideal, was the 63 Table6.'lheoretiml'mnstEfficianies mimental Results Mass TankPressure Thrust Specific CoolingAir Flow Regine Impulse Heat load (m/ S) (RN) (S) (W) ‘ . 436 High . 46 107 392 . 433 low 1 . 6 381 392 2 . 17 High 2 . 8 130 392 2 . 16 low 7 . 8 367 384 4 . 28 High 7 . 5 178 394 4 . 28 Low 14 . 7 350 376 Ideal Nozzle Flow Mass Tank Pressure Gas-Dynamic Predicted Specific Nozzle Flow Regnne Temperature Thrust Inpulse Efficiency (ms/S) (K) (mN) (S) (’3) .436 High 410 1.5 350 31 . 433 low 420 1. 5 355 107 2 . 17 High 400 7 . 4 347 37 2 . 16 low 405 7 . 4 349 105 4 . 28 High 395 14 . 5 346 51 4.28 low 395 14.5 346 101 Table 6 (caxt'd.) Mass Tank Pressure Predicted Specific Regine Flow (119/S) .436 .433 2.17 2.16 4.28 4.28 64 We High Low High Low High Low Thrust Impulse (mN) (S) 18.7 4338 18.7 4354 42.7 1999 41.4 1954 60.1 1426 58.7 1395 eoretica - Real Nozzle Efficiency (’6) 2.5 8.8 6.5 19.8 12.5 25. Mass Tank Pressure Predicted Specific Efficiency Flow Regine Thrust Impulse (DU/S) (mN) (S) (*5) .436 High 5.6 1326 8.1 .433 Low 20.0 4672 8.2 2 . 17 High 15. 9 749 17 . 2.16 TOW“ 43.6 2055 18. 4 . 28 High 30.8 734 24 . 4.28 Low 59.1 1441 25. 65 same in the theoretical system as in the real one. Since these effi- ciencies are the ratios of the thrusts, and the energy involved is proportional to the square of the thrust, squaring the given efficiencies will give the efficiencies in energy terms. Estimating temperatures from these theoretical performances, by solving Equation 3a for To, gave values ranging from 6500 K at the highest flow to 63,000 K at the lowest, with no dissociation. Under the conditions of the experiment (Table 5), adiabatic thermal dissociation would bring the atom-molecule mixture to equilibrium at about 30 per cent dissociated and 3100 K at the highest flow rate/pressure, and full dissociation, 100 per cent atoms, at 34,000 K at the lowest flow and pressure. This would reduce the ideal nozzle specific impulse, to 1070 seconds in the first case and to 3800 seconds in the second. These extreme predicted temperatures, particularly at low flow rates, demon- strate the "overkill" inherent in the experimental arrangement. The microwave flux incident on the cavity is so high that, especially at low flow rates, there is no practical mechanism by which the working fluid can absorb more than a small fraction of the incident microwave energy. . This may explain, to some degree, the need for cooling the plasma cavity, by showing how much of the microwave power absorbed by the cavity cannot be taken off by the working fluid without drastic effects. However, it seems more complicated than that, since plasmas have been maintained for long periods at low flow rates, and incident power levels similar to those used in the experiments described above, without any cooling of the plasma Chamber. On the contrary, the cooling air jacket can be sealed off and evacuated, in an attempt to insulate the plasma 66 chamber and the recombination zone from their surroundings. Even so, with that path by which energy might leave the plasma cut Off, and all the energy absorbed by the plasma supposedly remaining in the working fluid, the performance never approached the theoretical levels predicted above. There was, at best, a negligible improvement in the performance over that in cooled plasma tubes. How the energy is lost from the plasma gas remains unexplained. The observed heating of the outer wall of the vacuum jacket and inside the resonant cavity suggests that the vacuum jacket is not the insulator it was expected to be. Perhaps there were radiative losses, though in what wavelength is unknown. There was not enough visible light given off by the plasma to account for more than a small fraction of the lost energy. Another consideration comes from the manner in which plasma tubes fail while holding plasmas. At low flow rates and pressures, even vacuum-jacketed tubes can hold a plasma for long periods, with no maximum time yet found, under some conditions. However, at high flow rates and pressures, the tube fails quickly. The power that must be absorbed by the walls is (very slightly) higher at low flows than at high, for the same absorbed power. Hence, it might be' expected that the tube walls might fail at least as fast in low flows as in high. The greater stress due to the greater pressure at high flows may lead to faster failures at high flow rates, however. Air-cooled tubes have also failed, though, at very high plasma pressures and flow rates, though the pressure was still very much lower than the over atmospheric pressure in the cooling air jacket. Since the differential pressure across the wall of the tube is actually less at higher plasma pressure, when the jacket is open to the atmosphere, the stress on the wall in a cooled tube is less, hence it should not fail so fast. What 67 appears to happen is that the plasma shrinks at higher pressures, and the energy transfer to the wall must pass through a smaller area. The quartz of the tube walls cannot conduct the heat away fast enough, without heating up at the contact point, and the smaller the area the hotter it gets, and, eventually, it melts. There are a number of ways to improve the experimental design to avoid the problem of overcooling the gas. The microwave power level must be better matched to the capacity of the gas. The cavity resonant mode and the gas flow pattern can be selected to reduce the danger to the tube walls. One potential method, transpiration flow, was tried and could not be made to work. The distance between the plasma and the nozzle can be reduced to reduce the downstream cooling, it is to be hoped, without reducing the accuracy of the thrust measurements. A possible alternative design, for both experiment and applications, is regenerative cooling, where the feed to the plasma is ducted past the nozzle to cool the nozzle and pre-heat the feed. This is a common design in more conventional chemical rocketry. Another possible experimental arrangement, under discussion when this work was being done, would be to mount the entire cavity, plasma tube, and nozzle in a vacuum tank, on a thrust stand. This would help diagnose the origin of the back-pressure dependence, by determining whether the effect is due to the interaction between the impinging jet and the vane. It would also help reduce the effect of cooling by reducing the distance that would be required between the plasma and the nozzle. Chapter 5 - VIBRATIONAL ENERGY MODELING An important consideration in evaluating the practical utility of the microwave-plasma electrothermal propulsion system is the rate of energy release in the recombination reaction. A side factor in this consideration is the form which the energy, when released, should take. On the microscopic scale, this energy can take the form of translational kinetic energy, vibrational energy, or rotational energy. The latter two are also called internal energy of the molecules. This internal energy is not to be confused with the thermodynamic internal energy of the gas, which includes the translational energy on the molecular level. The first is most useful, as what is usually called heat is the randomly directed kinetic energy of the molecules of a substance, and thermo- dynamic effects transform this random motion into directed motion in a nozzle, to produce thrust. The conventional assumption in chemical reaction engineering is to assume that the energy is immediately manifested as heat, and that the internal energy of the molecules is either negligible, or reaches equi-- librium, thermalizes, quickly. However, the exact pathway of the trans- formation of energy from chemical potential energy through internal energy to translational energy greatly influences the rate of the reac- tion, as well as the rate at which the reaction generates useful energy. Theoretical models (Whitlock, g£_gl., [1974]) of the hydrogen recombina— tion reaction propose that highly vibrationally-rotationally excited states may serve as transition states, and experimental reaction rate determinations (Trainor, et 81., [1973]) provide some support to the models. Simulations of the recombination reaction (Levitskii and Polak, 68 69 [1973]) suggest that, at very high temperatures, over 80 per cent of the recombination energy may be manifested as internal energy of the recom- bined molecules. In addition, experimental results (Villermaux, [l964a,b]) have been interpreted as showing a "high" degree of vibra- tional excitation, and "slow" relaxation to translational energy, in the recombined molecules. However, the extent of the excitation and the rate of relaxation were not clearly quantified. This could be significant, because, if the recombined molecules still possess a significant fraction of the energy of recombination as internal energy when they reach the nozzle exit, the performance of the rocket could be seriously degraded from that predicted based on the normal assumption of total thermalization of internal energy. A study was therefore done to examine the effects of vibrational excitation on the performance of a hydrogen-propellant, microwave-plasma, electrother- mal propulsion system. Rotational energy, since the levels are so much more closely spaced than the vibrational levels, was assumed to thermal- ize much more quickly than the vibrational energy, hence vibrational relaxation would provide an upper bound for the total relaxation Of' internal energy. One of the models of the recombination section, below the plasma, was modified to account for above-equilibrium vibrational excitation of the recombined atoms, and for collisional relaxation of this excitation. To model the vibrational excitation, the hydrogen molecular gas was modelled using two temperatures: the translational temperature T, which describes the distribution of translational kinetic energy; and the 70 vibrational temperature TV, which describes the distribution of vibra- tional energy, as described by a simple one-dimensional oscillator model. The vibrational energy Ev is given as a function of T by the following equation (Herzberg [1950]): h u / k T Eva) - NAv h u / (e - 1) (4) where NAv is Avogadro's number, h is Planck's constant, v is the ground- state frequency of the oscillator, and k is Boltzmann's constant. Differentiating this with respect to T gives a vibrational heat capacity CPV: 2 ch v / k T cPVm - (Eva) / T) / R (5) where R is the gas constant. The energy of recombination was computed based on both T and TV by: DHR(T,TV) - DHR(0.,0.) + 3/2 R T - EV(TV) (6) with DHR(0.,0.) computed so that DHR(298,298) agrees with the literature values for the standard heat of dissociation for hydrogen. If it is. assumed that the energy of recombination is distributed between transla- tional and vibrational energy so that Y DHR(T,TV) goes into vibrational energy, and (l - Y) DHR(T,TV) goes into translational energy, where Y is a parameter between 0 and 1, the mass, translational energy, and 71 vibrational energy balance equations in an adiabatic, gas-phase, plug- flow reactor can be written as follows: F dX/dz + r A - 0 (7) H2 -F [(1 - X) CPH 2 + 2 x CPH] dT/dz + {(1 - Y)(-rH2)[-DHR(T,TV)] + kVT [EV(TV) - EV(T)]} A - 0 (8) 'F (1 7 x) CPV dTv/dz + {Y ('rHZ) [-DHR(T’TV)] - kVT [Ev(Tv) - Ev(T)]} A - 0 (9) where F is the molar flow rate based on undissociated hydrogen, X is half the flow rate of hydrogen atoms divided by F (the fraction disso- ciated), A is the cross-sectional area of the reactor, 2 is the distance along the reactor, r is the rate of production of molecular hydrogen H2 per unit volume by recombination, and kvT is the rate of collisional vibrational-translational relaxation, per unit volume. The adiabatic inviscid-flow reactor model (Eqns. (13), (15), and (16), p.83, below), with Equation (9), properly transformed, added and. Equation (15) modified to: 2 {mo ub / gc [l / (l + X) dX/dz - dA/dz / A] + (1 - Y) DH dX/dz + A [ (T ) - (T)] / F) dT/dz - R kVT EV V EV (88) (mO ui / gc T + CP) was solved using Computer Program 1 in Appendix B, via a Runge-Kutta-4 method, as described in Appendix A. The parameter Y was treated as a constant, user adjustable parameter, over each run of the computer program. 72 Figure 13 shows the output from a number Of simulations using this model, using the parameters of the modelled run (Table 8, p. 85), with -5 3 initial temperature 300 K and values of Y of 0., 10 , and 10- . The Figure shows the translational temperature profile for all three cases (solid curve), indistinguishable at the scale of the Figure, and the vibrational temperature profiles for the three values of Y and corres- ponding initial values for vibrational temperature of 300, 400, and 600 K (shortest, medium, and longest dashed curves, respectively). The initial values of the vibrational temperatures were selected to minimize the steepness of the initial portions of the curves, since any value selected would very quickly reach the same profile, as seen in Figure 14. It is also expected that, with finite values of the Y parameter, there would be some population of excited molecules leaving the plasma. The dissociation within the plasma would be accompanied by recombina- tion, the plasma conditions merely throwing the equilibrium towards the atoms side of the reaction. Thus, at finite values of Y, there would be some vibrational excitation generated within the plasma. The shortest_ dashed curve is seen to drop below the solid curve, and so do both other dashed curves, where they meet, though the scale of the Figure is too large to show this. This is due to the fact that the heat of recombina- tion includes some vibrational energy when assuming equilibrium ther- malization, hence there is an equilibrium value of the parameter Y, and if the value of Y input to the program is less than this equilibrium value, the vibrational temperature will lag behind the translational temperature. Temperature (K) 73 700 , -— Translational i w. Y- 0. l j —- TV. Y= 1.8-5 5 -- N. Y- Lil-S i 800—4 .4 5004 .1 400-1 300-4 i J 200 *T ‘1— T 1 m ' _T_ r T II T T fi— fi— fi“ j 0.00 0.05 0.10 0.15 Reactor Length (mm) Figure 13. Modelled Translational and Vibrational Temperature Profiles — Low Temperature 74 -_... 1340 J 1320 -1 1800 a] 1280 -J p N G O l L I Temperature (K) I ,..--------------___ \ I 1240 4 l-~ w 1220 -—< i 1200 T; "' TV. No - 8000 x '7 TV. No -= 1200 K -- Translational T fi— 1— r °‘° 0.5 1.0 Reactor Length (le—B m) 1180 Figure 14. Modelled Translational and Vibrational Temperature Profiles - X0 - 0.99 - Y - 1. (Low Pressure) 75 Experimentation with the program, at very short reactor lengths, showed that TV very quickly, on the order of microseconds in residence time and microns of reactor length, reaches values close to the curves shown if the assumed initial value is varied. This is shown in Figure 14, where the two curves, (b) and (c), vary only in their initial value for TV. To save computational time, Tv0 should be selected so that dTV /dz is zero at the start of the reactor, which is an option in the program, selected by feeding an initial value of TV less than 0. Figure 14 shows the output for the modelled run (Table 8, p. 78), with initial dissociation 0.99 and parameter Y - 1.0, showing the translational temperature profile (solid curve), and vibrational temperature profiles for initial vibrational temperatures of 1200 K (longer dashed curve) and 3000 K (short dashed curve). The small increase in the translational temperature reflects the small degree of recombination that takes place in the very low residence time of the simulated length of reactor. The vibrational temperature has, very quickly, reached a very nearly constant value, about 40 K above the translational temperature. Figure- 15 shows the translational (solid curve) and vibrational (dashed curve) temperature profiles for the same parameters as in Figure 14, with initial vibrational temperature of 3000 K and at a pressure a hundred times higher, some 3.5 atm. Within the 1.5 nanometers of reactor length plotted in the Figure, the vibrational temperature drops to about 70 K above the translational temperature, and slowly continues to drop. Figures 14 and 15 show that, at high temperatures, the vibrational temperature very quickly reaches a value somewhat above the transla- tional temperature, even with the assumption of Y - 1.0, and thereafter Temperature (K) 76 -- Vibrational 1 —- Translational .-rc.——...'-_--.. -.- .._. _l_ N 5.5 O O L 1 I l I l i l l l l l . ' l l l l l i l l l l l l L “00.1,. T r . f 4 I 0.0 0.5 1.0 1.5 Reactor Length (nm) Figure 15. Modelled Translational and Vibrational Temperature Profiles - X0 - 0.99 - Y - 1. (High Pressure) 77 the difference between the two temperatures decreases very slowly. This relaxation takes place on a time scale on the order of .1 usec, compared to the msec time for the reaction and residence time in a recombination reactor. It appears that the vibrational temperature has reached a quasi-steady-state value, at which the rate of excitation from recom- bination is almost exactly balanced by an equal rate of collisional relaxation. At very much lower temperatures, the initial difference is much larger, as at small Y it is very noticeable, and the difference drops much more slowly to the steady-state value, on a msec time scale. This suggests that an important value is this steady-state value, or offset, that difference between the vibrational and translational tem- peratures at which collisional relaxation transfers energy out of the vibrational mode as fast as the recombination reaction adds energy to the vibrational mode. This is determined by setting dTV/dz in equation (9), above, to 0, and solving the following equation for T at a given V T, X, P, and Y: {Y (-rH2)[-DHR(T.TV)] - kVT [EV(TV) - EV(T)]} A - 0 (10) which is transformed to: {Y -rH2)[-DHR(T.0)] + kVT EV1 / ( Y (-rH2> + kvr’ - Ev(Tv> (11) Table 7 gives the results of solving this equation, based on cer- tain assumptions. The offsets T -T as given in the table are linearized, V assuming that: EV(TV) - EV(T) - CPV(T) (TV-T) (12) 78 Table 7. Temperature Offsets for Collisional Relaxation - Recombinational Excitation Vibrational Steady State \ T(K) 300 500 700 1000 Elf—T 1:12 __"I‘ TgiI' 11,11" Tgil‘ Tyil‘ Igir YP YP YP YP . YP YP YP YP (K/Pa) (K/Pa) (K/Pa) (K/Pa) (WHK/Pa) (K/Pa)(K/Pa) .01 7400. 1.8e-3 .020 5.1e-6 8.4e-5 3.3e-7 1.4e-6 3.6e-8 .15 1.7e+5 .033 .63 1.3e-4 2.7e-3 8.7e-6 4.6e-5 9.5e-7 .30 3.5e+5 .052 2.1 3.4e-4 9.0e-3 2.38-5 1.6e-4 2.5e-6 .50 4.8e+5 .048 6.2 6.6e-4 .027 4.5e-5 4.7e-4 5.0e-6 .70 5.6e+5 .031 17. 9.8e-4 .074 6.8e-5 1.3e-3 7.6e-6 .90 6.28+5 .011 65. 1.2e-3 .31 9.1e-5 5.5e-3 1.0e-5 1.00* 6.2e+10 0. 5.9e+7 O. 3.8e+6 O. 5.9e+5 O. * As X goes to 1.0, total dissociation, offsets becate independent of pressure. Table gives offsets ocuprted near X = 1.0 at P = 100,000 Pa. TvgoestoinfinityatfiniteTe-T. Notes: aen = a x 10n . Offsets linearized, assmning constant CPV at value at T. True TV values lower than table predicts. For exatrple, at P = 100,000 Pa, T = 300 K, X = .01, Y = 1.0, table predicts TV-T = 7.4 x 108 important at low tenperatures. K. Exact TV value 2550 K. Most 79 This is valid here only if kVT is much greater than Y (-r Under this ). H2 assumption, these offsets are directly proportional to [Y (-rH )] / kVT’ 2 hence the table gives values of (TV-T) / Y P. This linearization, par- ticularly at low temperatures and high values of Y, greatly overesti- mates the vibrational temperature at steady state. Also, since kVT is proportional to the molecular hydrogen density and to the total gas density, whereas the recombination reaction is three-body, at low atom concentrations and far from chemical equilibrium, the offsets are directly proportional to the pressure. Hence, the table gives values of (TV-T) / Y P. At total dissociation, though, kVT goes to 0, and the offsets become independent of pressure and of Y. The table gives linearized offsets (TV-T) computed for values of X close to 1.0 at 100,000 Pa, approximately 1 atmosphere pressure. Particularly at lower temperatures, the linearization with pressure may not exactly hold at lower dissociations, as well, and the numbers in the table were computed for a pressure of 100,000 Pa. The second column under each temperature. in the table gives an estimate, similarly linearized, as to the total energy tied up in vibrational excitation at the steady-state, in terms of a temperature Te’ which is the temperature the gas would have if all the excess vibrational energy were turned to translational energy. Since the offsets in Table 7 are given divided by Y, they can be interpreted as the offsets for Y - 1.0, and thus give the upper bounds for the offsets in any real case. The table shows that the offsets go up with fraction dissociated, as the recombination rate increases. However, 80 the Te-T Offset peaks, and drops to O. at full dissociation. This comes from the reduction of the population of hydrogen molecules as dissocia- tion nears completion, and the vibrational temperature must be infinite for there to be a non-zero excitation energy with no molecules present. The offsets drop sharply with temperature, due largely to the increase in the rate of exchange relaxation, the dominant mechanism for colli- sional relaxation when there are atoms present. In exchange relaxation, a free atom strikes a vibrating molecule, and knocks one of the atoms in the molecule free, forming a less excited molecule with the remaining atom. Additional reduction with temperature comes from the reduction of the recombination rate, in both the rate constants and the concentrations of the atoms and of the third bodies. In conclusion, the conditions necessary to reduce the influence of the internal energy Of the recombined molecules on the macroscopic energetics of recombination are near-equilibrium dissociation, high temperatures, and low pressures. Under these conditions, any excess vibrational energy will relax down to thermal equilibrium, or as near as_ the recombination excitation will allow, very quickly, and the steady- state value will be close to true thermal equilibrium. In a practical microwave-plasma electrothermal rocket, however, high operating pres- sures are likely to be the rule, to force the recombination reaction towards completion, and generate high temperatures. However, the pres- sure dependence is much less than the temperature dependence. Under conditions that are somewhat conservative, for these considerations, for a practical electrothermal propulsion system, 3000 K temperature, and 20 atmospheres pressure, the Te-T offset will be less than 0.3 K and the 81 relaxation to this value is likely to be on the order of microseconds, a negligible concern in the design of the system. The assumption of imme- diate thermalization will be sufficiently valid for accurate design. Villermaux's experiments [1964a,b], which reported high vibrational excitation and slow relaxation, were performed at very much the other end of this regime. The temperatures were low, around 300 K, and the pressures were relatively high, 30 to 300 Torr. Under these conditions, very sensitive instruments may detect the temperature Offsets of the magnitude calculated here. Chapter 6 - COMPUTER MODELLING In order to help design further experiments, and potentially to predict the performance of practical microwave-plasma electrothermal rocket systems, development of a numerical model of the recombination section, downstream of the plasma, Of such a system was initiated. The recombination section is modelled as a gas-phase chemical reactor, with the reactants being a mixture of atomic and molecular hydrogen, and the main reaction being the reversible recombination of atomic hydrogen. The reaction rate constants and other physical constants used and their sources are described in Appendix C. The model has gone through a lOng and somewhat evolutionary process of testing and relaxation of assump- tions, in attempting to fit the results computed from the model to the experimental data, to ascertain what factors need to be considered in the model, to verify the model's applicability, and to justify confi- dence in the model's predictions for further work. The reactor described by the model was an adiabatic, gas-phase, one-dimensional, plug-flow model at first. To make the model more realistic, and to correct the' discrepancies found between the model and the experimental results, cooling, surface reaction, radial diffusion, making the model two- dimensional, and last, axial dispersion were added in turn. The equations were solved, at first, as an initial value problem by a specially-written version of a fourth-order Runge-Kutta technique (see Appendix A), then, as difficulties with the dispersion model turned up, other implicit techniques were tried, culminating with the IMSL DGEAR subroutine. 82 83 Adiabatic, Gas-PhaseI Plug Elow Model As described in Chapman, et a1. (1982), at that point the model was of a one-dimensional, adiabatic, gas-phase, inviscid flow, constant pressure reactor. The equations for the model were as follows: dX/dz - A (k, [H] + k, [H2]) ([H2] - K [H12> / F <13) dT/dz - - DHR dX/dz / CP (14) Inviscid flow pressure effects were added, after Shapiro [1953], for some inconclusive investigation of nozzle flow. This changed equation (14), above, and added another equation, as follows: 2 {m / g [l / (l + X) dX/dz - dA/dz / A] + DH dX/dz} dT/dz - 0 uh c g (15) (m0 ufi / gc T + 6,) dP/dz - - mo ui / gC [dT/dz / T + dX/dz / (1 + X) - dA/dz / A] x P / [R T (l + X)] (16) However, the singularities in the flow equations at sonic velocity' proved intractable, so the investigation was abandoned. Viscous losses 84 were added to a model for a cooled reactor, changing the energy and mechanical energy balance equations to: 2 {m0 ub / gC [dX/dz / (l + X) - dA/dz / A] + DHR dX/dz - dT/dz - u: f / [01 R T (1 + X)] + 00 AS (TC - T) / F} (17) (m0 u: / 8c T + CP) dP/dz - {- mO u: / gC [dT/dz / T + dX/dz / (1 + X) - dA/dz / A] - f ufi / 8. 0,) P / [R T (1 + X)] (18) The model considering viscous losses predicted, in the straight—tube geometry of the experimental apparatus, pressure drops much less than one percent of the total system pressure over the length of the reactor, so the viscous loss terms were neglected in further investigation. This model, equations (13), (15), and (16), with parameters given in Table 8, yields the temperature profile shown in Figure 16. This modelled profile is so incompatible with the experimental profile shown in Figure 11 that is is Obvious that the assumptions used in developing. the model need to be reevaluated. 85 Table 8. Parameters for HOdelled Run Modelled Gas Hydrogen Initial Conditions Fraction Dissociated 0.10 Gas Temperature (K) 1200. Gas Pressure (Torr) 26.5 Coolant Temperature (K) 400. System Parameters Gas Flow Rate (mg/s) 0.77 Tube Internal Diameter (cm) 2.2 Coolant (Air) Flow Rate (g/s) 5.0 Coolant Flow Direction co-current Adiabatic Equilibrium Conditions Fraction Dissociated 0.02 Temperature (K) 2300. 86 mm thlrln nofloeom 00057006 0309084 0595 05508.38“. “0:000: .3 RENE A83 nausea 8503* on an FthPlFLPLLFL ON —lP[Ph b ma _ OH 0 FLLlI—LkhnLrh b L h l.‘r...i I“.ll..lli v.‘ I I! coo.“ loom“ T and.“ Tone“ loam” T88. CONN ()1) engenders], 87 Cooled Model The change in the experimental profile with the change in coolant flow geometry argues for the inclusion of cooling of the modelled reactor, changing the energy balance equation to: dT/dz - {1110 u: / gC [dX/dz / (l + X) — dA/dz / A] + DHR dX/dz + UO As (TC - T) / F} / (m0 ufi / 80 T + c?) (19) Simulations, by trial and error, assuming zero dissociation, constant coolant temperature and constant overall heat-transfer coefficient, produced some fairly good fits to the experimental results. In Figure 17, the solid curve, which corresponds to an initial gas temperature of 700 K, coolant temperature of 370 K, and overall heat-transfer coeffi- cient of l W/mzK, is seen to fit the data from 6 cm out fairly well. However, the initial temperature is much lower than that measured at the exit of the cavity. The dashed curve in Figure 17 assumes a more realis- tic gas temperature of 1200 K, and requires a heat-transfer coefficient_ of 3.7 W/mzK to bring the temperature down to match the experimental 600 K at 6 cm, and a coolant temperature of 400 K to keep the temperature at 36 cm comparable to the experimental values. However, in between 6 and 36' cm, the model predicts a lower temperature than the experiment recorded, and the 400 K coolant temperature is somewhat higher than that used in the experiments. To see whether these heat-transfer coefficients are indeed reason- able, the fluid-flow regime must be determined. Assuming the viscosity of the plasma is roughly that of normal, equilibrium-dissociated 88 ebuaouomldoz 0595 ouuuuuomaoa voaovoz .5 0.5th A83 new—Hog acaouom on , 8 mm em 3 S m o ~bPF~brthanrb—LbF~I~rL~b—bel—|—LIL£ :::::::::::::::::::::::::::::: :o: 7 o I’m/f Jana o hunch II a u DD .1... III... ‘li‘l'i'zlllivu‘IlI‘i‘b. l":‘ III! } lull. I, I’....ul|l.|ll I’lai IIIIIIII‘II.‘IIIIIIII.’.II. ‘ can.“ 69: ()1) amquzedmal 89 hydrogen, reasonable at low dissociations, as the viscosity of atomic hydrogen is 0.7 times that of molecular hydrogen (Browning and Fox (1964)), close enough to 1.0 for these purposes, at the highest flow rate used in the experiments, 0.77 mg/s, that used in the sample run, the Reynolds number is about 5 at 300 K, and about 12 at 1200 K. This put the flow clearly in the laminar regime. The question of entrance effects must be considered, as the flow profiles in and exiting the plasma are unknown. Entrance effects in heat transfer are correlated by the Graetz number, the product of the Reynolds number, the Prandtl number, and the ratio of the diameter to the length of the tube. Theo- retical analyses indicate that the Nusselt number reaches a limiting value of 3.66, in the case with constant wall temperature, for Graetz numbers less than 4.0 (Perry and Chilton, [1973], pp. 10-12,10-13). Since the Prandtl number for simple gases is about 0.7, at 1200 K, entrance effects are smoothed out after about 2 tube diameters, 4 cm in the experimental system, so the flow may be assumed to be fully devel- Oped laminar flow, with this constant Nusselt number, with reasonable accuracy. The Graetz number dependence is such that the Nusselt number, and hence the heat-transfer coefficient, is infinity at the entrance to the tube, and decreases with tube length. This may suffice to produce the shift in heat-transfer coefficient predicted by the fitting of the model to the data. However, the entrance zone is very short, and the flow through the plasma is probably laminar already, hence, this effect is probably minor. At a constant Nusselt number and tube diameter, the heat-transfer coefficient is directly proportional to the thermal con- ductivity. Literature values for the thermal conductivity of hydrogen 90 (Lbig., p. 3—215) show a marked increase as the temperature goes up, by about a factor of 2.5 over the range from 400 to 1200 K. However, after correcting for an apparent error in the table in Perry and Chilton ([1973], p. 3—215) (see Appendix C), the internal heat-transfer coeffi- cient, in the experimental system, for a Nusselt number of 3.66 should range from 37 W/mzK at 400 K to 88 W/m2K at 1200 K. From annulus flow correlations (lbid., p. 10-15) and the flow rate and properties of the air used as coolant, the heat-transfer coefficient in the cooling jacket is calculated to be around 75 W/mzK. Thus, the overall heat-transfer coefficient should range from about 25 W/mZK at a gas temperature of 400 K to about 40 W/mZK at a gas temperature of 1200 K, well above the values that produced an apparently good fit. The temperature profile for the modelled run at zero dissociation and with a constant Nusselt number of 3.66 and assuming normal hydrogen thermal conductivity is shown in Figure 18. Although assuming normal hydrogen thermal conductivity is only completely appropriate when zero or equilibrium dissociation is also assumed, for this work the assumption of a constant thermal conductivity with dissociation was made. First, it was made to simplify the computa- tion, however, it can be readily justified on other grounds (see Appendix C). Relaxing the initial assumption, of constant coolant temperature, by including the energy balance equation for the coolant in the model: ch/dz - [-UO As (TC - T) + 110 A0 (To - Tc)] / Fc CPC (20) Shannan. «mom ofinmuom l upwaonoMIdoz 0595 3303.589 vozovo: .3 9»qu A53 dawned ~30qu 91 mm on on om ma 3 m o FLLhLPLPIPPrFHLherLIPIFFI—IhPFHL—Flplblr—Lph cam o o Tog m E o m o r m lccc T 18 noon: coma Jena 0 oodnsz 1... 1 4: 63% ()1) emquzodtna; 92 showed that the original assumption was largely justifiable. At the flow rates used in the experiments, the total heat content of the gas would not raise the coolant temperature more than around 5 kelvin. The second term in equation (20), above, allows for natural convective and radia- tive losses to the environment, with coefficients from Perry and Chilton [1973], p. 10-11. Inclusion of these effects has a larger effect on the coolant temperature than does the heat transfer from the gas, cooling the coolant by some 10 kelvin over the length of the reactor. The environment of the experimental system was fairly complex, and there was a fan blowing air across the outside of the cooling jacket between the cavity and the Vvacuum tank. The presence of this forced convective cooling was found to make a difference of about 5 kelvin in the outlet temperature of the cooling air. However, this additional cooling could not be modeled, but a 10 kelvin temperature drop for the coolant is probably conservative. All in all, although the measured coolant tem- peratures did not reach 400 K, assuming the coolant temperature to be a constant at this value is probably not an unreasonable simplification of. the model. Even better is that used in the sample runs for co-current flow situations, with this temperature being the initial, cavity exit, temperature for the coolant. The maximum coolant temperature possible in the experiments was 420 K, which assumes that all the 600 W of microwave power absorbed by the cavity in a typical experiment was taken off by the coolant flow normally used. As Figures 17 and 18 show, the good fit values for heat-transfer coefficients are much less than those calculated based on the properties of the gas. Thus, the simple, no reaction, heat exchanger model for the experimental system is inadequate, and reaction must be considered. 93 Assuming a ten per cent initial dissociation, and gas-phase reaction only, in a cooled reactor, as per the modelled run, the model produced the temperature profile of Figure 19. The agreement with experiment is still poor, though improving, and trial and error on the initial condi- tions ought to produce even better fits. Before this was done, though, further testing of the assumptions in the model was done, to judge whether the further trial and error would be meaningful. Surface Reaction The next assumption to be checked was that of a gas-phase reaction only. The model was modified to include a second-order, reversible, surface recombination reaction, with rate constants after Wood and Wise [1962] (see Appendix C). The mass balance equation was modified to include surface recombination as follows: dxv /dz - A (k, [H] + k, [32]) <[H21 - K [312) / F (21) de/dz - As kw ([112]w - Kw [1113) / F (22) dX/dz - dXv/dz + de/dz (23) At first, it was assumed that the energy released in surface recombina- Temperature (K) 94 — Gas -- Coolant o Expt. 1 2 3 4 5 6 Reactor Length (cm) Figure 19. Modelled Temperature Profiles for Modelled Run Without Surface Recombination 95 tion would remain in the surface, and be taken up by the coolant, making the energy balance equations: dT/dz — {mo u: / gC [dX/xz / (l + X) - dA/dz / A] + DHR dXv/dz + 2 U0 AS (Tc - T) / F) / (m0 “b / 8c T + GP) (24) [- U0 A§ (TC - T) + Uo A0 (To - TC) - F DHRw dX /dz] dTC/dz - w (25) Fc CPc Figure 20 shows the temperature profile for the sample run, under these assumptions, integrating equations (16), (21), (22), (23), (24), and (25), and it is virtually identical to that in Figure 18, which assumes no dissociation. The lower curve in Figure 21 shows, to a greatly enlarged axial scale, the simulated extent of dissociation under the assumptions of Figure 20. The upper curve in Figure 21 shows the same profile without surface recombination, replacing equations (21), (22), (23), and (24) with (13) and (20). Under the assumptions used here, the model predicts that surface recombination should consume nearly all the. free atoms in the first tenth of a millimeter of the reactor. Under the assumption of total transfer of surface recombination energy to the coolant, little recombination energy remains to heat the gas. The lowest curve in Figure 22 is the initial section of the tempe- rature profile from Figure 20, expanded to the same axial scale as Figure 21. The top curve in Figure 22 is the temperature profile for a similar run, under the assumption that the surface recombination energy is taken up by the gas, then lost to the coolant by convection ~- integrating equations (l6), (19), (20), (21), (22), and (23). The middle curve is an intermediate case, in which the surface recombination energy Temperature (K) 96 — X0 = 0.1 0 X0 = 0.0 o Expt. I I r I I I T Reactor Length (cm) Figure 20. Modelled Tem erature Profile Reactive Lossy urface 97 00326 £85 do obsooom homo; 3:3:— 823835 823°: .E 955 3:5 damned acuooom 0.“ ad ad v.0 06 md To 0.9 Nd ad cd - — b b - — P — p F p — P — . — p — _ l°°.° 25.—coca II J fins - .. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII lead pezoioossia 110110915 98 I.“ II... a.“ l | 'l '- ‘l 'I ‘I‘ ' "I '9' |' ‘l | II' ‘|'I ‘l“ | " lll. " |‘ ||' |' I. " Ill'l " udSoou \ saw 3. ........................... OGU u..- uuuuuuuuuuuuuu s I x \\1.\.\\\\\!\ coca. 99 is transferred to both the gas and the coolant by convection. This was done by estimating a wall temperature Tw such that: F DHRw de/dz - {hi (Tw - T) + hc (Tw - Tc)} AS (26) and rewriting the energy balance equations as follows: dT/dz - {1110 u: / gC [dX/dz / (l + X) - dA/dz / A] + DH dXv/dz + R hi As (Tw - T) / F) / (m0 ufi / gc T + GP) (27) dTC/dz - [-hc AS (Tc - Tw) + Uo Ao (To - Tc)] / Fc CPC (28) The middle curve in Figure 22 is, therefore, the result of integrating equations (16), (21), (22), (23), (27), and (28), such that equation (26) is satisfied. Although both the top and the middle curves in the Figure show a sharp peak in gas temperature, to above 2000 K, after the peak, the gas still cools far more quickly than the experiments reported. So far, the surface reaction rate constants had been computed based on the bulk temperature of the gas. However, they should properly be computed based on the temperature of the surface. Since the wall tem- perature computed to satisfy equation (26) was usually much higher than the bulk temperature, when there was surface reaction, there should be a significant difference, by reducing the surface, and hence the total, recombination rate. A trial and error approach (see Appendix A) was implemented to solve equation (26) for the surface temperature at which the rate of energy generation balanced the rate of energy removal. Figure 23 shows the temperature profile under this "hot wall" assump- tion, for the gas, lower curve, and the wall, upper curve, to the scale of Figure 22. Figure 24 shows the gas temperature profile from the same 100 DA “3395—34. :3» uom «058m ouaaouonaoa. vozovoz .mm 9533 E25 5934 uoaooom ad ad to ad ad so a... No no cd new In a; .. .. .. T T83 r T 183 a ,,,,-,-, s 1: ))))) :::: r :::::::::::::: 1 .. .. 18mm uuuuuuuu m. ()1) eanoJedmeJ, Temperature (K) 101 2000 1800-1 1300-] 1400-4 1200 A 1000 -4 800 4 600 —-] O 400 -4 ‘ o Expt. --— Gaa - Model 200 r T r T T T T r frr I 0 1 2 3 4 5 Reactor Length (cm) Figure 24. Modelled Temperature Profile Hot Wall Assumption r I 6 7 102 run to the scale of Figure 20. The peak in the gas temperature is moved downstream, but not sufficiently far to keep the gas from being cooled very quickly to the coolant temperature. A further, similar consideration, is whether the surface reaction is diffusion controlled, under the experimental conditions. That is, it was to be considered whether the rate of recombination at the surface might be limited by the rate of transport of atoms to the surface. This was checked by calculating the value for the concentration of atoms at the wall such that the rate of convective mass transport of atoms to the wall would equal the rate of consumption of atoms in surface recombina— tion, again, by trial and error (see Appendix A). The equation that had to be satisfied was as follows: hIn As ( [H] - [H]w) / [1 + Xw/(l + Xw)] - F de/dz (29) The convective mass transfer coefficient was computed using a Sherwood number, analogous to the Nusselt number in heat transfer, and so, by a conventional analogy between heat, mass, and momentum transport (Bird, et a1. [1960], p. 645) also assumed constant at 3.66. The diffusion coefficient and density needed for the Sherwood number computation were computed under bulk conditions. The gas temperature profile derived from integrating equations (16), (21), (22), (23), (27), and (28), such that equations (26) and (29) are satisfied, for the modelled run is shown, to the scale of Figure 21, in Figure 25, and, with the wall and coolant temperature profiles, to the scale of Figure 20, in Figure 26. The peak is shifted yet further downstream, very slightly, as the surface con- centration, corresponding to 0.8 per cent dissociation against the bulk 10 per cent, is lowered enough to have an impact, however slight, on the 103 CA ad vozohnoo down: an £05955: Sch 3595 9:595 Bop. nozouoi .mm on: E 0.9 L F 325 among uoaooom b6 96 _ . h b 06 r F Tc _ b om p unseen—:82 ..i flOwnfiuzn— ll ()1) eanozedmel Temperature (K) 104 2400 — Bulk Gas Temperature 4 --- Surface Temperature i —— Coolant Temperature 2200 a] 0" IF 0 O l 600 -] 0 400 _],___ ...... :.'.:.=.r.- 200 # T I I f T T T I T ‘I T ‘r O 1 2 3 4 5 6 Reactor Length (cm) Fi ure 26. Modelled Tem erature Profiles for Mode ed Run - Pltg‘gli-lflow odel - One Dimensional Y Developed 105 surface reaction rate. The surface reaction consumes the free atoms very quickly, even so, as shown by Figure 27, which compares the profiles of dissociated fraction for volume recombination only (upper curve) and for volume with surface recombination (lower curve), over the first 7 cm of the recombination zone. Figure 28 shows the gas and coolant temperature profiles for the modelled run, but with counter-current flow geometry. This was simulated by making the coolant flow rate negative when it was input to the program. The initial coolant temperature was assumed to be 305 K, to allow for the expected heat load and the roughly 300 K inlet air tem- perature, at the far end of the recombination zone. The profile is nearly identical to that of Figure 26, the main difference arising from the different coolant, and hence ultimate, temperatures. The simulation agrees fairly well with the experiment, in this case, predicting that the gas temperature should reach equilibrium with the coolant very quickly. The effect of the mass flow rate through the system is shown in Figure 29, which shows gas temperature profiles for runs at the lowest hydrogen flow rate, and the corresponding plasma pressure, used in the axial dependence studies, .055 mg/s and 4.26 Torr (short dashed curve), and at twice the modelled run's flow rate and pressure, 1.54 mg/s and 52 Torr (long dashed curve), alongside the modelled run (solid curve). The Figure shows that the temperature profile becomes independent of the mass flow rate, beyond 5 cm below the plasma, over the range of flows studied in the experiments, which agrees with the experimental results. However, the model predicts that this is due to the gas temperature Fraction Dissociated 106 0.10-1x —- Active Surface -- Inert Surface 0.00 Reactor Length (cm) ure 27. Modelled Dissociation Profiles Effect o Surface Reaction — Fully Developed Model 107 1600 — Gas j - - Coolant n t. 1400 -4 Exp 1200 -* 1000 -J 800 -— Temperature (K) 600 -+ 400 e 200 .T. # T T I T Reactor Length (cm) Figure 28. Modelled Temperature Profile Counter—Current Coolant Flow 44L.-.J._--_.-_-_. Temperature (K) 108 --- High Flow — Standard —- Low Flow o Expt. I r fl I I I l O 1 2 3 4 5 Reactor Length (cm) Figure 29. Modelled Temperature Profiles Flow Rate Dependence 109 coming into equilibrium with the coolant, and approaching the coolant temperature. Since the coolant, even if it absorbed all the incident microwave energy, could not exceed 420 K, the model does not explain the 600 K reported gas-dynamic temperature. Figure 30 shows the simulated gas temperature profiles assuming initial dissociations of 0., .l, .2, and .3, with other parameters as in the modelled run. The higher the initial dissociation, the higher and, very slightly, the broader the initial peak in the temperature profile. The peak still disappears by 5 cm, however. A higher initial dissocia- tion is probably not realistic, and would not, probably, improve matters enough. Figure 31 shows the effect of raising the initial temperature to 1500 and to 1800 K, from the modelled run's 1200 K, holding all other parameters as in the modelled run. As with increasing the initial dis- sociation, the increased initial temperature raises the peak and broad- ens it slightly. The net initial rise actually decreases with increasing temperature, however, and the peak disappears within 4 cm of reactor length, anyway. Figure 32 shows the gas temperature profiles for a non-dissociated and 0.1-dissociation run, with the heat-transfer coefficients multiplied by 0.1 from the literature values, and all other parameters as in the modelled run. Comparing these curves with the two lowest curves in Figure 30 shows a dramatic difference. The cooling is much slower, the peak temperature is higher, and the peak is broader. It seems to suggest that the correlations used to estimate the heat-transfer coefficients have given values that are too high, compared to the real values, since these profiles are well above the experimental profiles, and adjusting the initial conditions and the heat-transfer coefficients may serve to 110 Cc 0.553an: 7. 3.210. a 0000 _.______ F6 0000 XXXX T -—. _.n_ 1.5 r 14 T 1.3 [2 [1 \‘I lllllllllllllllllllll I 4J+4+.44HJlmuJH-j_.n A _ n O O O 0 O O O 0 O O O O O O O O O 0 O O O 0 0 2 O B 6 4 2 O 8 6 4 2 2 2 1. 1 1. 1. 1 Reactor Length (cm) Figure 30. Modelled TemBerature Profiles ependence Initial Dissociation Temperature (K) lll he 3% O O l r0 - 1000 d -- 'ro - 1500 — 1'0 - 1200 i -——. m.~ .o..<.... -‘ I Reactor Length (cm) (II—1 G-‘i 4.4... . -_ Figure 31. Modelled Tem erature Profiles Initial Temperature ependence Temperature (K) 112 2400 2200—J 2000 -4 1800 4 I \ 1600 _ [I \\ pa 3 O _L / [I l -- X0 ==OJl ——- XD‘= 01) 200 “—1— T' r l r r ‘r T 'Ti r l j 0 l 2 3 4 5 6 Reactor Length (cm) I Figure 32. Modelled Temperature Profiles - Conductivity 0.1 Hydrogen Initial Dissociation Effect 113 bring the simulation into line with the experimental results. Yet, the adjustment on the heat-transfer coefficients is not reasonable. Engi- neering correlations, as were used in the simulations, tend to be con- servative. In heat transfer work, conservative values for heat-transfer coefficients tend to be underestimates, to ensure the heat-transfer apparatus will continue to perform its designed task, even if it is somewhat fouled. Another use of this adjustment is to investigate, to a first ap- proximation, the behavior of other gases. If the product of the flow rate and the heat capacity of the second gas is the same as the product of the input flow rate to the simulation and the heat capacity of hydrogen, and the heat—transfer coefficients multiplied by the ratio of the thermal conductivity of the second gas to that of hydrogen, the output of the model will approximate the performance of the second gas. This assumes that the same constant Nusselt number applies, making the internal heat-transfer coefficient solely dependent on thermal conduc- tivity. This also involves neglecting the reaction, or requiring zero dissociation, and neglecting the differing dynamics and thermodynamics of compressible fluid flow for the different gases. In argon, the ther- mal conductivity is about one-tenth that of hydrogen (Perry and Chilton [1973], p. 3-215). The argon flow rate used in the experiments whose results are plotted in Figure 13 corresponds to the hydrogen flow of the modelled run, when corrected for the differing heat capacities. Since argon is a monatomic gas, the zero-dissociation curve in Figure 32 compares rather well with Figure 13. The thermal conductivity of nitrogen is about 0.15 that of hydrogen (Ibid.) over the range 300 to 1200 K. Simulating the experimental runs 114 in nitrogen by the modelled run at zero dissociation and with the heat- transfer coefficients reduced by this factor produced the gas tempera- ture profile shown by the upper curve in Figure 33. Also shown on the Figure are the experimental gas-dynamic temperatures determined in nitrogen flows, with co-current coolant flow. The agreement between the computed profile and the measured one is remarkable. However, this was done for a flow rate corresponding to the highest nitrogen flow used in the experiments. Reduce the flow rate by a factor of five, to match the lowest experimental flow, and the model produces the temperature profile shown by the lower curve in the Figure. Since the experimental tempera- tures were very close together over the whole range of flow rates, as shown by the cluster of experimental points, the model obviously does not accurately reproduce this lack of dependence on flow rate. The counter-current coolant flow experiments, as shown in Figure 10, show a much larger dependence on flow rate, much more as a simple heat-ex- changer model, such as this model at zero dissociation, would predict. Only the point at 6.4 cm below the cavity showed any residual heating in counter-current flow. Two-Dimensional Model Since this modelling has shown that consideration of surface recom- bination necessitates sharp radial gradients in temperature and con- centration, under the conditions of the experimental system, a two- dimensional model was developed to further investigate the system and, it was hoped, to model the results better. The Peclet numbers, relating the bulk flows of mass and energy to the diffusion and conduction, 115 mufiofiwuoaxm cowouuaz .m> nowoup%: ma. I hua>uuoapcoo .0>Huomom-coz - oHHmoum ouauouonaoe noHHovoz .mm madman A88 finned youoeom on on mm on ma o“ #17. _L1__ptlu|L[L...l—L[Pb.bft—FLLC F40 P Jana o Form #3 .. I to: sea .I ()1) emissadmej, 116 respectively, are 0.88 and 1.2, also respectively, using the transport properties as calculated at 1200 K and the tube diameter as the defining dimension. This further suggests that radial diffusion and conduction should be important. The energy and material balances for a two- dimensional model are as follows: First, defining dX/dzR - (k1 [H] + k2 [H21> ([H21 - K [H12> / G (30> dX/dz - dX/dzR + D d(r d[H]/dr)/dr / G(r) (31) H-H2 dT/dz - [DHR dX/dzR + A d(r dT/dr)/dr / G(r)] / CP (32) with boundary conditions in the radial direction: from symmetry: dX/dr|r._0 - 0 (33) dT/dr|r_O - 0 (34) from mass and energy conservation in wall recombination: d[HJ/drw - kw ([HZJW - Kw [H12> [1 + xw/<1 + xw>1 / D (35) H-H2 A dT/drw - DHR kw ([112]w - Kw [012) / As - q (36) where q - - hc (TC - Tw) (37) Additional terms, as in equation (15), were included in the energy balance equation (32) to account somewhat for inviscid flow effects. These inviscid flow terms were based on the bulk velocity, computed from the mass-average temperature and composition. Similarly, equation (16) was used to give the pressure dependence, assumed constant across the 117 reactor, based likewise on the mass-average state. That is, the bulk temperature was computed assuming the heat capacity per unit mass was independent of the variation of composition or temperature across the tube. The mass flux down the tube was assumed constant with length, and parabolic in cross-section, that of an isothermal laminar flow. The radial diffusion term was modified to ensure this, transforming equation (31) to: dX/dz - dX/dzR + D d(r d[H]/dr)/dr / G(r) / [1 - X/(l + X)] (38) H-H2 This provides for the equal-mass counterflow, also apparent in equation (36), the surface boundary condition providing for conservation of mass in the surface reaction. These. equations were solved by the method of lines, that is, each of the transverse derivatives was written in its finite-difference form, rendering the axial derivative at each station across the tube as a function of the conditions at the stations on either side, as well as those at the station itself. To help ensure the conservation laws were followed, the transport properties used in the transverse derivatives were computed at the arithmetic average of the temperatures of the two stations involved. In the finite-difference form used, the mass flux through each radial segment was taken to be the average value based on the parabolic profile, and the conditions were assumed to be uniformly that of the station within each segment, centered except for the outer- most one, which was taken to be at the wall. The wall boundary condi- tions, equations (35) through (37), were not explicitly invoked in solving the equations. Instead, the outermost section was assumed to be at the wall conditions, and that section's material and energy balances 118 included the wall recombination reaction and losses through the wall, along with flow through that section. The first station was located on the center line of the reactor. Experimentation with the number of stations used showed that there was little difference in the radial profiles when five or nine stations were specified, under the standard set of parameters. Hence, for com- putational economy, five stations were used in the model computations. This model produced results much like the earlier ones (Figure 34). Though the bulk, or mass-average, temperature profile is a little closer to the experiment, the improvement is not significant. This figure also shows a temperature profile that makes allowance for the assumptions involved in the gas-dynamic temperature measurement (dotted curve) by allowing, to a first approximation, for the change in molecular weight and heat capacity due to dissociation. This curve shows roughly the gas- dynamic temperature profile that would be estimated under the assump— tions made in the experimental program if the model accurately reflected the experimental system. This predicted gas-dynamic temperature Tg(e) is computed from the bulk temperature T by Tg(e) - T (7 + 3 X) (l + X) / (5 + X) / 1.4 (39) The closeness of the predicted gas-dynamic temperature curve to the bulk temperature partially validates using the experimental gas-dynamic temperatures as a reasonable estimate for the gas temperature, ignoring the inaccuracy implicit in the assumptions used in computing the ex- perimental values. Figure 35 shows the bulk temperature profile at the lowest flow from the experimental series, and the model continues to Temperature (K) 119 7 T" I To s Tue)! oExpt.§ i i i 3 I l ; °| i 3 l i zooJi‘fiFrrrr'rTrvrri 0 1 2 3 4 5 6 7 . Reactor Length (cm) Figure 34. Modelled Temperature Profiles for Modelled Run Two-Dimensional Model - Five Stations Temperature (K) 120 1600 - _, —— Model . i o Expt. 3 1400 i 1200 --4 1000—1 800 '1 soul 0 a 400 -J A: 200 i I fl r f T 'r r r T“ I T 0 1 2 3 4 5 6 Reactor Length (cm) Figure 35. Modelled Average Temperature Profile - Low Flow Two-Dimensional Model - Five Stations - Cool Sheath 121 predict a much higher dependence on flow rate than the experiment showed. Axial Dispersion Model The next modification to the model that was considered was the inclusion of axial dispersion. Since the first data point in the experi- mental system was located about three tube diameters downstream from the plasma, the Peclet numbers using this distance as the dimension were about 3. This supports the idea that axial dispersion should be consid- ered, as diffusional transport is thus a significant fraction of the bulk flow. Adding axial thermal conduction (through the gas) to the two- dimensional model changes the energy balance equation (32) to: CP dT/dz - A d2T/dz2 / G(r) - [DHR dX/dzR + A d(r dT/dr)/dr / G(r)] (40) and adds an additional boundary condition: dT/dz|z_co - o (41) Analogous modification to the material balance equation (21) allows consideration of axial diffusion, as well. Solution of these equations by reducing the second-order differential equations to first-order equations by using the first axial derivatives as dependent variables and using a Runge-Kutta marching technique (see Appendix A) to solve the coupled system proved to yield unstable solutions (See also Coste et a1. [1961]). Coste's solution technique was adopted to find stable solu- tions. The non-dispersion equations were solved, first, by using an implicit cell model, which, as Coste shows, is roughly equivalent to a dispersion model if the cell length is chosen appropriately, based on 122 the relative magnitude of velocity and dispersion coefficient. The cell length used was computed based on the thermal conductivity at bulk (mass-average) temperature and bulk flow velocity, computed at (mass-) average temperature and concentration. Once the end of the reactor was reached, the reduced, first-order, dispersion equations were integrated backwards, using the Runge-Kutta techniques on an arc-length method to reduce the problems from stiff equations. In addition, the program was modified to explicitly solve the boundary conditions at the walls (Equa- tions (23)-(25)), by the trial and error method described in Appendix A, in hopes of reducing the instabilities further. This changed the dis- tribution of stations within the reactor, so that the first station was effectively half a cell-width off the center line of the reactor, and the last station an equal distance from the wall. The two-dimensional model, in this version, could be used to fairly mimic the most highly developed one-dimensional model by using one radial station. In effect, this use would be equivalent to a one-dimensional, axial-dispersion model with surface reaction and heat and mass transfer to the walls using constant Nusselt and Sherwood numbers of 4.0. This is not out of line with actual Nusselt and Sherwood number values for well-developed laminar flow (Perry and Chilton, [1973], pp. 10-12,10-13). Comparing Figures 34 and 35 with Figures 26 and 29, respecively, shows that the two-dimensional model actually made little difference in the profile in the plug-flow case, as well. Runs with one station were therefore used to test-run the model and conserve computation time. On the forward, implicit, cell-model integration, this technique generated temperature profiles little different from the non-dispersed model. The temperature was somewhat higher around 5 cm. downstream, 123 though not significantly, and the initial peak had been passed at the first cell. However, this technique, on the backwards integration, still produced instability, of the sort reported by Coste, g§_§l. [1961]. Coste's method probably only allows for stable integration of equations within a certain range of parameters, such as the relations among flow rate, dispersion coefficient, and reaction rate. It seems likely that the system considered here does not fit in this range, and so cannot be integrated stably either backwards or forwards. Therefore, a new in- tegration technique was adopted, an implicit version of the Runge-Kutta method, as described by Cash [1975] (see Appendix A). This technique proved to consume excessive amounts of computer time, as developed from scratch, so an implicit-method integrator routine from the IMSL library, their routine DGEAR, was selected and the program rewritten to access the integrator. At first, the initial derivatives of temperature and concentration were assumed to be those computed without consideration of axial disper— sion. This resulted in instabilities, including instances of predicted negative concentrations of hydrogen atoms, and temperatures oscillating up and down along the tube. Assuming zero initial derivatives produced better results, though with its own inconsistencies. The program has never modeled more than the first five mm of the reactor under these assumptions, due to some uncorrectable difficulty in solving the surface equations. Beyond a certain point, perhaps in a certain range of near- surface conditions, possibly high dissociations at low temperatures, the solution techniques seem to fail to find a solution for the boundary conditions, and the program runs indefinitely. It is unclear what the conditions that cause this are, since it seems that the program had 124 little problem solving the equations under very similar conditions. Various techniques for solving the equations describing the boundary conditions have been attempted, with various degrees of success but no breaking through to model the full length of the reactor. Even these short runs, though, produce some results that might be significant in considering the value of the models. In the pseudo-one- dimensional case, using a single station, the bulk dissociation rises briefly, then drops off slowly, much more slowly than in the non-dis- perse model. The bulk temperature drops off as sharply as in a non- disperse case, reaching 400 to 500 K in the first few mm below the atom source. However, the surface temperature and dissociation both rise sharply, to around 3000 K and over 0.5 dissociation. The surface concen- tration, however, still remains lower than the bulk, permitting mass transfer to the wall, despite the higher degree of dissociation, due to the increase in surface temperature. It is rather startling that the bulk gas temperature in the model drops as the surface temperature rises. Apparently the axial dispersion of energy is computed to be so much greater than the radial dispersion that radial conduction effects haven't shown up as far as the model has been able to go. Using multiple radial stations has produced much the same results, and, in addition, has predicted that radial gradients seem to persist and even intensify down the pipe. Differences between two stations seem to increase with distance, though radial transfer should damp out such differences. Besides the added computer time to compute the additional stations, the problem of the surface equations locking up seems to arise even faster when using a number of radial stations, so this is based on even shorter reactor sections than the quasi-one-dimensional results are. 125 One possible source for these results, counter-intuitive to the point of being non-physical, is the model itself, and its technique for handling the second-order equations. To solve a second-order differen- tial equation, such as (40) above, the model rearranges the equation to the form of: A d(dT/dz)dz / G(r) - cP dT/dz - [DHR dX/dzR + A d(r dT/dr)/dr / G(r)] (42) Then, the coupled system of equations, with dT/dz as a separate depend- ent variable, is solved by using first-order methods. If it is assumed that the last two terms on the right side of equation (42) can be neglected, it reduces to‘ A d(dT/dz)dz / G(r) - CP dT/dz (43) which can be integrated analytically to dT/dz - dT/dzlz.O exp(G(r) CP 2 / A) (44) T - TIz—O + dT/dzlz_o A / G(r) / CP [exp(G(r) CP 2 / A) - l] (45) Equations (44) and (45) tend to blow up with increasing 2, tending to infinity, either positive or negative depending on the initial deriva- tive. Hence, it is to be hoped that the terms in equation (42) that were neglected should not be. Analysis assuming that those terms are constant only adds a linear term to the final equation for T, which will be swamped by the exponential term, given sufficient distance down the reactor. In addition, the terms enter the equation in a manner that seems contrary to physical expectation. It seems unreasonable that the generation term, the second term on the right side, should have a nega- tive influence on the second derivative, and through it, on the first. Likewise, the radial conduction term has an influence opposite to that 126 it would have in the non-disperse case. In the program, the subroutine that computes the dispersion terms does so by calling the subroutine used to compute the first derivatives for a non-disperse model, and subtracting the outputs from the subroutine from the dependent variables that are used as the first derivatives, and applying the appropriate multipliers, of flow rate, capacity, and dispersion constant. Perhaps this subtraction should be reversed, so that the system won't blow up, and will yield more physically reasonable results, yet repeated analysis of the system yields the model as given. It is unclear what the answer to this problem is. This might be one reason that Coste, g§_gl., [1961] said that forward marching methods are unusable for this sort of second- order differential equation problems. The problems with the implementation of the Coste method program were eventually traced to incorrect initial conditions on the reverse integration. A more generalized program, using an arbitrary step size and a search method to find the downstream conditions to initiate the reverse integration, was developed, and it proved to yield results comparable to the Gear method. However, the reverse integration method required much more computer time, so the Gear method integrator was selected for further work. It was eventually realized that randomly selected initial condi- tions would not necessarily yield results that would converge on the infinity boundary conditions (Equation 41). Indeed, a marching integra- tion on the second-order equation (42) will produce, for example, values of temperature that will be driven away from the neighboring tempera- tures, downward if the temperature is less than the neighbors, upward if greater. Therefore it seemed reasonable and experimentation verified 127 that properly selected initial values for the derivatives may generate results that are realistic and reasonable. Since the IMSL DGEAR inte- grator proved as effective and efficient as any, it was selected and a search algorithm designed to find the appropriate initial derivatives (see Appendix A). However, the search technique failed to find a set of initial derivatives, for the standard model, that would allow the integration to proceed beyond a short reactor length without yielding results that were assumed to be incompatible with the boundary conditions at infinity. Inconsistent results, to the point of having two runs with the same set of initial derivatives invoking contradictory indicators, and a rela- tively strong dependence of the point towards which the search would converge on the value of the integration and surface boundary conditions convergence parameter indicated a possible numerical instability in the integration. A possible way around this would be to search for initial conditions that would allow integration to some preset length, back up to some point where the variation in the search is expected to be small, and search for the derivatives at that point that would permit integra- tion to proceed further. This might cause a slight discontinuity in the derivatives at the transition points, however. Investigation of this technique is left for future work. Before the instabilities appeared, though, the simulation reached a point at which most of the changes within the modelled reactor seemed to have completed, and the search had reached the point where this regime had been fairly well established, with differences from one run to the next being smaller than the precision of the output. The movement of the search convergence point with the convergence parameter was within this 128 range. The temperature profile from the longest run, standard condi— tions, quasi-one-dimensional, using the diffusion coefficient and con- ductivity for the values of the mass and thermal dispersion coefficients respectively, is shown in Figure 36 (solid curve), along with the same profile from the plug-flow model (dashes). Comparing the curves shows that the temperature peak when dispersion is considered is much lower than that in the plug-flow situation, and the peak is much broader. The profile of fraction dissociated is shown in Figure 37, again paired with the plug-flow profile, and the decay of the free atoms is much slower. However, the temperature peak is largely dissipated at 5 cm. below the reactor. The model is approaching the results of the experiment, but has not yet reached them. A fully two-dimensional dispersion model might provide results more in agreement with experiment, and the program developed in this work could be used provided sufficient computation time is available. The search for the quasi-one—dimensional case involved large amounts, on the order of hours, of computer time, and a n-station case requires a 2n-dimensional search. Or else, the initial conditions might be adjusted to fit the experiment better. Figure 38 shows the temperature profiles calculated at the low flow rate and pressure (solid line), and at the standard flow rate (long dashes), from the dispersion model, and the low flow rate from the plug- flow model (short dashes). Inclusion of dispersion has a much greater effect at the lower flow rate than at the higher, in broadening the peak, and the dependence on flow rate is much lower in the dispersed than in the plug-flow. Thus, the lack of flow-rate dependence seen in the experiments is seen as an indication that axial dispersion needs to be considered. The modelled results do not yet match the experiment in Temperature (K) 129 1600 V 1 4' ‘A 1“. 1400—l: \, u , l l i o Expt. ; ‘ -— Axial Dispersion Model 'e 200 -- Plug—Flow Model J "r—‘v—‘I‘VI'IYTttIvT I 0 1 z 3 4 5 6 7 Reactor Length (cm) Figure 36. Modelled Temperature Profiles Axial Dispersion vs. Plug—Flow Fraction Dissociated 130 0.10 P o 0: l — Axial Dispersion Model -- Plug-Flow Model Reactor Length (cm) Figure 87. Modelled Dissociation Profiles Axial Dispersion vs Plug-Flow ad’s. on. ~—_.- ----'-—.~--.-. . Temperature (K) 131 an O O L -— Low Flow Rate. Dispersed -— Standard Flow Rate. Dispersed --- Low Flow Rate. Plug Flow 0 Expt. Reactor Length (cm) ‘ ure 38. Modelled Temperature Profiles ispersion Model — Flow Rate Effect 132 every detail, but further adjustment of initial conditions and the dispersion coefficients may provide much better fits. Chapter 7 - CONCLUSIONS A spacecraft propulsion and maneuvering system has been proposed that would use electrical energy, from a solar cell or other source, to heat the working fluid, or propellant, through a microwave-frequency plasma. The heated gas would then flow through a nozzle to generate thrust, much as in a conventional rocket. The main difference is that a chemical rocket heats the propellant with energy liberated by chemical reactions within it, whereas this system heats the propellant with energy obtained from outside it. A number of engineering problems remain to be solved before this microwave-plasma electrothermal rocket becomes practical. Among these problems is the basic one of how efficiently the energy is transferred, or coupled from the microwave-frequency electro- magnetic fields to the gas. More specifically, it is desired to generate translational kinetic energy of the molecules of the propellant. This motion, random at first, can be turned to directed flow and thrust by a nozzle. One potential mechanism for this coupling is the dissociation- recombination cycle of a diatomic gas in the working fluid. In the plasma, the diatomic gas would dissociate, to recombine and release the energy absorbed in the dissociation downstream of the plasma. Hydrogen has a number of advantages as a propellant, particularly if this cycle is important, and this work examines several aspects of the use of hydrogen in this proposed system. The work reported here is divided into three distinct sections. First, the paper describes a program of preliminary experimental work, intended to demonstrate the capability of the microwave-plasma to couple energy to the gas from the electromagnetic fields, to measure the thrust 133 134 and temperature profiles generated in the gas flowing through a micro- wave-plasma, and to test various experimental designs for their effi- ciency in coupling energy to the plasma. A modelling study was under- taken to evaluate the importance of vibrational and rotational excita— tion on the hydrogen recombination reaction and energy transfer in the reaction, and its effects on a practical microwave-plasma electrothermal rocket. A series of computer models have been developed to explain the experimental results and to provide a basis for the design of future experiments and practical microwave-plasma electrothermal propulsion systems. The three sections of this research each lead to their own independent conclusions, though there are also certain common conclu- sions from comparison of the results from different phases of the work. The experimental work mainly involved relatively low-pressure, up to 100 Torr, plasmas in hydrogen, with flows to 2 mg-mol/s, and incident microwave powers between 300 and 600 W. Work was also done, for compar- ison and safety considerations, in nitrogen and argon, under similar conditions. The experiments showed that the microwave-plasma system can pump energy into the flowing gas. Thrusts were measured directly, using a vane-type thrust stand, and temperature profiles obtained, through a gas-dynamic technique, that were in rough proportion. That is to say, the thrusts obtained, at the same flow rate, were approximately propor- tional to the square root of the estimated temperature, in configura- tions of the apparatus in which the thrust measurements were most reli- able. With the plasma, when the gas-dynamic temperature estimate was 400 K, the thrust was approximately 17 per cent higher than the thrust measured by flowing the same rate of cold, room temperature, roughly 135 300 K, gas through the apparatus. A predicted mode, TE012, of the microwave resonant cavity was verified, for use at higher plasma pressures. However, the experimental arrangement pumped the absorbed microwave energy out of the gas again as fast. The experiments' poor design forced this loss of most (60 to 80 per cent, and more) of the incident micro- wave energy to the coolant, as the gas could not carry off the "over- kill" of microwave energy by any available mechanism. In addition, the cooling needed to protect the plasma chamber walls also cooled the region below the plasma, so the gas retained but a small fraction of the energy fed into it in the plasma at any distance downstream. If the temperature was estimated far enough downstream, though still within the apparatus, the gas and the coolant were at the same temperature. With the air coolant and the plasma gas flowing counter-currently, in hydro- gen, the gas-dynamic temperatures were 300 K as close as 6 cm below the plasma. With the coolant flow reversed, that temperature had dropped to 600 K, from 1200 and more at the plasma exit, at 6 cm and had cooled to match the coolant, 400 K, by the end of the plasma tube, 36 cm below the plasma. Accurate thrust measurements required that the nozzle be at the end of the plasma tube, for free expansion into the vacuum tank housing the thrust stand, and some 36 cm from the plasma. This was so far down- stream from the plasma that the gas and the coolant were nearly equili- brated. Though the values of the measured thrust were not in agreement with theoretical values, as computed from gas-dynamic theory, the measured values being from 25 to 110 per cent of the ideal, depending largely on the pressure around the thrust stand, their relative values 136 were in rough accord with the temperature estimations, as described above. The temperature profile, in hydrogen, especially, was very nearly independent of the flow rate. The experiments demonstrated that the cavity resonant mode was a controlling factor on the plasma behavior, that different modes could sustain a plasma in different pressure regimes. Attempts to cool the wall by transpiration cooling proved unsuccessful, due to failure of the porous wall and migration of the plasma out of the transpiration region. An attempt to create high-pressure plasmas by a capillary tube design also was unsuccessful. The experimental work was basically preliminary, but yielded a number of curious and potentially significant results for further exper- imentation. For example, there was an as yet unexplained extreme back- ground pressure dependence of the measured thrusts, that had also been reported elsewhere, though not, I believe, to the degree reported here. There are a number of ways to improve the experimental design to avoid the problem of overcooling the gas. The microwave power level must be better matched to the capacity of the gas. The cavity resonant mode and the gas flow pattern can be selected to reduce the danger to the tube walls. One potential method, transpiration flow, was tried and could not be made to work, due to failure of the porous wall and migra- tion of the plasma out of the transpiration region. The distance between the plasma and the nozzle can be reduced to reduce the downstream cool- ing, it is to be hoped, without reducing the accuracy of the thrust measurements. A possible alternative design, for both experiment and applications, is regenerative cooling, where the feed to the plasma is ducted past the nozzle to cool the nozzle and pre-heat the feed. This is 137 a common design in more conventional chemical rocketry. Another possible experimental arrangement, under discussion when this work was being done, would be to mount the entire cavity, plasma tube, and nozzle in a vacuum tank, on a thrust stand. This would help diagnose the origin of the back-pressure dependence, as well as reducing the distance for cooling between the plasma and the nozzle. A potential sink for the recombination energy was thought to be the internal, vibrational and rotational energy of the recombined molecules. If the recombined molecules retained a significant fraction of the energy of recombination in these modes for a significant length of time before the energy would be transferred to translational kinetic energy, the performance of a rocket that relies on the dissociation-recombina- tion cycle to transfer energy to the working fluid might be seriously degraded. A computer modelling study, based on a two-temperature model of vibrational excitation, was done to estimate how important this internal energy of the molecules would be in a practical system. This study showed that it is most likely that molecular internal energy considerations will not be important in practical electrothermal rocket design. At temperatures and pressures likely to be used in such designs, any internal energy excess will largely thermalize very quickly, within microseconds, compared to the millisecond time scale of the recombina- tion reaction. The energy that would be stored in excess vibrational excitation is likely to be but a small fraction of the total energy released in the recombination reaction. Under practical conditions, this loss of energy would reduce the temperature of the recombined gas on the order of a few kelvin compared to the 3000 K operating temperature. 138 However, low temperature, high pressure work with atomic hydrogen may require consideration of internal energy. Attempts to correlate the experimental results with the macroscopic computer models have verified that wall cooling must be considered in any model of the experimental arrangement as used, as the profound influence of the cooling geometry on the experimental temperature profile shows. In addition, the models show that the surface reaction is important, under the experimental conditions. This reaction creates conditions at the surface that imply large thermal and concentration gradients across the tube, hence two-dimensional models are required. Though many of the qualitative characteristics of the experimental results have been reproduced, quantitative agreement remains elusive. The absence of dependence on flow rate has been seen, but only as a result of total equilibration with the coolant. The reaction is all but complete and the temperature equilibrated by the time the model reaches the first experimental data point. However, the experiments show no flow rate dependence at a temperature too high to be explained by equilibra- tion at that first point. In addition, some experimental results in nitrogen parallel those in hydrogen almost exactly, while the models suggest that, at some of the experimental flow rates in nitrogen, equilibration should not have occurred at that first point. Further improvement and complication of the model are required, and one including axial dispersion has been developed to investigate the matter further. A solution technique involving a search for the initial derivatives is developed. Preliminary results indicate that dispersion may play a role in the performance of the experimental arrangement, 139 though instabilities limit the utility of the model. The work has indi- cated the directions further develoPment should take, coupling the two- dimensional and dispersion models together, and trial and error or theoretical studies to determine better values for the initial conditions and the dispersion coefficients. The experiments have shown that the microwave-plasma electrothermal rocket may be feasible, various diagnostic techniques have been demon- strated, and some of the constraints on an experimental arrangement have been identified. The vibrational-energy modelling study has shown that one possible constraint on the performance of a microwave-plasma elec- trothermal rocket is unlikely to be an important factor in practical design of the system. The further modelling has shown that cooling, surface reaction, and axial mixing must be considered under the condi- tions of the experiment to adequately model the experiment. It has also indicated that these factors are likely to be important in other config- urations. Flow rates and pressures must be much higher than those in the experiment before axial mixing and surface reaction can be neglected. Though much work remains before the microwave-plasma electrothermal rocket takes flight, this work, at least, has demonstrated its scien- tific feasibility. The major problem to be solved is the question of surface reaction, and keeping the energy in the gas. APPENDICES APPENDIX A Appendix A - COMPUTER PROGRAM DESCRIPTIONS All inputs, outputs in Systeme Internationale units. If 2 input less than 0., ends computer run. Inputs echoed. Inputs free format. Changes From Cyber to Vax The computer work was begun on the MSU Computer Laboratory Control Data 6500 computer, and some work was done using a NASA-Lewis IBM 360. The 6500 was replaced with a Control Data Cyber 750-175, and then, for reasons of computer access, the programs were transfered to a MSU College of Engineering Digital Equipment VAX 11/750. This produced a few changes in the code, some mandated by the new computer, others facilitated by the VAX/VHS operating system. Mandatory Changes: 1. To maintain consistency with the Cyber 64-bit word, double precison (Real*8) was required on the 32-bit VAX. 2. Code structures of the form A-B-C, permitted on the Cyber, are non-standard, and forbidden on the VAX. 3. Certain statement function structures would not work, and those statement functions had to be transformed into Function Subprograms. 140 141 Other Changes: 1. The program was redesigned to save sufficient data to permit a limited interrupt-restart capability. 2. The code was modularized, permitting more efficient debugging and modification. 3. A prompt for the initial inputs was added. The Cyber had a prompt utility that the IBM and VAX lack. In Appendix B: Program 1: Vibrational Energy Model. One-Dimensional, Plug Flow, Adiabatic, Gas-Phase Reactions. Inputs: Initial values for z, X, T, P, and T Hydrogen flow rate V! to plasma, Y, convergence value, area of reactor cross-section, nozzle throat area, position of end of reactor, printing interval. If input TV less than 0, uses steady-state value as initial value. Outputs: At initial 2 and print intervals, prints 2, X, T, P, TV, reactor flow area, bulk velocity, local sonic velocity, and Mach number. Program 2: One-Dimensional Model. One-Dimensional, Plug Flow, Cooled, Reactive Surface. Inputs: Initial values for z, X, T, P, and Tc. Mass flow rate in kg-mols molecular hydrogen fed to plasma per second, convergence control 142 parameter, the two area parameters, constant and throat areas, final reactor length, print interval, cooling air flow rate in kg-mols per second, a real-valued heat transfer-coefficient control parameter, U0, two integer-value mode control parameters, MODE and M2, and a real- valued viscous-flow control parameter, FR. Control Parameters: U0: If U0 greater than 0, internal heat transfer coefficient com- puted from conductivity and Nusselt number, then multiplied by U0, external coefficient is constant 47.06 * U0 W/m2 K. If U0 equal to 0, internal heat transfer coefficient computed as though U0 - 1.0, external coefficient 0. Adiabatic with reacting walls. If U0 less than 0, then the absolute value of U0 in W/m2 K is used for both internal and external heat transfer coefficients. MODE: Value Result <-0 Wall Temperature based on convective heat transfer alone, with surface recombination energy transferred to cooling air. 1 Non-reactive walls 2 Wall Temperature based on convective heat transfer alone, with surface recombination energy transferred to gas. 3, >4 Wall Temperature based on energy balance between convec- tion to gas and to coolant, and surface recombination rate. 4 As MODE - 3, but wall reaction rates computed from bulk temperature, rather than wall temperature. 143 M2: If M2 - 0, Dissociation at surface assumed bulk value. Other- wise, computed based on mass balance at surface, with convective mass transfer balancing surface reaction. FR: If FR - 0., inviscid flow. Otherwise, Fanning friction factor, from standard correlations, multiplied by FR to compute viscous losses. Outputs: At initial 2 and at print intervals thereafter, prints 2, X, T, P, Tc, Bulk Velocity, Flow area, Sonic Velocity, Estimated Corrected Gas Dynamic Temperature, and the Wall Temperature (at downstream print stations only). Program 3: Two-Dimensional Model. Two-Dimensional, Plug Flow, Cooled, Reactive Surface. Inputs: Initial Values for z, , , P, Tc’ Mass flow rate, convergence limit, area parameters, reactor length, print interval, cooling air flow, real valued heat transfer control parameter U0, inte- ger valued, 1 to 4, initial profile parameter INPRO, and the number of radial stations, at least 2, but not more than 10, with arrays dimensioned as given. Control Parameters: U0: As in One-dimensional model, but controls only external heat transfer coefficient. (Note: due to error in initial evaluation of 144 jacket heat transfer coefficient, U0 should be 1.6 for most realistic simulation. Not a strong influence on results.) INPRO: Value Initial profile used 1 Uniform to wall at O, o 2 Center station at o, 0, others at X - 0., T - TcO' 3 Arbitrary Profile: Prompts for dissociation profile, from center to wall, then Temperature profile. Recomputes averages. 4 Cool sheath profile. As 1, but at wall station, X - 0., T - TcO' Averages recomputed. Other Error, run aborts. Outputs: At initial axial position and at print intervals there- after, prints 2, , , P, Tc’ Bulk Velocity, Flow area, Sonic Velocity, and Estimated Corrected Gas Dynamic Temperature, the last four based on the average composition and temperature, then the radial profile of dissociated fraction, from the center out to the wall, then the radial temperature profile in the same order. The local state variables should stack one above the other. Program 4: Dispersion Model.(Gear method, search for initial deriva- tives, on VAX) Two-Dimensional, Dispersion, Cooled, Reactive Surface. Inputs and Control Parameters: In file startg.dat: as Two-Dimen- sional Model, (one station valid) plus, IMSL DGEAR routine control parameters METH and MITER (see IMSL documentation), real-valued mass and 145 energy dispersion control parameters, and a real-valued debugging print control parameter. If dispersion control parameters positive, dispersion coefficient is diffusion coefficient or thermal conductivity time con- trol parameter; zero, plug flow in mass or energy; negative, dispersion coefficient is constant at absolute value in SI units. Debugging print will not print if parameter set to 0. In file search.dat: From central to outer, for dissociation at all stations then for temperature -- a low value, a high value, and the current value of a search parameter, and two integer valued search control parameters. Last, on a separate line, a print control parameter. The first search control parameter controls whether the search will be bounded if the associated variable is the next to invoke one of the runaway conditions. The second, if 1, will make the search bounded regardless of the history, and, if 0, the search will run freely. The first is set to the value of the second if the corresponding condition is not the first to reach runaway and to 1 if it is. The second search control value should be set to 1 only if you are sure, from watching the behavior of the search, that any invocation of the runaway conditions associated with that variable will set an absolute, global minimum or maximum of that parameter. In the sample case, quasi-one-dimensional, calling for a two-dimensional search, it was noticed that, as the search narrowed, the temperature runaway would only show up along a very nar- row, diagonal region, between regions where the dissociation runaway dominated. To keep the search from looping, the temperature search was strictly bounded. In addition, it proved faster to search more or less parallel to that diagonal, rather than strictly orthogonally. Hence, the transformed search coordinate matrix was devised (vide infra). If either 146 the high or low value of the search parameter is equal to the current value, the search will step the difference between the high and low if the the runaway indicates that the search should continue in the direc- tion corresponding to whichever is equal to the current value, regard- less of the value of the search control variables. If the first search control value is O, the search steps to the other value, if the control value is 1, it steps to half-way between, if the direction of the search is the other direction. If the print control parameter is not 0, the program will print to the screen or to the *.1og file after every print step. If it is l, the program will stop once it reaches the target reactor length. In file saves.dat: If defaulted, at least 4n+9 commas, will begin from startg.dat starting point. Else will start at last print position from last run. In file hmat.dat: Matrix for linear transformation of search para- meters to initial derivatives. Default, 4n2 commas, is identity matrix. Output: Before each attempt and restart, prints initial deriva- tives. When print control parameter non-zero, as Program 5, with surface conditions on all print intervals, and beneath each local condition, first and second derivatives. Program 5: Subroutine RK4: Final form of fourth-order Runge-Kutta in- tegration subroutine. Additional features: generalized through use of parameter to specify function to be integrated, cascading return capa- bility, integration to dependent variable value. 147 Area Function Program models as 0.025 m straight, circular tube, with flow area AK, first area input parameter, followed by a straight-cone converging section 0.025 m long to a throat with area A2, second area input para- meter, and another straight-cone diverging section 0.25 m long, back to flow area AK, and a straight tube with flow area AK thereafter. Solution Methods Differential Equations For Programs 1-3, the systems of coupled differential equations were solved using a standard fourth-order Runge-Kutta method. Program 5 is the latest version of the integrator. With bold face representing vector notation, for a series of coupled differential equations: dY(X)/dx -f(y.X) (A1) with initial conditions: Y(x0) - YO (A2) 148 y(xO + h) may be estimated by: k - f(yO + k1 / 2, x0 + h/2) (A4) k - f(y0 + k2 / 2, x0 + h/2) (A5) k - f(yO + k + h) (A6) 4 3’ x0 y(xO + h) - yo + (k1 + 2 k2 + 2 k3 + k4) h / 6 (A7) This approximation is improved by taking smaller values of h for the integration interval. In these programs, the integration is done over the print interval, then the step size is halved, and the integra- tion repeated. If the sum of the squares of the ratio of the difference between the values of each dependent variable after the two integrations to the change of the value over the final integration, summed over the dependent variables, is larger than the input convergence value, the integration interval is halved, and the integration is repeated until convergence. The stiffness of the equations is controlled by having the integrator check to see if the change of any of the temperatures over the integration interval will be greater than one-tenth of the value of the temperature. If it is, then the integration interval is halved until it is not, and the final value of the interval is used over the rest of the print interval. There is a floor, so that if the integration inter- val becomes less than one-millionth of the print interval, the computer prints out "INTEGRATION DIVERGED" and the integrator subroutine returns to the main program, with the final value of the overall integration interval made negative. The program may then prompt for new input (Pro- gram 1), or reduce the interval fed to the integrator, and repeat until 149 that becomes too small, at which time it will prompt for new input (Programs 2 and 3). Program 5 includes two new features that earlier versions of did not include. Instead of using the final length as the error flag, it calls for another integer input variable as an error flag, initially 0, returning it set to 1 if the integration diverged. The print statement was also removed. If that variable, through a common block, is set non- zero by the subroutine called to compute the derivatives, the integrator will return at once to the calling routine, with the error flag still set. Second, the subroutine requires input of an integer variable, zero or the index for one of the elements in the Y array and a target value for that element. If that variable is non-zero, the integration will run until the final value of the independent variable or until the desig- nated element passes the target value, whichever occurs first. If the index is 0, integration is strictly controlled by the independent vari- able. This was implemented because the backwards integrations, using an arc-length method, would sometimes overshoot z-O. The variable selected should be monotonic over the integration range, of course, though the integrator can handle either increasing or decreasing target variables. It was thought at first that an extension of the Runge-Kutta method, with dT/dz being used as a new dependent variable, might be applicable to solve the dispersion equations. However, this method provided unstable solutions, and a survey of the literature (Coste, gt a1. [1961]) showed that that was to be expected from forward-marching methods for solving the dispersion equations. Hence, Coste's method for solving such equations has been adopted. The non-dispersed equations are 150 solved first via a cell model, an implicit method, which, for appropri- ately chosen cell lengths, depending on the relationship between the dispersion coefficient and flow rate, approximates a dispersed flow reactor. Then, once the end of the reactor is reached, the values at the end of the reactor are used as input, and the exact equations integrated backwards along the reactor to give the exact profiles. To solve the cell model (Program 4), which entails the solution of a large number of non-linear equations, first, the cell length was chosen to fit Coste's parameter based on the bulk or average conditions at the input to the reactor. Then, the conditions at the end of the cell were computed based on the axial derivatives of the conditions computed using the initial conditions. Certain limits were applied, temperatures and dissociations were kept positive, for instance, to maintain realism, and the cell length and conditions at the end of the cell recalculated. Then, for each dependent variable, the difference between the guessed value and the recalculated value was computed for both cases, and, if the difference was too large, a new value of the dependent variable was extrapolated or interpolated, as appropriate, linearly. Then, the limits were applied, and the process repeated until the sum of squares converged. 151 Arc-Length Method For the reverse integration and for the later Gear model integra- tion, an arc-length method (after Bhatia and Hlavacec [1983]) was imple- mented to avoid stiffness problems. This method involves a new independ- ent variable 5 and redefines the equations as follows: 2)1/2 dz/ds - 1 / ( 1 + |df/dz| (A8) df/ds - df/dz / ( 1 + |df/dz|2)1/2 (A9) where |df/dz|2 is the sum of the squares of the elements of df/dz. This system of equations is integrated with respect to s, assuming as in- tegrating intervals over s the print interval, and printing after the value of' 2 has passed the print interval, regardless of whether the print interval has been hit exactly. Since dz/ds is always less than 1, it will print at least once in each print interval. To help limit the running time, the step size in s is adjusted after each call to the integrator to be the remaining length to the next print point divided by 16 dz/ds, as computed by finite difference over the last integration, or the print step, whichever is larger. This may negate the advantages of the arc-length method, and the adjustment, originally the full length divided by dz/ds, had to be reduced by a factor of 16 to avoid stepping too far past zero in some reverse integration cases. 152 Implicit Runge-Kutta Method When the dispersion model equations could not be integrated stably on Coste's backwards integration, using the standard, explicit Runge- Kutta integrator described above, an implicit Runge-Kutta method, as described by Cash [1975] was implemented. This method would change the equations (A3)-(A7) above to: kl - f(y0, x0) (A10) k2 - f(y(x0 + h) - k4 / 2, x0 + h/2) (A11) k3 - f(y(xO + h) - k2 / 2, x0 + h/2) (A12) k4 -f(y(x0 + h), x0 + h) (A13) y(x0 + h) - yo - (k1 + 2 k2 + 2 k3 + k4) h / 6 (A14) These equations were solved by assuming that, as an initial guess, y(xO + h) - yo + h k and using equations (All)-(A13) to compute the other 1 . k vectors, then substituting them into (A14) to obtain a new guess, and repeating until the new guess was close enough to the last one. When the integrators specially developed for this modeling proved intractably impractical, due to excessive computation time for appro- priate reactor lengths, an appropriate "canned" integrator was sought. The subroutine DGEAR from the IMSL library was selected. This subroutine uses, at operator option, either implicit Adams methods or Gear's stiff, backward differentiation methods, with a variety of user-selectable corrector iteration methods to solve the system of algebraic equations generated. 153 Algebraic Equations The surface conditions were estimated by solving the appropriate material and energy balances for the surface by a trial and error method. The first guess was that the surface conditions were those of the bulk gas or of the station closest to the wall, in a two-dimensional model. Then, the surface temperature needed to satisfy the energy balance, with generation at these conditions, was estimated, and the equilibrium conversion at this temperature was calculated, for the next guess. The signed error of closure in the material and energy balance equations was computed for each case, and the next guess taken from a bounded linear inter- or extrapolation, as appropriate, with bounds selected to maintain realistic values and to avoid runaway. Initial Conditions (Dispersion Model) After each print length integration is completed, the results are checked to see if each station's conditions, temperature and concentra- tion, were likely to be compatible with the infinity boundary condition. The following conditions were originally considered to be incompatible with the boundary conditions, suggesting the condition would run away, eventually going to plus or minus infinity: 1) The variable, if the first derivative did not change, would become negative in the next print length. 2) The first and second derivatives of the variable were both positive. 3) The first, second, and third derivative were of the same sign. The third derivatives were calculated by a finite-difference 154 method. 4) The square of the first derivative, divided by the product of the value and the second derivative, was positive and less'than unity. The first and third, if the derivatives were negative, were assumed to mean a negative runaway, so that the initial guess of the derivative involved had to be increased. The others, and the third with positive derivatives, were taken to indicate positive runaway, requiring a lower initial guess. Initial derivative values were assumed, and the equations integrated until one of the runaway conditions turned up for one of the values. Then the initial value of the derivative of that value was adjusted appropriately, and the integration repeated, until the integra- tion went the full length of the modelled reactor. The initial values were stepped along until both positive and negative runaway were indi- cated, then the point half-way between was selected. The variables would reach runaway in turns, though not strictly alternating, and they were clearly linked, so that finding a better value for one would shift the best value for another. This produced long running times, since the search is literally looking for an infinite integration length, even for a pseudo-one-dimensional, one station case, and the program was redesigned to allow for interruption and restarting. Inclusion of the last listed runaway condition produced dramatic changes in the performance of the search algorithm from that seen with- out that condition. That condition would occur much earlier than the others, and therefore a shorter integration would be needed to determine which section of the search area a point was in. The boundaries of the sections would change, as that condition might be satisfied earlier than one of the other conditions, and define a point as belonging in a dif- ferent section than the other condition would. However, that condition 155 produced false positive responses in certain regions. Inclusion of that condition would send the search one way when a longer integration would send the search in exactly the opposite direction. In a preliminary, one-dimensional, two-variable modeling, part of one section disappeared entirely, since the runaway conditions on one variable would never be satisfied before one or the other of the runaway conditions on the other were. The transition from one runaway to the other occurred at a finite reactor length. This forced the search technique to a halt, with a false convergence on the transition from positive to negative runaway. Hence, condition 4 was dropped from further searching, except as a means to establish lower bounds. Likewise, on experimentation with a two-dimen- sional model, condition 2 and the negative runaway from condition 3 were seen to create false results, so they were also dropped. This left only condition 1, with the coolant temperature as the floor for the tempera- tures, and zero, or the wall dissociation, if the dissociation in question was higher than that at the wall, for the dissociation values, to mark negative runaway and the positive side of condition 3 to mark positive runaway. u o d estart C a it Due to long normal running times, plus a number of bugs in the program that caused or required the program to stop, the program was equipped with a limited interrupt and restart capability. At each print point, the reactor state array was written into a special file called saves.dat. Each time the search algorithm was implemented, the values required for the algorithm were likewise written onto their own file, 156 search.dat. Each of these files was rewound before the new information was written. On starting the program anew, these files were read, and, if saves.dat was not defaulted, the integration resumed at the last print point before the previous run was halted. Search.dat had to be initialized with the intended start point and step sizes, and saves.dat could be set to default to begin the integration from z-0.0. Two problems with the program remain unresolved, though the program has been modified to recover from them. First, the search for the solu- tion to the surface boundary conditions may not converge. Second, the computer would report an "undefined exponentiation" error, and halt. It was discovered that a simple halt and restart, as described above, could usually integrate right over the problem. I presume that the DGEAR algorithm chooses the points at which it calls the differential equation computation subroutines in such a way that which specific points are evaluated changes with where the DGEAR initiation takes place. The program was restructured so that if the surface conditions search did not converge or the exponentiation probably responsible would be unde- fined the program would perform the equivalent of a restart from the last print location. 157 Transformed Search Coordinates The program calls a file hmat.dat, and reads off it a 2n by 2n matrix called HMAT, which holds the matrix for transforming the vector of search parameters, by matrix multiplication, to the initial deriva- tives used in the calculation. If the file is defaulted, with nothing but 4n2 commas, HMAT is assumed to be the identity matrix. A(3,3) Al A2 ALFA AMO AQ AR(Z) AS AVG B BQ BRAVO C(3.3) CHARLIE CP(I,T) CPM(X,T) CPV D2Y(47) 158 VARIABLE DICTIONARY (AVG) Dimensionless Outer Radius of Cell (DELTA, TVSS) Array Containing Constants for Rate Constant Calculation (AVG, DELTA2) Real Constant for Fewer Mode Problems Throat area for Area Function. (DELTA2) Real Value of Loop Counter For Fewer Mode Problems. Constant Area for Area Function Common Block Name Molecular Weight of undissociated gas (Prog. 4, DELTA) Parameter For Quadratic Fit In Temperature Search Area Computation Function Subroutine Surface Area In Square Meters Per Meter Length (Perimeter) Assumes Circular Cross Section Function Subroutine to Compute Mass-Averages (AVG) Dimensionless Inner Radius of Cell (Prog. 4, DELTA) Parameter for Quadratic Fit During Temperature Search Common Block Name Array Containing Constants for Heat Capacity Calculations Common Block Name Heat Capacity at Constant Pressure calculated as quadratic function of temperature T (I-: 1, Molecular Hydrogen; 2, Atomic Hydrogen; 3, ACP of reaction) Constants in Array C Heat Capacity of Hydrogen Mixture with Dissociation X at Temperature T (Prog. 1) Vibrational Heat Capacity of Molecular Hydrogen (Frog. 4) Array of derivatives from a DELTA2 call D3Y(47) DAR(Z) DELTA DELTA2 DGO DGEAR DH(T) DHO DH2HP(T) DHHZP DHS DHT DHV DISP(2) DISPM DR DS DTS DTSZ DTSO 159 (Prog. 4) Array of finite-difference computed second and third derivatives. Function Subroutine for finite-difference computation of derivative of area with respect to 2. Step size 1 pm. Subroutine that Computes Local Derivatives Based on First-Order Conservation Equations Subroutine to Compute Derivatives Including Second Derivatives using Arc-Length Method Standard Free Energy of Reaction at 298.15 K (Prog 4) ISML Gear Method Integration Subroutine Heat of Reaction Calculated From AC and DHO For Temperature T P Heat of Reaction Calculated at O K Diffusion Coefficient between Molecular and Atomic Hydrogen as Function of Temperature T Mass Dispersion Coefficient Energy Released By Surface Recombination (Prog. l, TVSS) Heat of Reaction Considering Vibrational Excitation (Prog. 1, DELTA) Heat of Reaction Considering Vibrational Excitation (Prog. 4, ROCKET) Two element array containing Mass and Energy Dispersion Control Variables. (Frog. 4, DELTA2) Mass Dispersion Coefficient Control Variable Radial Spacing Between Stations, Dimensionless Value of Coefficient for Conversion From Direct to Arc-Length Method Difference Between Actual Surface Temperature and That Needed to Balance Heat Flux Storage Variable for Surface Temperature Difference During Search Storage Variable for Surface Temperature Difference During Search DUMMY DX DXO DXL DXN(2) DXO DXR DXS DXSO DXW DY DY(N) DY4 DZ ECHO EDISP EGRESS EK(Y,Z) ERR ERRX EV(T) EVTV 160 (Frog. 4) Subroutine required for DGEAR integration (see IMSL documentation) (Prog. l, TVSS) Variable used to Compute Reaction Rate (RK4) Initial Overall Step Size (RK4) Dimensionless Step Size (Prog. 4) Array for Projection phase of Solution of Wall Boundary Conditions (RK4) Initial Interval for Control of Stiffness (DELTA) Change in Dissociation Based on Reaction and Axial Flow Alone Temporary Value, and Difference Between Real Surface Dissociation and That Necessary to Maintain Flux To Surface Storage Variable for DXS During Search Change in Dissociation From Surface Reaction (DAR) Step-size for finite-difference derivative calculation Array of Derivatives, parallel local Y(N) array (Frog. 1) Temporary Variable in computation of Vibrational Temperature Change (ROCKET) Interval over which Integration Occurs. Common Block Name (Frog. 4, DELTA2) Energy Dispersion Coefficient Control Variable (Prog. 4) Entry point for error recovery restart Function Subroutine Computes Kinetic Energy of Flow at length 2 and Conditions described by vector Y. Y(l)-X, Y(2)-T, Y(3)-P, assumes uniform (ROCKET, DELTA) Relative Error Tolerance. (RK4) Relative Difference Between Successive Integrations. (Prog. 4, DELTA) Error Tolerance on XS search (Prog. 1, DELTA, TVSS) Vibrational energy as function of Temperature T (Prog. l, TVSS) Temporary Variable in computation F0 FCN FD FFF(Y,Z) FI FOXTROT FR FXT(4) G2 GAMMA GOLF HI HMAT HO HOTEL 11 12 13 14 IB ICON ID 161 Molar Flow Rate of Gas into Plasma (Prog. 5) Dummy variable name for function to be integrated (Prog. 2) - FFF / Di (Prog. 2) Fanning Friction Factor Function Subroutine (Prog. 2, FFF) Temporary Variable (Prog. 2) Common Block Name (Prog. 2) Viscous Flow control variable Array for Values in Simplex Method Search Debugging Print Control Variable (Prog. 1) Fraction of Recombination Energy in Vibrational Excitation. (Prog. 4) Debugging Print Control Variable. (Prog. l, TVSS) Temporary Variable in Computation (Frog. 4) Common Block Name (Frog. 4) DGEAR Control Step Size (see IMSL documentation) Heat Transfer Coefficient From Outermost Station To Wall (Prog. 4, ROCKET) Matrix for Transformation of Search Parameters Coolant Heat Transfer Coefficient Common Block Name Loop Counter, Used Repeatedly Loop Counter. (DELTA2) Locates Dissociation Loop Counter. (DELTA2) Locates Temperature (DELTA2) Locates Temperature Derivative (DELTA2) Locates Dissociation Derivative (Frog. 4, ROCKET) Array of Second Search Control Parameters (Prog. 4) Print, stop control variable (Prog. 4) Flag for Dispersion check in search Algorithm IER IFl IF2 IJ IMAX IMIN INDEX INDIA INPR INPRO IS ISN IT IWK J1 J2 J5 J6 K(I,T) KEQ(T) KOUNT 162 (Prog. 4) DGEAR Error Flag Variable (see IMSL documentation) (Prog. 4) Flag Variable Used to Check for Non-Realistic Results (Prog. 4) Flag to tell whether to cascade return or call entry EGRESS on error result. (Frog. 5) Error Flag in Runge-Kutta integrator. Loop Counter Marker For Maximum Value During Simplex Search Marker For Minimum Value During Simplex Search (Prog. 4) DGEAR Control Parameter (see IMSL documentation) (Prog. 4) Common Block Name (Prog. 4) Default value for INPRO - 3 Initial Profile Control Variable (1-5) (Prog. 4) Array of First Search Control Parameters (Prog. 4) Temporary Variable in Search Algorithm (Frog. 5) Dummy Parameter for index of dependent variable that controls length of integration (- 0 for none) (Prog. 4) Array for DGEAR integration controls (see IMSL documentation) Loop Counter Index For Dissociation Value Index For Temperature Value Loop Control Variable in Wall Boundary Condition Computation Loop Control Variable in Wall Boundary Condition Computation (RK4) Loop Counter Reaction Rate Constants as Function of Temperature T. I- 1) Hydrogen Molecules as Third Body, 2) Hydrogen Atoms as Third Body, 3) Surface Reaction. Constants in Array A. Form k(T) = a Tb ec/T Equilibrium Constant at Temperature T (RK4) Loop Counter 163 KTH(T) Thermal Conductivity of Hydrogen as Function of Temperature T KTHY Heat Dispersion Coefficient KVT(Y) (Prog. l) Vibrational-Translational Relaxation Rate Constant KW Mass Transfer Coefficient From Gas To Wall L Loop Counter, Used Repeatedly Ll Leap Counter L2 Loop Counter L3 Loop Counter L5 Loop Control Variable in Wall Boundary Condition Computation L7 Loop Control Variable in Wall Boundary Condition Computation L8 Loop Control Variable in Wall Boundary Condition Computation LAST Marker For Last Point Changed in Simplex Method Search LCOUNT (Prog. 4, DELTA) Loop counter LT Loop Control Variable in Wall Boundary Condition Computation LX Loop Control Variable in Wall Boundary Condition Computation M2 (Prog. 2) Control Variable METH (Prog. 4) DGEAR Control Variable (see IMSL documentation) MIKE (Prog. 2) Common Block Name MITER (Prog. 4) DGEAR Control Variable (see IMSL documentation) MODE (Prog. 2) Control Variable MODEL (Frog. 4) Main Program Name N (DELTA, DELTA2, RK4, AVG, DUMMY) Length of Y Array Submitted N1 (ROCKET, RK4) Number of Radial Stations N2 (Prog. 4, ROCKET) - Nl*2+5, Location of Last Non-Derivative Element in Y Array. (DELTA, DELTA2) Number of Radial Stations. (RK4) Location of First Temperature Value N3 N4 N5 N6 N7 NCOUNT NF PD PI POPPA QI QO R R1 R2 RE RK4 RNHl RNH2 ROCKET RQ1 164 (Frog. 4, ROCKET) - 4*Nl+7, Total Number of Elements in Y Array. (DELTA) Location of First Temperature Value. (RK4) Location of Last Temperature Value. (Prog, 4, ROCKET) - N2+1, Number of First Derivative Element () in Y Array. (DELTA) - Number of Radial Stations Minus One, for Loop Counter that Computes all but Outermost Station (Prog. 4, ROCKET) - N4+1, Location of First dX/dZ Element in Y Array (Prog. 4, ROCKET) - N3-1, Location of Last dT/dZ Element in Y Array (Prog. 4, ROCKET) - N1*2, For Search Counters Loop Control Variable for Debugging Printing (Prog. 4, DELTA2) Length of Y Array Submitted to Subroutine Delta From Delta2 (Prog. 4, DUMMY) Array required for DGEAR integration (see IMSL documentation) Mathematical constant, a - 3.14159... Common Block Name Heat Flux From Gas To Wall Heat Flux From Wall Into Coolant Universal Gas Constant Dimensionless Inner Radius of Radial Cell Dimensionless Outer Radius of Radial Cell (Frog. 2, FFF) Reynolds Number Runge-Kutta-4 integration subroutine Radial Diffusion Mass Flux Into Cell From Next Inner Cell Radial Diffusion Mass Flux Outward From Cell (Prog. 1-3) Main Program Name (Prog. 4) Main Driver Subroutine Name Radial Thermal Flux by Conduction Into Cell From Next Inner Cell RQ2 SRC SUM(2) T0 TAIR TANGO TB TG THETA T00 TOL TS TSl T82 TSA TSI TSMIN TSN(2) TSO TSX(2) TVSS 165 Radial Thermal Flux by Conduction Outward From Cell (Prog. 4) Array of search parameters (Prog. 4, DELTA) Array Used For Computation of New Point in Simplex Method Search Standard Temperature for Thermodynamic Data, 298.15 K Molar Flow Rate of Cooling Air (Prog. 1) Common Block Name Average Temperature Between Radial Stations, Used to Compute Fluxes and Transport Properties Estimate of Gas-Dynamic Temperature based on Mass-Average conditions (Prog. l) Fundamental Vibrational Frequency of Hydrogen Molecule Expressed as a Temperature Far Upstream Temperature for Calculation of Danckwerts Boundary Conditions (RK4) Error Tolerance Surface Temperature Non-Reaction Surface Temperature Based on Flux Balances Storage Variable for Surface Temperature During Search (Prog. 4, DELTA) Surface Temperature From Flux Balances and Reaction at TS Storage Variable for Surface Temperature During Search Minimum Surface Temperature Value, Arbitrarily One-Tenth of Coolant Temperature Array for Projection phase of Solution of Wall Boundary Conditions Storage Variable for Surface Temperature During Search Array for Projection phase of Solution of Wall Boundary Conditions (Prog. 1) Function Subroutine to estimate initial value of TV based on dTV/dz - 0. at z - 0. U0 UB(Y,Z) UI VB VIS 166 (Prog. 2) Temperature From Which Wall Reactions Were Computed Coolant Heat-Transfer Coefficient Control Variable Function subroutine Computes Bulk Velocity at length 2 and Conditions described by vector Y. Y(l)-X, Y(2)-T, Y(3)-P, assumes uniform, Ideal Gas Law. Coolant Heat-Transfer Coefficient Control Variable Bulk Velocity (Used when use of DB seemed precluded) (Prog. 2, FFF) Viscosity of Reacting Mixture, computed from Prandtl number relations VMACH(Y,Z) Mach Number at length 2 and conditions described by vector VS(Y) VSS Y. Y(l)-X, Y(2)-T, Y(3)-P, assumes uniform, Ideal gas Function Subroutine (or statement function) Computes Velocity of Sound at Conditions described by vector Y. Y(l)-X, Y(2)-T, assumes uniform, Ideal gas Speed of Sound (Used when use of VS seemed precluded) W(47,15) Working Array X X0 X1 XB XC XEQ(T.P) XF XIT XL XS XSl (RK4) Initial Value of Independent Variable (Frog. 4, DUMMY) Dummy Parameter Far Upstream Dissociation for Calculation of Danckwerts Boundary Conditions (RK4) Intermediate Value of Independent Variable Average Dissociation Between Stations, Used For Calculating Fluxes (RK4) Counter for Debugging Output Equilibrium Dissociation at Temperature T and Pressure P (RK4) Final Value of Independent Variable (Progs. 1-3, RK4) Also Used as Flag for Divergent Integration (Prog. 5) Dimensionless Value of change of controlling dependent variable (RK4) Dimensionless Value of Independent Variable Surface Dissociation Surface Dissociation To Sustain Mass Flux To Surface XSN(2) XSX(2) XTN(2,4) Y0(47) Y1(47) YMIN YO 21 22 22D Z3 24 167 Array for Projection phase of Solution of Wall Boundary Conditions Array for Projection phase of Solution of Wall Boundary Conditions Array of X8 and TS Values used for Simplex Method Search Storage Variable for Wall Dissociation Value During Search Storage Variable for Wall Dissociation Value During Search Working Array: (Prog. 1) 1) X; 2) T; 3) P; 4) TV. (Prog. 2) l) X; 2) T; 3) P; 4) Tc (Prog. 3) 1) ; 2) ; 3) P; 4) Tc; 5 - 4+2n) odd elements X, even T. (Frog. 4, ROCKET, DELTA2) 1) Z; 2) ; 3) ; 4) P; 5) Tc; 6 - 5+N) X; 6+N - 5+2N) T; 6+2N) ; 7+2N - 6+3N) dX/dZ; 7+3N - 6+4N) dT/dZ; 7+4N)
. (Prog. 4, DELTA), elements 2 through 5+2N, renumbers. (RK4, AVG, VS, EK, UB, TVSS, DUMMY, Prog. 5): Dummy Name for Submitted Array. (Prog. 4) Y Array at initial location (Prog. 4) D2Y Array at previous print point. (Prog. 4) Minimum value for variables to invoke negative runaway. For dissociations - 0. Temperatures - coolant temperature Length at Call to Integrator routine (Progs. 1-3, DELTA, UB, EK, AR, DAR) Length Down Reactor. (Prog.4, ROCKET, DELTA2) Independent Variable for Arc-Length Method Integration. (AR) Length at start of Area Defining Function (ROCKET) Final Value of Independent Variable for Integration. (AR) Second defining length for Area Function. Length at Start of Integration Step (AR) Third Defining Length for Area Function (AR) Fourth Defining Length for Area Function (AR) Dummy Variable for Area Computation ZB ZL ZP ZPR ZQ 168 (AR) Dummy Variable for Area Computation Total Reactor Length Print Interval Print Length (Prog. 5) Target value for dependent variable Y(IT) (Prog. 4) - -l., for early restarting procedure APPENDIX B Program Program Program Program Program APPENDIX B COMPUTER PROGRAM LISTINGS Contents Vibrational Modelling Reactive Surface Two-Dimensional Model Axial Dispersion Model Runge-Kutta-4 Integrator 169 170 175 182 189 206 170 Program 1 Vibrational Modelling One—Dimensional, Plug Flow, Adiabatic, Gas-Phase Reaction, Inviscid Flow 00000 000 0000 0' 103 20 C 171 PROGRAM ROCKET(INPUT.OUTPUT) PROGRAM TO MODEL HYDROGEN ATOM RECOMBINATION IN A PLUG-FLOW REACTOR ONE°DIMENSIONAL ADIABATIC MODEL HITH INTERNAL VIBRATIONAL ENERGY CONSIDERATION. COMMON/ALFA/A2.AK COMMON/BRAVO/AMO,FO.R COMMON/CHARLIE/DH0.0GO.TO.C(3.3).GAMMA DIMENSION v(4).w(4.9) DATA R.TO/8314.34.298.15/ DATA C/26880..4.347.-2.645E-04.2O7ee..o..o./ DATA DHO.DCO/4.3141an.4.0643E08/ DATA AMO/2./ OH(T)-DHO+T#CP(3,T) CP(I.T)=C(1.I)+T‘(C(2.I)+TPC(3.I)) VELOCITY OF SOUND COMPUTED FROM IDEAL GAS RELATIONSHIPS VS(Y)'SORT((1.0+Y(1))*R‘Y(2)/AMO/(1.-R*(1.+Y(1))/ +(2.‘Y(1)‘CP(2.Y(2))+(1.-Y(1))’CP(1.Y(2))))) VMACH(Y,Z)=UB(Y.2)/VS(Y) INITIALIZATION INCLUDING EVALUATION OF DELTA CP OF REACTION AND STANDARD SET OF INITIAL CONDITIONS AND SYSTEM PARAMETERS. DD 5 1:1.3 C(I.3)‘(2.*C(I.2)-C(I.1))/I Zao. Y(1)=o.1 Y(2)'1200. v(3)=3530. Y(4)=1200. FO=3.82E-7 ERR-1.E-6 AK=A2=3.aE-4 GAMMA=O. CONTINUE READ t.z.v.FO.CAMMA.ERR.AK.A2.ZL.2P IF(Z.LT.O.)GO TO 102 PRINT t,Z.Y,FO.GAMMA,ERR,AK,A2.ZL.ZP IF(V(4).LT.O.)v(4)-Tvs5(v) PRINT 20.Z.Y.UB(Y.Z).AR(Z).VS(Y).VMACH(Y.Z).TVSS(Y) FORMAT(3X.1OG12.5) C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C 101 102 0000000 ZF-AMIN1(Z+ZP.ZL) CALL RK4(z.zE.v.w.4.ERR) IF(ZF.LT.O.)GO TO 103 PRINT 20.Z.Y.UB(Y.Z).AR(Z).VS(Y).VMACH(Y,Z).TVSS(Y) IF(Z.LT.ZL)GO TO 101 GO TO 103 CONTINUE END FUNCTION AR(Z) REACTOR CROSS-SECTIONAL AREA AS FUNCTION OF LENGTH MODEL USED HAS 2.5 CM STRAIGHT LENGTH OF AREA AK (USER SELECTED) STRAIGHT CONICAL CONVERGING SECTION 2.5 CM LONG TO AREA A2 (ALSO USER SELECTED) AND STRAIGHT CONICAL DIVERGING SECTION 25 CM LONG BACK TO AREA A2. AK. A2 DEFAULTS ARE 3.8 E-4 SO M. 172 COMMON/ALFA/A2.AK DATA 21.22.23.24/.O25..05..05..30/ AR-AK IF(AK.EO.A2)RETURN IF((Z.LE.Z1).DR.(Z.GE.24))RETURN 100 IF(Z.LT.Z2)GO TO 101 IF(Z.GT.23)GO TO 102 AR-A2 RETURN 101 ZA'Z1 ZB=22 GO TO 103 102 2A'24 28.23 103 AR=(SORT(AK)+(SORT(AK)-SQRT(A2))‘(Z-ZA)/(ZA-ZB))'*2 000 0000 0000 0000 000 RETURN END FUNCTION DAR(Z) DERIVATIVE OF AREA VERSUS LENGTH BY FINITE DIFFERENCE METHOD STEP SIZE 1.E-6 M Dv=1.E-6 DAR-(AR(z+Ov)-AR(Z))/Dv RETURN END FUNCTION UB(Y.Z) BULK VELOCITY BY IDEAL GAS LAW COMMON/BRAVD/AMO.FO.R DIMENSION Y(A) UB=FO*(1.+Y(1))*R*Y(2)/Y(3)/AR(Z) RETURN END SUBRDUTINE DELTA(Z.N.Y.DY) COMPUTATION OF FIRST-ORDER DERIVATIVES OF CONCENTRATION TEMPERATURES AND PRESSURES FROM CONSERVATION RELATIONSHIPS IMPLICIT REAL (K) COMMON/aRAvo/AM0.FO.R COMMON/CHARLIE/DHO.DCO.TO.c(3.3).GAMMA COMMON/TANGO/THETA.A(3.2) DIMENSION Y(N).DY(N) RATE CONSTANTS DERIVED FROM SHUI AND APPLETON AND SHUI DATA LEAST-SQUARES FIT TO FORM OF A(I.1)*T*tA(I.2)*EXP(A(I.3)/T) DATA A/2.014E11.-O.7497.-20.43.8.816E11.-O.7219.'92.09/ DATA THETA/5990./ RATE AND EQUILIBRIUM CONSTANT EXPRESSIONS K(I,T)-A(1.I)tTttA(2.I)fiEXP(A(3.I)/T) KEO(T)-1.0132SEOSPEXP((~DGO/TO+C(1.3)‘ALOG(T/TO)+(DHO/T/TO+ +C(2.3)+C(3.3)*(T+TO)/2.)t(T-T0))/R) DH(T)-DHO+TtCP(3.T) CP(I.T)cc(1.I)+T*(C(2.I)+T*C(3.I)) KINETIC ENERGY OF FLOH AS FUNCTION OF X. T. P. AND 2 173 c EK(Y.2)=1./(1./(AMO‘UB(Y.Z)**2)-1./R/Y(2)/(1.+Y(1))) C c EXPRESSIONS FOR VIBRATIONAL ENERGY. HEAT CAPACTIY. RELAxATION c RATE AND HEAT OF REACTION c EV(T)-RtTHETA/(EXP(THETA/T)-1.) CPV(T)'(EV(T)/T)t'2tEXP(THETA/T)/R KVT(Y)=((1.-Y(1))*25300.*EXP(-100./(Y(2)"'(1./3. )))+2*Y(1)* +6.60E09tEXP(-3800./Y(2))/Y(2))*Y(3)/(1.+Y(1)) DHV=4.318852£OB+1.5thY(2)-EV(Y(4)) c C MATERIAL BALANCE EOUATION INCLUDING RATE EXPRESSION C DY(1)-AR(Z)/FO*(K(2.Y(2))E2.*Y(1)+K(1.Y(2))*(1.-Y(1)))#(KEO(Y(2))* +(1.-Y(1)*t2)/Y(3)-4.*Y(1)‘Y(1))-(Y(3)/R/Y(2)/(1.+Y(1)))t'3 c C ENERGY BALANCE EQUATIONS INCLUDING HEAT OF REACTION. KINETIC ENERGY C AND VIBRATIONAL ENERGY TERMS c DY4t-((1.-Y(1))tKVT(Y)t(EV(Y(4))-EV(Y(2)))/UB(Y.Z)+ +OHVtGAMMAtDY(1)) DY(2)=((EK(Y.Z)/(1.+Y(1))+DHV)*DY(1)-EK(Y.Z)’DAR(Z)/AR(Z)+ +DY4)/(EK(Y.Z)/Y(2)+R¢(7.+3.‘Y(1))/2.)*(-1.) DY(3)'-EK(Y.Z)‘(DY(2)/Y(2)+DY(1)/(1.+Y(1))-DAR(Z)/AR(Z))* +Y(3)/R/Y(2)/(1.+Y(1)) DY(4)-(DY4)/(1.-Y(1))/CPV(Y(4)) RETURN END SUBROUTINE RK4(x.xr.Y.V.N.TOL) DIMENSION Y(N).V(N.9) c C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C on-xE-x DXL-i. OO 6 I-1.N 6 w(I.9)=Y(I) 100 xL-O. DO 5 1:1.N S V(I.1)=Y(I) 101 CALL OELTA(x+chon.N.V(T.1).V(1.2)) C C CHECKS FOR STEEP GRADIENTS. AND IF STEP SIZE BECOMES TOD SMALL. C RETURNS TO MAIN PROGRAM MITH NEGATIVE FINAL VALUE FOR INDEPENDENT c VARIABLE C 105 IF(ABS(U(2.2)*DXLRDXD).GT.O.1*U(2.1))GO TO 106 IF(ABS(U(4.2)*DXL‘DXO).LE.0.IRV(4.1))GO T0 T04 106 DXL-DXL/2. GO TO 105 104 IF(DXL.GT.2.**(-20))GO TO 114 PRINT -.- INTEGRATION DIVERGED' xr--xr RETURN 114 XL-XL+DXL/2. XTIX+XL*DXO Dx-DXLtho DO 10 I-1.N 1O M(I.6)-H(I.1)+Dx-V(I.2)/2. CALL DELTA(x1.N.M(1.6).w(1.3)) 2O 3O 4O 60 50 103 70 174 DO 20 I=1.N H(I.7)=H(I.1)+DX*H(I.3)/2. CALL OELTA(x1.N.V(1.7).M(1.4)) DO 30 I=1.N U(I.8)=W(I.1)+DX*U(I.4) XL=XL+OXL XO=X+XL*DXD CALL DELTA(XO.N.M(1.B).V(1.5)) DO 40 I=1.N V(I.1)'V(I.1)+(U(I.2)+2.‘(V(I.3)+U(I.4))+H(I.5))*Dx/6. IF(XL.LT.1.)GO TO 101 DXL-DXL/2. ERR=O. DO so 1-1.N IE(ABs((w(I.1)-Y(I))/Y(I)).LT.TOL)GD TO 60 ERR=ERR+((U(I.9)-H(I.1))/(U(I.1)-Y(I)))*‘2 U(I.9)=V(I.1) CONTINUE IF(ERR.GT.TDL)GO T0 100 x-xF DO 70 1:1.N Y(I)-M(I.9) RETURN END FUNCTION Tvss(Y) C C COMPUTATION OF STEADY-STATE VALUE OF VIBRATIONAL TEMPERATURE C BALANCES EXCITATION AND RELAXATION RATES. C C IMPLICIT REAL (K) COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.DGO.TO.C(3.3).GAMMA COMMON/TANGO/THETA.A(3.2) DIMENSION Y(4) C RATE AND EQUILIBRIUM CONSTANT EXPRESSIONS C K(I.T)-A(1.I)tTttA(2.I)tEXP(A(3.I)/T) KEO(T)=1.01325E05tEXP(('DGO/TO+C(1.3)*ALOG(T/TO)+(DHO/T/TO+ +C(2.3)+C(3.3)#(T+TO)/2.)t(T-TO))/R) CP(I,T)=C(1.I)+Tt(C(2.I)+T*C(3.I)) EV(T)-RtTHETA/(EXP(THETA/T)-1.) KVT(Y)=((1.-Y(1))R25300.*EXP(-1OO./(Y(2)¢*(1./3.)))+2*Y(1)t +6.60E09-EXP(-3800./Y(2))/Y(2))*Y(3)/(1.+Y(T)) DHT=4.318852E08+1.5thY(2) DX-(K(2.Y(2))*2.*Y(1)+K(1.Y(2))‘(1.-Y(i)))*(KEO(Y(2))* +(1.-Y(1)"2)/Y(3)-4.*Y(1)*V(1))*(V(3)/R/Y(2)/(1.+Y(1)))#t2 GOLFIGAMMAth/KVT(Y)/(1.-Y(1)) EVTVI(EV(Y(2))-GDLF*DHT)/(1.-GOLF) TVSSITHETA/ALOG(1.+R¢THETA/EVTV) RETURN END 175 Program 2 Reactive Surface One-Dimensional, Plug Flow, Cooled, Surface Reaction, Viscous Flow 00000 000 (II 0000 103 00000 20 C 176 PROGRAM ROCKET(INPUT.OUTPUT) PROGRAM TO MODEL HYDROGEN ATOM RECOMBINATION IN A PLUG-FLOW REACTOR ONE-DIMENSIONAL. COOLED REACTOR WITH VISCOUS DRAG OPTION INCLUDES OPTIONS FOR SURFACE REACTION CDMMDN/ALFA/A2.AK COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.TD.C(3.3).GAMMA COMMON/ECHO/TAIR.UO COMMON/POPPA/PI CDMMDN/FOXTROT/FR COMMON/MIKE/MODE.M2.TS DIMENSION Y(4).w(4.9) DATA R.TO/8314.34.298.15/ DATA C/2GBBO..4.347.-2.645E-04.20786..O..o./ DATA DHO.DGo/4.3141EOB.4.0648EOB/ DATA AMO/2./ DH(T)=DHO+T#CP(3.T) CP(I.T)=C(1.I)+T*(C(2.I)+T‘C(3.I)) VELOCITY OF SOUND BY IDEAL GAS RELATIONS VS(Y)-SORT((1.O+Y(1))*R*Y(2)/AMO/(1.-R*(1.+Y(1))/(2.PY(1)# +CP(2.Y(2))+(1.-Y(1))*CP(1.Y(2))))) VMACH(Y.Z)'UB(Y.Z)/VS(Y) INITIALIZATION INCLUDING EVALUATION OF DELTA CP OF REACTION AND STANDARD INITIAL CONDITIONS AND SYSTEM PARAMETERS. DO 5 121.3 C(I.3)=(2.*C(I.2)-C(I.1))/I UO=1.6 FR=O. PI=4.-ATAN(1.O) z=o. Y(1)=o.1 Y(2)-12oo. Y(3)-3530. Y(4)-4oo. FO=3.82E-7 ERR-1.E-6 AK=A2-3.BE-4 TAIR=1.7E-4 CONTINUE READ t.Z.Y.FO.ERR.AK.A2.ZL.ZP.TAIR.UO.MODE.M2.FR IF(Z.LT.O.)GD TO 102 PRINT *.Z.Y.FO.ERR.AK.A2.ZL.ZP.TAIR.UO.MODE.M2.FR ESTIMATION OF THE GAS-DYNAMIC TEMPERATURE THAT THE EXPERIMENTAL ASSUMPTION OF ZERO DISSOCIATION WOULD REPORT UNDER THE CONDITIONS OF THE MODEL. TG=Y(2)‘5.*(7.+3.*Y(1))*(1.+Y(1))/7-/(5.+Y(1)) Ts-Y(2) PRINT 20.2.Y.UB(Y.2).AR(z).YS(V).TG FORMAT(3X.10612.5) C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C 101 104 ZE-AMIN1(z+zP.2L) Dz-zP 22-AMIN1(Z+DZ.ZF) 177 CALL RK4(Z.Z2.Y.W.4.ERR) IF (22.LT.O) DZ=DZ/16. IF (DZ.LT.1.E-7*ZP) GO TO 103 IF (Z.LT.ZF) GO TO 104 C C ESTIMATION OF THE GAS-DYNAMIC TEMPERATURE THAT THE EXPERIMENTAL C ASSUMPTION OF ZERO DISSOCIATION WOULD REPORT UNDER THE CONDITIONS C OF THE MODEL. C TG=Y(2)t5.*(7.+3.tY(1))t(1.+Y(1))/7./(5.+Y(1)) PRINT 20.Z.Y.UB(Y.Z).AR(Z).VS(Y).TG.TS IF(Z.LT.ZL)GO TO 101 GO TO 103 102 CONTINUE END FUNCTION AR(Z) C C CROSS-SECTIONAL AREA OF REACTOR AS FUNCTION OF LENGTH C MODEL USED WAS 2.5 CM LONG STRAIGHT TUBE OF AREA AK (USER SELECTED) C THEN STRAIGHT CONICAL CONVERGING SECTION 2.5 CM LONG TO AREA A2 C (ALSO USER SELECTED) FOLLOWED BY A 25 CM LONG STRAIGHT CONE DIVERGING C SECTION BACK TO AREA AK. DEFAULTS AK. A2 3.8E-4 SQ M. c COMMDN/ALFA/A2.AK DATA 21.22.23.24/.025..Os..05..ao/ AR=AK IF(AK.EQ.A2)RETURN IF((2.LE.21).OR.(2.GE.24))RETURN 100 IF(z LT.22)GO T0 101 IF(Z.GT.23)GO TO 102 AR-A2 RETURN 101 zA-z1 ZB=22 GO TO 103 102 ZA-24 23:23 103 AR-(SORT(AK)+(SORT(AK)-SORT(A2))‘(Z-ZA)/(ZA-ZB))"2 RETURN END FUNCTION DAR(Z) DERIVATIVE OF AREA WITH RESPECT TO LENGTH FINITE DIFFERENCE METHOD. STEP 1.E-6 M 0000 DY-1.E-06 DAR=(AR(2+DY)-AR(2))/DY RETURN END FUNCTION UB(Y.z) BULK VELOCITY BY IDEAL GAS LAW 000 COMMON/BRAVO/AMO.FO.R DIMENSION Y(4) UB-F0t(1.+v(1))PROY(2)/Y(3)/AR(Z) RETURN END SUBROUTINE DELTA(Z.N.Y.DY) 00 COMPUTATION OF FIRST-ORDER DERIVATIVES BY CONSERVATION RELATIONS 178 IMPLICIT REAL (K) COMMON/BRAVO/AMO.FD.R COMMON/CHARLIE/DHO.TO.C(3.3).GAMMA COMMON/ECHO/TAIR.UI COMMON/POPPA/PI COMMON/MIKE/MODE.M2.TS DIMENSION A(3,3).Y(N).DY(N) DATA UGO/4.0648E08/ RATE CONSTANTS DERIVED FROM SHUI AND APPLETON AND SHUI DATA LEAST-SQUARES FIT TO FORM OF A(I.1)‘T‘*A(I.2)*EXP(A(I.3)/T) 0000 DATA A/2.014E11.-O.7497.-20.43.8.816E11.-O.7219.-92.09.3.855E18. +-3..~2900./ c C RATE AND EQUILIBRIUM CONSTANT EXPRESSIONS C K(I.T)-A(1.I)FT"A(2.I)*EXP(A(3.I)/T) KEQ(T)=1.O1325605FEXP((-DGO/TO+C(1.3)FALOG(T/TO)+(DHO/ +T/TO+C(2.3)+C(3.3)‘(T+T0)/2.)‘(T-TO))/R) DH(T)=DHO+T*CP(3.T) CP(I.T)=C(1.I)+TF(C(2.I)+T*C(3.I)) C C KINETIC ENERGY OF FLOW AS FUNCTION OF X. T. P. AND 2 c EK(V.Z)'1./(1./(AMOFUB(Y.Z)‘*2)-1./R/Y(2)/(1.+Y(1))) C C MATERIAL BALANCE EQUATION INCLUDING RATE EXPRESSION C AND COMPUTATION OF SURFACE REACTION OPTIONS CONTROLLED BY C PARAMETERS MODE AND M2. LINEAR INTER- AND EXTRAPOLATION TO 0 FIND SURFACE CONDITIONS. C A(1.3)-3.855E1B LOOP-1 C C ASSUME CIRCULAR REACTOR CROSS SECTION TO RELATE SURFACE AREA AND C CROSS-SECTIONAL AREA c AS=2.*SORT(PI#AR(Z)) IF(MODE.EO.1)A(1.3)-O. FD=FFF(Y.Z)*PI/AS C C DETERMINATION OF HEAT-TRANSFER COEFFICIENT c IF UI LT 0.. USE ABSOLUTE VALUE IN SI UNITS C IF UI - o.. AOIABATIc c IF UI GT 0.. MULTIPLY BY NUSSELT NUMBER COMPUTATION C IF(UI)111.112.113 112 Ho-O. HIsT. GO TO 114 111 HI=HO=ABS(UI) GO TO 110 113 HOsUI:47.Oe HI-UI 114 HIFHIF3.66F1.8436-3*Y(2)*‘.8#PI/AS 110 TST-Ts-(HO#Y(4)+HI*Y(2))/(HO+HI) KW83.66tPIt1.264E-3'Y(2)“O.72/(1.+Y(1))/R c C MATERIAL BALANCE BASED ON FLOW ANO GAS-PHASE REACTION 150 120 0000000 122 121 130 151 C 179 DXVIAR(Z)/FO*(K(2.Y(2))*2.‘Y(1)+K(1.Y(2))'(1.'Y(1)))* +(KEO(Y(2))*(1.-Y(1))/Y(3)*(1.+Y(‘))’4-*Y(1)*Y(1))‘ +(Y(3)/R/Y(2)/(1.+Y(1)))*’3 DXSO=Y(1) XWO=O. XW=Y(1) XWI-XW TS=AMIN1(IS,Y(2)+10000.) TS=AMAX1(TS.Y(4)/2.) TSI=TW=TS IF(MODE.EO.4)TW-Y(2) DXW=AS/FO*K(3.TW)‘(KEO(TW)/Y(3)*(1.-XW*XW)- +4.*XW"2)#(Y(3)/R/TW/(1.+XW))*‘2 DY(1)=DXV+DXW ENERGY BALANCE EQUATION INCLUDING HEAT OF REACTION. KINETIC ENERGY AND THERMAL ENERGY TERMS ENERGY TRANSFER FROM SURFACE REACTION CONTROLLED BY PARAMETERS SURFACE CONDITIONS DETERMINED BY SOLUTION OF CONSERVATION LAWS BY LINEAR INTER- AND EXTRAPOLATIONS. DHS=DXW*DH(TS)*FO IF(MODE.GE.3)TS=TS1-DHS/AS/(H0+HI) DTS=TS-TSI IF(ABS(DTS).LE.1.)GO TO 130 GO TO (121.122).LO0P IF(TSI EQ.TSO)GO TO 130 IF(DTS.EQ.DTSO)GO TO 130 TS=TSI+DTS*(TSO-TSI)/(DTS-DTSO) LOOP=2 TSD=TSI DTSOsDTS GO TO 120 DXSI(Y(1)/(1.+Y(1))-XW/(1.+XW))+DXWFFO/KWF(1.+XW/(1.+XW)) IF(M2.EQ.0)GO TO 151 IF(DXS.EQ.DXSO)GO TO 151 XWFXW+DXS*(XWO-XW)/(DXS-DXSD) XWIAMAX1(XW.O.) IF(ABS(XWO-XW).LE.XW*O.1)GO To 151 XVOsXVI DXSOsst GO TO 150 OI-(Y(2)-TS)*HI*AS IF(MOOE.EQ.2)QI-QI+OHS QO=QI-OHS DY(2)=((EK(Y.Z)/(1.+Y(1)))‘DY(1)+DH(Y(2))‘DXV-EK(Y.Z)*DAR(Z)/ +AR(2)+OI/FO+FD‘UB(Y.Z)*‘4/R/Y(2)/(Y(1)+1.))/(EK(Y.Z)/Y(2)+ +(1.-Y(1))*CP(1.Y(2))+2.*Y(1)*CP(2.Y(2)))*(-1.) DY(3)--EK(Y.Z)#(DY(2)/Y(2)+DY(1)/(1.+Y(1))-DAR(Z)/ +AR(Z)+UB(Y.Z)#‘2RFD)‘Y(3)/R/Y(2)/(1.+Y(1)) C ENERGY BALANCE EQUATION ON COOLING AIR C INCLUDES TERM FOR RADIATIVE-CONVECTIVE LOSS TO ENVIRONMENT C USES DIATOMIC KINETIC THEORY VALUE FOR AIR HEAT CAPACITY C DY(4)-(OO+1.64*(300.-Y(4)))/TAIR/3.5/R RETURN END SUBROUTINE RKA(X.XF.Y.V.N.TOL) 180 C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C 100 5 101 C DIMENSION Y(N).W(N.9) DXO=XF-x DXL-T. CALL DELTA(X.N.Y.W(1.8)) KOUNT=-1 KOUNT-KOUNT+1 XC=.1 XL=O. DO 5 I=1.N V(I.2)-V(I.B) H(1.1)=Y(I) GO TO 105 CONTINUE CALL DELTA(X+XL#DXO.N.W(1.1).W(1.2)) C CHECKS FOR STEEP GRADIENTS. REDUCES STEP SIZE IF FOUND C IF STEP SIZE REDUCED TOO FAR. RETURNS. WITH NEGATIVE FINAL C X-VALUE AS FLAG C 105 106 104 110 10 2O 30 40 1000 1001 51 50 102 60 103 IF(ABS(W(2.2)*DXL*DXO).GT.O.1RW(2.1))60 To 106 IF(ABS(W(N.2)*DXLFDXO).LE.O.1FW(N.1))GO TO 104 DXLxDXL/z. GO TO 105 IF(DXL.GT.2.**(-20))GO TO 110 PRINT ~.- INTEGRATION DIVERGED" XFx-XF RETURN XLFXL+DXL/2. X1=X+XLFDXO DO 10 I-1,N W(l.6)=W(I.1)+DXLtDXOFW(I.2)/2. CALL DELTA(X1.N.W(1.6).W(1.3)) DO 20 I=1,N W(I.7)-W(I.1)+DXLFDXO*W(I.3)/2. CALL DELTA(X1.N.W(1.7).W(1.4)) 00 so I-1.N W(I.7)8W(I.1)+DXLFDXOFW(I.4) XL-XL+OXL/2. CALL DELTA(X+XL*DXO.N.W(1.7).W(1.5)) DO 40 I-1.N W(I.1)=W(I.1)+(W(I.2)+2.*(W(I.3)+W(I.4))+W(I.5))‘DXL¢DXO/6. IF (XL-xc) 1001.1000.1000 xc-XC+1./10. CONTINUE IF(XL.LT.1.)GO T0 101 DXL-DXL/z. IF(KOUNT.EQ.Q)GO TO 102 ERRsO. DO so I=1.N IF(ABS((W(I.9)-Y(I))/W(I.9)).LE.1.E-8)GO TO 51 ERR-ERR+((H(I.9)-H(I.1))/(w(1.9)-Y(I)))**2 CONTINUE CONTINUE IF(ERR.LE.TOL)GO TO 103 DO 60 I-1.N v(I.9)-M(I.1) GO TO 100 x-XF DO 70 I-1.N 70 000000 101 181 Y(I)=W(I.9) RETURN END FUNCTION FFF(Y.Z) COMPUTATION OF FANNING FRICTION FACTOR ITERATIVE SOLUTION OF RECURSIVE EQUATION FROM PERRY FDR SMOOTH PIPE TURBULENT FLOW LAMINAR FLOW USES F'64/RE COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.TO.C(3.3).GAMMA COMMON/POPPA/PI COMMON/FOXTROT/FR DIMENSION Y(4) CP(I.T)=C(1.I)+Tt(c(2.I)+TcC(a.I)) FFFcFR IF (FR.EQ.O.)RETURN HIF1.843E-3*Y(2)**.8 VIS=O.7*HI/(CP(1.Y(2))F(1.-Y(1))+2.tY(1)*CP(2.Y(2)))’AMO REFFOFAMO/VIS/SORT(PIFAR(Z)) FFF=FR¢64./RE IF(RE.LT.2100.)RETURN FI=FFF/FR FFF=(2.*ALOGTO(REFSQRT(FI))-O.8)**(-2) IF(ABS((FFF-FI)/FFF).GT.1.E-6)GO TO 101 RETURN END 182 Program 3 Two-Dimensional Model Two-Dimensional, Plug Flow, Cooled, Surface Reaction, Inviscid Flow 183 PROGRAM ROCKET(INPUT.OUTPUT) c C PROGRAM TO MODEL HYDROGEN ATOM RECDMBINATION IN A PLUG-FLOW REACTOR c TWO-DIMENSIONAL PLUG-FLOW MODEL C COMMON/ALFA/A2.AK COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.To.c(3.3).GAMMA COMMON/ECHD/TAIR.UQ CDMMON/POPPA/PI DIMENSION Y(24),W(24.9) EXTERNAL DELTA DATA R.To/8314.34.298.15/ DATA C/2BBBO..4.347.-2.645E-04.20786..o..O./ DATA DHO.DGO/4.3141E08.4.0648E08/ DATA AMo/2./ DH(T)-DHO+T#CP(3.T) CP(I.T)=C(1.I)+T‘(C(2.I)+T*C(3.I)) C c VELOCITY OF SOUND FROM IDEAL GAS RELATIONS C VS(Y)=SORT((1.0+Y(1))*R*Y(2)/AMO/(1.-R*(1.+Y(1))/(2.*Y(1)* +CP(2.Y(2))+(1.-Y(1))‘CP(1.Y(2))))) VMACH(Y.2)-UB(Y.z)/VS(Y) INITIALIZATION INCLUDING EVALUATION OF DELTA CP OF REACTION AND STANDARD SET OF INITIAL CONDITIONS AND SYSTEM PARAMETERS. DO 5 1:1.3 C(I.3)=(2.¢C(I.2)-C(I.1))/I 2:0. Y(1)=o.1 Y(2)=12oo. Y(3)=3530. Y(4)=4oo. Fo-3.B2E-7 ERR=1.E-6 AK=A2-3.BE-4 TAIR-1.7E-4 P184.*ATAN(1.0) 01 0000 UO=1.6 103 CONTINUE C C INPUT OF INITIAL CONDITIONS AND SYSTEM PARAMETERS C READ *.Z.(Y(L).L'1.4).FO.ERR.AK.A2.ZL.ZP.TAIR.UO.INPRO.N1 IF(Z.LT.O.)GO TO 102 PRINT *.Z.(Y(L).L'1.4).FO.ERR.AK.A2.ZL.ZP.TAIR.UO.INPRO.N1 N2=N1t2+4 GO TO (111.112.113.111).INPRO C c IF INPRO = 1. INITIAL PROFILES ARE UNIFORM ACROSS REACTOR c 111 OD 110 I-5.N2.2 Y(I)=Y(1) 110 Y(I+1)-Y(2) GO TO (150.150.150.114).INPRO c c IF INPRO - 4. INITIAL PROFILES ARE UNIFORM EXCEPT THAT OUTER-MOST C STATION IS AT COOLANT TEMPERATURE AND ZERO DISSOCIATION C 114 Y(N2-1)=o. Y(N2)-Y(4) ‘00000 12 120 d000000 13 150 00000 30 20 C 184 GO TO 150 IF INPRO = 2. INITIAL PROFILE IS ZERO-DISSOCIATION AND COOLANT TEMPERATURE ACROSS REACTOR EXCEPT FOR CENTRAL STATION. WHICH IS AT INPUT TEMPERATURE AND DISSOCIATION. Y(s)=v(1) Y(6)=Y(2) DO 120 I=7.N2.2 Y(I)=o. Y(I+1):Y(4) GO TO 150 IF INPRO'3. PROGRAM PROMPTS FOR ARBITRARY INITIAL PROFILES. IF CONTINUATION OF A PREVIOUS RUN IS OESIRES. DEFAULTS TO LAST PROFILE. SO STRING OF '."S WILL CONTINUE RUN. NO DEFAULT ON INITIALIZATION FOR FIRST RUN. PRINT *.' X PROFILE“ READ ‘.(Y(L).L85.N2.2) PRINT ¢.(Y(L).L=5.N2.2) PRINT i.” T PROFILE“ READ ‘.(Y(L).L=6.N2.2) PRINT *.(Y(L).L=6.N2.2) N3=N2-5 Y(1)'AVG(Y(5).N3) Y(2)=AVG(Y(6).N3) ESTIMATION OF THE GAS-DYNAMIC TEMPERATURE THAT THE EXPERIMENTAL ASSUMPTION OF ZERO DISSOCIATION WOULD REPORT UNDER THE CONDITIONS PREDICTED BY THE MODEL TG-Y(2)*5.t(7.+3.tY(1))¢(1.+Y(1))/7./(5.+Y(1)) PRINT 20.2.(Y(L).L-1.4).UB(Y.2).AR(2).VS(Y).TG PRINT 30.(Y(L).L-5.N2.2) PRINT 30.(Y(L).L-6.N2.2) FDRMAT(8X.10611.4) FORMAT(3X.1OG12.5) C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C 101 104 00000 00000 ZF-AMIN1(Z+ZP.ZL) DZ‘ZP Z2'AMIN1(Z+DZ.ZF) CALL RK4(DELTA.Z.Z2.Y.W.N2.ERR) WHEN THE INTEGRATION SUBROUTINE REPORTED THAT IT DIVERGED. THIS VERSION WOULD CUT THE INITIAL STEP SIZE BY A FACTOR OF 16 AND REPEAT. UNTIL IT WAS REDUCED BY A FACTOR OF 10 MILLION IF (22.LT.O) Dz-DZ/16. IF (DZ.LT.1.E-7*ZP) GO TO 103 IF (z.LT.zF) GO TO 104 ESTIMATION OF THE GAS-DYNAMIC TEMPERATURE THAT THE EXPERIMENTAL ASSUMPTION OF ZERO DISSOCIATION WOULD REPORT UNDER THE CONDITIONS PREDICTED BY THE MODEL TG'Y(2)‘5-*(7.+3.‘Y(1))‘(1.+V(1))/7./(5.+Y(1)) PRINT 20.2.(Y(L).L-1.4),UB(Y.Z).AR(z).VS(Y).TG PRINT ao.(Y(L).L-5.N2.2) 185 PRINT 30.(Y(L).L'6.N2.2) IF(Z.LT.ZL)GO TO 101 GO TO 103 102 CONTINUE 000000000 END FUNCTION AR(Z) REACTOR AREA FLOW PROFILE SUBROUTINE PRESENTLY CONSISTS OF A 25 MM. STRAIGHT LENGTH OF TUBE WITH CROSS-SECTIONAL AREA AK (OPERATOR INPUT. DEFAULT 3.8E-4 SQ. M.) FOLLOWED BY A STRAIGHT CONICAL TAPER OVER ANOTHER 25 MM. TO A THROAT WITH AREA A2 (ALSO INPUT. SAME DEFAULT) WHICH IS FOLLOWED BY ANOTHER STRAIGHT CONICAL OIVERGENT TAPER OVER 25 CM. LENGTH BACK TO AREA AK. CDMMON/ALFA/A2.AK DATA 21.22.23.24/0...ozs..os..30/ AR-AK IF(AK.EQ.A2)RETURN IF((z.LE.21).OR.(2.GE.Z4))RETURN 100 IF(Z.LT.Z2)GO TO 101 IF(Z.GT.23)GO TO 102 AR=A2 RETURN 101 ZA=ZT 28322 GO TO 103 102 2A824 28823 103 AR-(SQRT(AK)+(SQRT(AK)-SQRT(A2))*(Z-2A)/(ZA-ZB))**2 000 00000 000 RETURN END FUNCTION DAR(2) COMPUTATION OF DERIVATIVE OF CROSS-SECTIONAL AREA WITH LENGTH USED IN INVISCID FLOW RELATIONS FINITE DIFFERENCE CALCULATION. WITH STEP-SIZE 1 MICRON. DY-1.E-06 DAR-(AR(z+DY)-AR(Z))/DY RETURN END FUNCTION UB(Y.Z) CALCULATION OF BULK VELOCITY BASED ON IDEAL GAS LAW COMMON/BRAVO/AMO.FO.R DIMENSION Y(4) UB=FO#(1.+Y(1))FRFY(2)/Y(3)/AR(Z) RETURN END SUBROUTINE DELTA(Z.N.Y.DY) CALCULATION OF MATERIAL AND ENERGY BALANCES FOR FIRST DERIVATIVES IMPLICIT REAL (K) COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.TO.C(3.3).GAMMA COMMON/ECHO/TAIR.UI COMMON/POPPA/PI DIMENSION A(3.3).Y(N).DY(N) 186 DATA UGO/4.0648E08/ c C RATE CONSTANTS DERIVED FROM SHUI AND APPLETON. SHUI. AND WOOD AND WISE C DATA C LEAST-SQUARES FIT TO FORM OF EXPRESSION K - A'(TF*B)*EXP(C/T) C DATA A/2.014E11.-O.7497.-20.43.B.B16E11.'0.7219.-92.09.3.855E18. +-3..-2900./ RATE AND EQUILIBRIUM CONSTANT EXPRESSIONS 000 K(I,T)-A(1.I)*TFFA(2.I)*EXP(A(3.I)/T) KEQ(T)-1.01325E05*Exp((-DGO/TO+C(1.3)FALOG(T/TO)+(DHO/ +T/TO+C(2.3)+C(3.3)‘(T+TO)/2.)‘(T-TO))/R) DH(T)=DHO+T*CP(3.T) CP(I.T)=C(1.I)+T*(C(2.I)+T#C(3.I)) DH2HP(T)=1.264E-3tTtt1.72 KTH(T)=1.843E-3*T**.8 KINETIC ENERGY OF FLOW AS FUNCTION OF X. T. P. AND Z 000 EK(Y.Z)=1./(1./(AMO*UB(Y.Z)'F2)-1./R/Y(2)/(1.+Y(1))) AS=2.*SORT(PI*AR(Z)) IF(UI)111.112.113 112 HO=O. GO TO 110 111 HO=ABS(UI) GO TO 110 113 HO=UI-47.OB 110 CONTINUE DR=2./(N-6.) N3=N-3 DO 160 U=5.N3.2 R2=DR*(d-4.)/2. R1=AMAX1(R2-DR.O.) C C MATERIAL BALANCE EQUATION INCLUDING RATE EXPRESSION C AND RADIAL DIFFUSION TERM C DXRaAR(Z)/FOF((K(2.Y(J+1))¢2.*Y(d)+K(1.Y(d+1))*(1.-Y(d)))‘ +(KEQ(Y(d+1))‘(1.-Y(J))/Y(3)‘(1.+Y(d))‘4.*Y(d)‘Y(d))‘ +(Y(3)/R/Y(d+1)/(1.+Y(d)))¥*3)/(2.-(R2tR2+R1FR1)) DY(d)-DXR-PI/FO/R/DRF(DH2HP((Y(J+1)+Y(d-1))/2.)FRTF +(Y(u)/(Y(d)+1.)/Y(d+1)-Y(d-2)/(Y(d-2)+1.)/Y(d-1))+ +DH2HP((Y(J+1)+Y(d+3))/2.)‘92‘ +(Y(d)/(Y(d)+1.)/Y(d+1)-Y(d+2)/(Y(d+2)+1.)/Y(d+3)))* +(1.+Y(d))/(1.+2.#Y(d))/(2.-R1tR1-R2tR2) c c ENERGY BALANCE EQUATION INCLUDING HEAT OF REACTION. KINETIC ENERGY C AND THERMAL ENERGY TERMS C DY(d+1)'(((EK(Y.Z)/(1.+Y(d)))*DY(J)+DH(Y(d+1))#DXR-EK(Y.Z)* +DAR(Z)/AR(Z))+2.tPI/DR/FO/(RZFR2-R1RR1)/(2.-R2#R2-R1tR1)* +(KTH((Y(d+1)+Y(J-1))/2.0)*R1‘(Y(d+1)-Y(d-1))+KTH((Y(J+1)+Y(d+3))/ +2.)#(R2*Y(d+1)-R2*Y(d+3))))/(EK(V.Z)/Y(J+1)+ +(1.-Y(d))*CP(1.Y(d+1))+2.#Y(J)FCP(2.Y(J+1)))F(-1.) 160 CONTINUE C c COMPUTATION OF BOUNDARY CONDITIONS AT OUTERMOST STATION C QO=(Y(N)-Y(4))*HO*AS 187 DXR=AR(Z)/FO*(AS/AR(Z)*K(3.Y(N))+(1.-R2*R2)*Y(3)/R/Y(N) +(1.+Y(N-1))F(K(2.Y(N))t2.*Y(N-1)+K(1.Y(N))t(1.-Y(N-1))) +(KEQ(Y(N))*(1.-Y(N-1)"2)/Y(3)-4.*Y(N-1)*Y(N-1))F +t(Y(3)/R/Y(N)/(1.+Y(N-1)))**2/(1.-R2*R2*(2.-R2*R2)) DY(N-1)FDXR-PI/FO*DH2HP((Y(N)+Y(N-2))/2.)/R/DR/(1.+Y(N-1)/ +(Y(N-1)+1. ))*(R2*Y(N-1)/(1.+Y(N-1))/Y(N)-R2*Y(N-3)/ +(1.-Y(N-3))/(N-2))/(1.-R2*R2) t / ) C C ENERGY BALANCE EQUATION INCLUDING HEAT OF REACTION. KINETIC ENERGY C AND THERMAL ENERGY TERMS c DY(N)=(((EK(Y.Z)/(1.+Y(N-1)))tDY(N-1)+DH(Y(N))*DXR-EK(Y.Z)* +DAR(Z)/AR(Z))+(KTH((Y(N)+Y(N-2))/2.)/DRF(R2*Y(N)-R2# +Y(N-2))*2.FPI+QD)/Fo/(1.-R2FR2)**2)/(EK(Y.z)/Y(N)+ +(1.-Y(N-1))FCP(1.Y(N))+2.#Y(N-1)tCP(2.Y(N)))t(-1.) DY(1)-AVG(DY(5).N-5) DY(2)-AVG(DY(6).N—5) DY(3)=-EK(Y.Z)*(DY(2)/Y(2)+DY(1)/(1.+Y(1))-DAR(2)/AR(Z))F +Y(3)/R/Y(2)/(1.+Y(1)) DY(4)=(QO+1.64t(300.-Y(4)))/TAIR/3.5/R RETURN END SUBROUTINE RK4(DELTA.X.XF.Y.W.N.TOL) C C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA RUNGE-KUTTA 4 METHOD C DIMENSION Y(N).V(N.9) DXO=XF-x DXL=1. KOUNT-~1 100 KOUNTsKOUNT+1 XC=.1 XL=O. DD 5 1-1.N 5 H(I.1)=Y(I) 101 CONTINUE CALL DELTA(X+XL‘DXO.N.W(1.1).W(1.2)) C c PROVIDES CHECK AGAINST STIFF EQUATIONS. IF CHANGE OVER INTEGRATION STE c STEP Is TOO LARGE. REDUCES STEP SIZE. IF STEP SIZE REACHES FLOOR C RETURNS WITH FLAG SET TO TELL THAT THE INTEGRATION DIVERGED. C DD 120 K=2.N.2 105 IF(ABS(W(K.2)*DXLFDXO).LE.O.1*W(K.1))60 TO 114 DXL=DXL/2. GO TO 105 114 CONTINUE 12o CONTINUE IF(DXL.GT.2.**(-20))GD TO 110 PRINT .,- INTEGRATION DIVERGED” XFa-XF RETURN 110 XL-XL+DXL/2. X1-X+XL*DXO DO 10 I-1.N 10 W(I.6)-W(I.1)+DXL*DXORW(I.2)/2. CALL DELTA(X1.N.W(1.6).W(1.3)) DO 20 I-1.N 20 V(I.7)=w(I.1)+DXLtDXOtw(I.3)/2. CALL DELTA(X1.N.W(1.7).W(1.4)) DO 30 I=1.N 30 40 1000 1001 51 50 102 60 103 70 C 188 W(I.8)=W(I.1)+DXL*DXO*W(I.4) XL=XL+DXL/2. CALL DELTA(X+XL*DXO.N.W(1.8).W(1.5)) DO 40 1:1.N W(I.1)=W(I.1)+(W(I.2)+2.F(W(I.3)+W(I.4))+W(I.5))*DXL*DXO/6. IF (XL—XC) 1001.1000.1000 XC=XC+1./10. CONTINUE IF(XL.LT.1.)GO TO 101 DXL=DXL/2. IF(KOUNT.EQ.O)GO TO 102 ERR=O. DO so I-1.N IF(ABS((W(1.9)-Y(I))/W(I.9)).LE.1.E-8)GO TO 51 ERR=ERR+((W(I.9)-W(I.1))/(W(I,9)-Y(I)))**2 CONTINUE CONTINUE IF(ERR.LE.TOL)GO TO 103 OD so I-1.N W(l.9)=W(I.1) GO TO 100 X=XF DD 70 I-1,N Y(I)-V(I.9) RETURN END FUNCTION AVG(Y.N) C COMPUTES BULK AVERAGES BASED ON PARABOLIC MASS FLOW PROFILE C 10 DIMENSION Y(N) A1=NT1 N2=N-2 A3=A1/2 AVG=Y(1)/A1**2*(2.-1./A1*‘2)+Y(N)*(1.-(N2/A1)‘*2*(2.'(N2/A1)"2)) DO 10 I33.N2.2 A=I/A1 B-A-1./A3 AVG=AVG+Y(I)#(A*A*(2.-A*A)-B*B*(2.-B‘B)) RETURN END 189 Program 4 Axial Dispersion Model Two-Dimensional, Axial Dispersion, Cooled, Surface Reaction, Inviscid Flow 190 PROGRAM MODEL C C DUMMY MASTER PROGRAM THAT CALLS MAIN SUBROUTINE C CALL ROCKET STOP END 191 SUBROUTINE ROCKET IMPLICIT REAL‘B (A-H.K.O-Z) c c PROGRAM TO MODEL HYDROGEN ATOM RECOMBINATION IN AN AXIALLY DISPERSED- c FLOW REACTOR C Two-DIMENSIONAL MODEL WITH IMSL DGEAR INTEGRATOR c SEARCH TO FIND INITIAL DERIVATIVES FOR DISPERSION MODEL C MODIFICATION BEGUN 5-27-86 C MODIFIED INTO SUBROUTINE TO PERMIT ERROR RECOVERY TRANSFER VIA CALL OF c SECOND ENTRY POINT c COMMON/ALFA/A2,AK COMMON/BRAVO/AMO,FO.R COMMON/CHARLIE/DHO.DGO.TO.C(3.3) COMMON/ECHQIERR.TS.xs.TAIR.Uo COMMON/GOLF/GAMMA COMMON/HOTEL/DISP(2) COMMON/POPPA/PI COMMON /INDIA/IF1.IF2 EXTERNAL DELTA2.DUMMY DIMENSION Y(47).w(47.75).IwK(47) DIMENSION VO(47).Y1(47).DZY(47).03V(47) DIMENSION SRC(2O.3).IS(20).IB(2O),HMAT(2O.20) DATA R.TO/8314.3400.298.1SDO/ DATA C/26880.DO.4.34700.-2.6450-04.20786.DO.O.DD.D.DO. +0.D0.0.D0.0.DO/ DATA DHO.DGO/4.3141008.4.0648008/ DATA AMO/2.DO/ DH(T)=DHO+T'CP(3.T) KEQ(T)=1.01325005‘DEXP((’DGO/TO+C(1.3)‘DLOG(T/TO)+(DHO/ +T/TO+C(2.3)+C(3.3)‘(T+TO)/2.DO)‘(T-TO))IR) XEQ(T.P)=1.DOIDSQRT(1.DO+4.DO‘P/KEQ(T)) VMACH(Y.Z)=UB(Y.Z)/VS(Y) DH2HP(T)=1.2640-3‘T“1.7ZDO KTH(T)=1.843D-3‘T°'.BDO c C INITIALIZATION INCLUDING EVALUATION OF DELTA CP OF REACTION C AND STANDARD SET OF INITIAL CONDITIONS AND SYSTEM PARAMETERS C OPEN (UNIT=34.FILE='STARTG.DAT'.STATUS='OLD') OPEN (UNIT=44.FILE=’SYSSOUTPUT'.STATUS=‘NEW’) OPEN (UNIT=45.FILE=‘SAVES.DAT'.STATUS='OLD’) OPEN (UNIT=SS.FILE=‘SEARCH.DAT'.STATUS='OLD') OPEN (UNIT865.FILE=’HMAT.DAT'.STATUS='OLD') INPR=3 ZQ=-1.DO INPRO=INPR DO 5 I=1.3 5 C(I.3)=(2.DO’C(I.2)-C(I.1))II Y(1)=O.DO Y(2)=O.1DO V(3)=1200.DO V(4)=3530.DO Y(5)=400.DO FO=3.BZD-7 ERR=1.D-6 AK=3.BD-4 A2=3.BD-4 TAIR=1.7D-4 PI=4.DO‘DATAN(1.DOO) ICON=O UO=1.GDO DO 555 I=1.47 DO 555 J=1.75 555 W(I.J)=O.DO 103 CONTINUE 192 WRITE (44.3).' ENTER INITIAL CONDITIONS AND SYSTEM PARAMETERS' READ (34,.).(Y(L).L=1.5).FO.ERR.AK.A2.ZL.ZP.TAIR.UO.INPRO.N1. +METH.MITER.DISP.GAMMA IF(Y(1).LT.O.DO)GO TO 102 WRITE (44.‘).(V(L).L=1.5).FO.ERR.AK.A2.ZL.ZP.TAIR.UO.INPRO.N1. +METH.MITER.DISP.GAMMA N2=N1‘2+5 N7=29N1 N3=4‘N1+7 GO TO (111.112.113.111.115).INPRO C C IF INPROtI. INITIAL PROFILE IS UNIFORM ACROSS REACTOR C 111 DO 110 I=1.N1 Y(I+5)=v(2) 110 Y(I+N1+5)=v(3) IF(INPRO.NE.4)GO To 150 C C IF INPRO=4. UNIFORM PROFILE EXCEPT THAT OUTERMOST STATION IS AT C COOLANT TEMPERATURE AND ZERO DISSOCIATION C V(N2-N1)=O.DO Y(N2)=Y(5) GO To 150 c C IF INPRO=2. INITIAL PROFILE Is UNIFORM AT ZERO DISSOCIATION AND c COOLANT TEMPERATURE EXCEPT FOR CENTRAL SEGMENT. WHICH IS AT INPUT C CONDITIONS (Y(2) AND Y(3)). C 1 12 Y(6)=V(2) Y(6+N1)=Y(3) DO 120 I=2.N1 v(I+5)=O.DO 120 Y(I+N1+5)=v(5) GO TO 150 C C IF INPRO = 5. INITIAL PROFILE MUST BE INPUT. DERIVATIVES CALCULATED C FROM DANCKWERTS BOUNDARY CONDITIONS BASED ON UNIFORM TO AND X0. C 115 XO=Y(2) TOD-Y(3) C C IF INPRO = 3. INITIAL PROFILES MUST BE ENTERED AT PROGRAM PROMPTS. C 113 WRITE (44.‘).‘ X PROFILE' READ (34.‘).(V(5+L).L=1.N1) WRITE (44.‘).(V(5+L).L=1.N1) WRITE (44.‘).' T PROFILE' READ (34.9).(Y(L+N1+5).L=1.N1) WRITE (44.‘).(Y(L+N1+5).L=1.N1) 50 CONTINUE Y(2) AND Y(3) ARE THE BULK (MASS AVERAGE) DISSOCIATION AND TEMPERATURE RESPECTIVELV. EXCEPT WHEN INPRO=1. THEY MUST BE COMPUTED INITIALLV FROM THE ACTUAL INITIAL PROFILES 00000-fi Y(2)=AVG(Y(6).N1) Y(3)=AVG(V(6*N1).N1) READ IN SEARCH PARAMETERS AND CONTROL VALUES 000 READ (55.‘).(((SRC(I.J).J=1.3).IS(I).IB(I)).I=1.N7) READ (55.‘).ICON READS IN MATRIX FOR TRANSFORMED SEARCH COORDINATES DEFAULT IS UNTRANSFORMED 000 931 932 902 901 1001 000 0000 000 #00000 19 701 30 20 193 00 932 I=1.N7 00 931 J=1.N7 HMAT(I.J)=0.DO HMAT(I.I)=1.00 READ (65.‘).((HMAT(I.J).J=1.N7).I=1.N7) 00 901 I=1.N7 Y(I+N2+1)=0.Do Do 902 J=1.N7 Y(I+N2+1)=Y(I*N2+1)+HMAT(I.J)‘SRC(J.3) CONTINUE Y(N2+I)=AVG(V(N2+2).N1) Y(N3)8AVG(Y(N3-N1).N1) DO 1001 I=1.N3 Y0(I)=Y(I) TS=Y(N2) XS=V(N2-N1) IF1=0 IF2=0 ENTRY POINT FOR ERROR RECOVERY RESTART ENTRY EGRESS IF TWO ERRORS WITHOUT REACHING NEW PRINT POINT. RESET T0 INPUT START POSITION IF(IF1.GE.2) GO TO 401 READ IN SAVED LAST PRINT POSITION FOR RESTART REWIND 45 READ (45.‘).(Y(I).I=1.N3).TS.XS ESTIMATION OF THE GAS-DYNAMIC TEMPERATURE THAT THE EXPERIMENTAL ASSUMPTION 0F ZERO DISSOCIATION WOULD REPORT UNDER THE CONDITIONS PREDICTED BY THE MODEL TG=V(3)‘5.00‘(7.00+3.00'V(2))‘(1.00+Y(2))/7.00/(5.DO+Y(2)) WRITE (44.‘).(Y0(6+2‘N1+L).L=1.N1) WRITE (44.¢).(YO(L+30N1+6).L=1.N1) IF2=1 CALL DELTA2(N3.0.00.Y.02Y) IF(IF1.GT.O)THEN TS=Y(N2) XS=Y(N2-N1) IF1=0 GO TO 119 ENDIF IF2=O DO 701 I=2.N3 02Y(I)=02V(I)/02Y(1) 02V(1)=1.00 220=Y(1) VB=UB(Y(2).Y(1)) VSS=VS(Y(2)) IF (ICON.EQ.O)GO TO 1200 WRITE (44.20).(Y(L).L=1.5),VB.AR(Y(1)).VSS.TG WRITE (44.30).(Y(5+L).L=1.N1).XS WRITE (44.30).(02Y(5+L).L=1.N1) IF(DISP(1).NE.O.DO)wRITE (44.30).(02Y(6+2‘N1+L).L=1.N1) WRITE (44.30).(Y(L+N1+5).L=1,NT),Ts WRITE (44.30).(02Y(L+N1+5).L=1.N1) IF(DISP(2).NE.0.DO)WRITE (44.30).(02V(L+3‘N1+6).L=1.N1) FORMAT(BX.TOGTT.A) FORMAT(3X.SG12.5) 194 C C INTEGRATION OF DIFFERENTIAL EQUATIONS VIA IMSL DGEAR SUBROUTINE C 1200 Z=O.DO ZPR=220 INDEX=1 101 ZPR=DMIN1(ZPR+ZP.ZL) vo=v(1) DO 601 I=1.N3 601 v1(I)=DZV(I) 105 DZ=ZP C C IN ARC-LENGTH INTEGRATION. USES VARIABLE STEP SIZE AS INPUT T0 C INTEGRATOR TO CONSERVE COMPUTATIONAL TIME (MIGHT NEGATE ADVANTAGES C 0F ARC-LENGTH METHOD) C IF(VO.NE.Y(1))DZ=DMAX1(DABS((ZPR-Y(1))/(V(1)-VO))‘DZ/16.00.ZP) YO=V(1) on=Dz 104 22=Z+Dz H=DZ CALL DGEAR(N3.DELTA2.DUMMY.Z.H.Y.22.ERR,METH.MITER.INDEX.IWK.W. +IER) IF (IER.LE.127) GO To 106 WRITE (44.0) ’INTEGRATION DIVERGED' GO To 103 106 IF (V(1).LT.ZPR) GO TO 105 CALL DELTA2(N3.Z.Y.DZY) DO 702 I=2.N3 702 02Y(I)=02Y(I)/02Y(1) 02v(1)=1.00 DO 602 I=1.N3 602 DBY(I)-(DZY(I)-Y1(I))/ZP TG=Y(3)‘5.DO‘(7.DO+3.DO‘Y(2))'(1.00+Y(2))/7.00/(5.00+V(2)) IF1=O IF (ICON.NE.O)THEN WRITE (44.20).(V(L).L=1.5).U8(V(2).Y(1)).AR(V(1)). +VS(Y(2)).TG WRITE (44.30).(Y(L+5).L=1.N1).Xs WRITE (44.30).(DZY(5+1),L=1.N1) IF(DISP(1).NE.0.DO)WRITE (44.30).(02Y(N2+1+1).L=1.N1) WRITE (44.30).(Y(L+N1+5).L=1.NT).Ts WRITE (44.30).(02Y(L+N1+5).L=1.N1) IF(DISP(2).NE.0.DO)WRITE‘(44.30).(02Y(L+3‘N1+6).L=1.N1) ENDIF SAVE DATA FOR RESTART 000 REWIND 45 WRITE (45.‘).(Y(L).L=1.N3).TS.XS SERIES OF CHECKS FOR RUNAWAY CONDITIONS USABLE IN SEARCH TECHNIQUE TO ESTABLISH CORRECT VALUES FOR INITIAL DERIVATIVES 00000 DO 595 IN=6.N2 IN1=IN-5 IN2=IN+2‘N1*1 IF PLUG-FLOW 0N EITHER MASS 0R ENERGY. DOES NOT SEARCH 000 ID=1 YMIN=O.DO IF(Y(IN).GT.XS)YMIN=XS IF(IN1.GT.N1)ID=2 IF(IN1.GT.N1)YMIN=Y(5) 195 IF(DISP(ID).NE.O.DO)THEN c c DIRECT CHECK FOR NEGATIVE RUNAWAY c IF(V(IN)+ZP'Y(IN2).LT.YMIN)GO T0 123 c c INDIRECT CHECK FOR POSITIVE RUNAWAY BASED ON FIRST. SECOND. AND c THIRD DERIVATIVES BEING POSITIVE IMPLYING RUNAWAY C IF(Y(IN2).GT.O.DO.AND.D3Y(IN2).GT.O.DO.AND. +DZY(IN2).GT.O.DO)GO TO 133 ENDIF 59s CONTINUE IF (Y(1).LT.ZL) GO TO 101 IF (ICON.EQ.1)GO T0 103 ICON=1 Go To 400 C C SEARCH ALGORITHM TO FIND NEXT GUESS FOR INITIAL DERIVATIVES C 123 ISN=IS(IN1) DO 125 I=1.N7 125 IS(I)=IB(I) IS(IN1)=1 IF(SRC(IN1.1).EQ.SRC(IN1.2))THEN SRC(IN1.1)=SRC(IN1,1)-DMIN1(.O1D0..O1DO‘DABS(SRC(IN1.1))) SRC(IN1.2)=SRC(IN1.2)+DMIN1(.0100..01DO'DABS(SRC(IN1.2))) ENDIF IF(SRC(IN1.3).EQ.SRC(IN1.2))GO T0 126 IF(ISN.NE.0)GO T0 124 SRC(IN1.1)=SRC(IN1.3) SRC(IN1.3)=SRC(IN1.2) GO TO 400 124 SRC(IN1.1)=SRC(IN1.3) 154 SRC(IN1.3)=(SRC(IN1.1)+SRC(IN1.2))/2.DO GO TO 400 126 SRC(IN1.3)=2.DO‘SRC(IN1.2)-SRC(IN1.1) SRC(IN1,1)=SRC(IN1.2) SRC(IN1.2)=SRC(IN1.3) Go To 400 133 ISN=IS(IN1) 00 135 I=1.N7 135 IS(I)=IB(I) IS(IN1)=1 IF(SRC(IN1.1).EQ.SRC(INJ.2))THEN SRC(IN1.1)=SRC(IN1.1)-DMIN1(.O1DO..01DO‘DABS(SRC(IN1.1))) SRC(IN1.2)=SRC(IN1.2)+DMIN1(.01DO..01DO'DABS(SRC(IN1.2))) ENDIF IF(SRC(IN1.3).EO.SRC(IN1.1))GO TO 136 IF(ISN.NE.0)GO To 134 SRC(IN1.2)=SRC(IN1.3) SRC(IN1.3)=SRC(IN1.1) GO TO 400 134 SRC(IN1.2)=SRC(IN1.3) GO TO 154 136 SRC(IN1.3)=2.DO'SRC(IN1.1)-SRC(IN1.2) SRC(IN1.2)=SRC(IN1.1) SRC(IN1.1)=SRC(IN1.3) 400 REWIND 55 WRITE (55.‘).(((SRC(I.J).J=1.3).IS(I).IB(I)).I=1,N7) WRITE (55.‘).ICON DO 911 I=1.N7 YO(I+N2+1)=D.DO DO 912 J=1.N7 912 Y0(I+N2+1)=Y0(I+N2+1)+HMAT(I.J)'SRC(J.3) 911 CONTINUE 401 1002 102 196 v0(N2+I)=AVG(YO(N2+2).N1) VO(N3)=AVG(V0(N3-N1).N1) DO 1002 I=1.N3 Y(I)=V0(I) GO To 119 CONTINUE IF(IF1.EQ.0)STOP RETURN END 0000 0000 140 142 143 141 150 152 153 151 100 131 130 C 197 SUBROUTINE DELTA2(N.Z.Y.DY) IMPLICIT REAL‘B (A-H.K.O-Z) COMMON/BRAVOIAMO.FO.R COMMON/CHARLIE/DHO.DGO.TO.C(3.3) COMMON/HOTEL/DISPM.EDISP COMMON /INDIA/IF1.IF2 DIMENSION Y(N).DY(N) DH2HP(T)=1.2640-3‘T“1.7200 KTH(T)=1.843D-3‘T“.ODO N2=(N-7)/4 NF=N-2‘N2-3 A1=N2 CALL DELTA (NF.Y(1).Y(2).DY(2)) PROVISION FOR CASCADING ERROR RECOVERY RETURN IF USING OWN INTEGRATOR SUBROUTINE IF(IF1‘IF2.NE.O)RETURN IF DISPERSION COEFFICIENTS SET TO 0.. CIRCUMVENT DISPERSION CALCULATION TO AVOID INFINITIES IF(EDISP.EQ.O.DO.AND.DISPM.EQ.O.DDD)GO TO 130 R2=O.DO DO 100 I=1.N2 I1=I+5 12=I+N2+5 IG=NF+N2+2+I A2=I R1=R2 R2=DMIN1(1.DO.A2/A1) IF(EOISP)140.141.142 KTHY=DABS(EDISP) GO TO 143 KTHY=EDISP‘KTH(Y(IZ)) DY(13)=(Y(13)-DY(12))‘(EK(Y(2).Y(1))/Y(12)+CPM(Y(II).Y(12)))‘FO’ +(2.DO-R2‘R2-R1‘R1)IAR(Y(1))/KTHY DY(IZ)=V(13) 14=NF+2+I IF(OISPM)150.151.152 DHH2P=DABS(DISPM) GO TO 153 DHH2P=DISPM‘DH2HP(Y(I2)) DY(IA)=Y(I1)‘(V(I1)+1.DO)‘((Y(I4)-DY(I1))‘R‘Y(12)‘(Y(I1)+1.00)‘ +F0'(2.DO-R2‘R2-R1'R1)/Y(11)IAR(V(1))IDHH2P+DY(13)/Y(12)+ +2.00'(Y(14)"2/(Y(I1)+1.DO)“2/Y(IT)+DY(4)‘DY(12)/Y(4)/Y(IZ)+ +Y(I4)'DY(I2)/Y(I2)/Y(I1)/(Y(I1)+1.DO)-DY(4)‘Y(I4)/Y(4)/V(I1)l +(Y(I1)+1.DO))) DY(II)=Y(I4) CONTINUE CONTINUE IF(EDISP.EQ.0.DO)GO T0 131 DY(N)=AVG(DY(N-N2).N2) DY(3)=V(N) IF(DISPM.EQ.0.DO)GO TO 130 DY(N-2'N2-1)=AVG(DY(N-N2‘2).N2) DY(2)=Y(N-2'N2-1) Dv(1)=1.00 C COMPUTATION OF TERMS FOR ARC-LENGTH METHOD INTEGRATION C 110 DS=0.DO DO 110 I=1.N DS=DS+DY(I)‘DY(I) DS=SQRT(DS) DO 120 I=1.N 198 120 DY(I)=DY(I)/Ds RETURN END SUBROUTINE DUMMY(N.X.Y.PD) C C DUMMY SUBROUTINE REQUIRED FOR IMSL DGEAR COMPUTATION C DIMENSION Y(N).PD(N.N) RETURN END 199 SUBROUTINE DELTA(N.Z.Y.DY) COMPUTES FIRST’ORDER DERIVATIVES 0F TEMPERATURE. CONCENTRATION. COOLANT TEMPERATURE. AND PRESSURE FROM MATERIAL AND ENERGY BALANCES. 00000 IMPLICIT REAL‘B (A-H.K.O-Z) COMMON/BRAVO/AMO.FO.R COMMON/CHARLIE/DHO.DGO.TO.C(3.3) COMMON/ECHOIERR.TS.XS.TAIR.UI COMMON/GOLF/GAMMA COMMON/POPPA/PI COMMON /INDIA/IF1.IF2 DIMENSION A(3,3).Y(N).DY(N) DIMENSION XSN(2).TSN(2).XSX(2).TSX(2).DXN(2) DIMENSION XTN(2.4).FXT(4).SUM(2) RATE CONSTANTS DERIVED FROM SHUI. APPLETON AND SHUI. AND SANCIER AND WISE DATA. LEAST-SQUARES FIT TO FORM OF A‘T“B‘EXP(C/T) 0000 DATA A/2.014011.‘O.7497DO.'20.43DO.8.816011.-O.721QDO.-92.0QDO. +3.855018.’3.DO.’2900.DO/ C C RATE AND EQUILIBRIUM CONSTANT EXPRESSIONS C K(I,T)=A(1.I)‘T"A(Z.I)‘DEXP(A(3.I)/T) KEQ(T)=1.O1325005‘DEXP(('DGO/T0+C(1.3)¢DL05(T/To)+(DHo/ +T/TO+C(2.3)*C(3.3)‘(T+T0)/2.DO)‘(T-TO))IR) XEQ(T.P)=1.DO/DSQRT(1.DO+4.DO‘P/KEQ(T)) DH(T)=DHO+T’CP(3.T) DHZHP(T)=1.2640-3‘T“1.7ZDO KTH(T)=1.343D-3‘T“.BDO DEBUGGING PRINT STATEMENT 000 IF(GAMMA.NE.D.DO)WRITE (44,9),z.v GZ=GAMMA AS=2.DO‘DSQRT(PI‘AR(Z)) IF(UI)111,112.113 112 HO=O.DO GO TO 110 111 HO=DABS(UI) GO TO 110 113 H0=UI‘47.0600 110 CONTINUE DR=2.00/(N-4.00) N2=(N-4)/2 N3=N2+5 N4=N2-1 CENTRAL BOUNDARY CONDITIONS 000 R2=O.DO RQZ=O.DO RNH2=0.DO L5=1 NCOUNT=1 LCOUNT=O C C OUTERMOST STATION WITH WALL BOUNDARY CONDITIONS COMPUTED SEPARATELY C SINGLE STATION COMPUTATION PSEUDO'ONE-DIMENSIONAL C IF(N4.EQ.O)GO T0 170 00 160 J=1.N4 THIS SAVES COMPUTATION AND ENSURES CONSERVATION 00 0000 00000 000 000 d000000‘ O) O ‘1 O 200 R18R2 RQ1=RQZ RNH1=RNH2 R2=DR‘J J1=J+4 J2=J1+N2 ERROR BAILOUT. IF OWN INTEGRATOR SUBROUTINE. CASCADING RETURNS. IF LIBRARY INTEGRATOR. CALL TO RECOVERY ENTRY POINT IN MAIN SUBROUTINE. IF(Y(J2).LE.0.00.0R.Y(J2+1).LE.0.00)THEN IF1I1+IF1 IF(IF2.NE.0)RETURN WRITE (44.-).'ILLEGAL TEMPERATURE' CALL EGRESS ENDIF RADIAL FLUXES BASED ON TRANSPORT PROPERTIES COMPUTED AT AVERAGE CONDITIONS TB-(Y(J2)+Y(J2+1))/2.DO XB=(Y(J1)+Y(J1+1))/2.DO RQZ=R2¢(-KTH(TB)¢(Y(J2+1)—V(J2))/DR) RNH2=R2‘(-DH2HP(T8)‘(Y(J1+1)/(V(J1+1)+1.DO)/Y(J2+1)-Y(J1)/ +(Y(J1)+1)/Y(J2))'2.DO/RIDR)‘(1.DO-XB/(1.D0+XB)) MATERIAL BALANCE BASED ON AXIAL FLOW AND REACTION ALONE DXR=AR(Z)/FO‘((K(2.Y(J2))‘2.DO‘Y(J1)+K(1.Y(J2))‘(1.D0~Y(J1)) ‘ ) +(KEQ(Y(J2))'(1.DO-Y(J1))/Y(3)‘(1.DO+Y(J1))'4.DO‘Y(J1)‘Y(J1))‘ +(Y(3)/R/Y(J2)l(1.DO+Y(J1)))“3)/(2.DO'(R2‘R2+R1‘R1)) MATERIAL BALANCE INCLUDING RADIAL MASS FLUX DY(J1)=DXR-2.DO‘PI/F0'(RNH2-RNH1)/(2.DO-R2‘R2-R1‘R1) DY(J2)=(((EK(Y,Z)/(1.D0+Y(J1)))‘DY(J1)+DH(Y(J2))‘DXR-EK(Y.Z)‘ +0AR(Z)/AR(Z))+2.DO‘PI/F0l(R2‘R2-R1‘R1)/(2.DO-R2‘R2-R1'R1)‘ +(R02-RO1))/(EK(Y.Z)/Y(J2)+ +(1.DO-Y(J1))'CP(1.Y(J2))+2.DO‘Y(J1)‘CP(2.Y(J2)))'(-1.00) CONTINUE SOLUTION OF WALL BOUNDARY CONDITIONS BY (1) BOUNDED LINEAR INTER— AND EXTRAPOLATIONS (2) PROJECTING ALONG LINEARIZED ZERO CURVE OF ONE BOUNDARY CONDITION TO FIND NEW POINT TO RESTART INTER- AND EXTRA- POLATIONS (3) SIMPLEX METHOD IF PROJECTION TAKES PAST BOUNDARY LX=1 LT=1 TSMIN=Y(4)‘0.1DO L7=1 J5=1 LAST=4 J6=1 L5=1 150 XWI=XS 120 Ts=DMIN1(Ts.Y(N)+1OODo.DO) 00000 TS 8 DMAX1(TS.TSMIN) ERROR BAILOUT. IF OWN INTEGRATOR SUBROUTINE. CASCADING RETURNS. IF LIBRARY INTEGRATOR. CALL TO RECOVERY ENTRY POINT IN MAIN SUBROUTINE. IF(Y(N).LE.O.D0.0R.TS.LE.0.DO)THEN IF1=1+IF1 130 201 IF(IF2.NE.O)RETURN WRITE (44.‘).’ILLEGAL TEMPERATURE’ CALL EGRESS ENDIF ERRX=ERR‘DMAX1(ERR.XS) HI=4.DO‘KTH((TS+Y(N))/2.DO)‘PI/AS/DR KW=4.DO‘DH2HP((TS*Y(N))l2.DO)‘PI/Y(3)/AS/DR TS1=(HO'Y(4)+HI'Y(N))/(HO*HI) TSI‘TS DXW=AS/(1.DO-RZ‘RZ)“2‘K(3.TS)‘(KEQ(TS)/Y(3)‘(1.DO-XS'XS)- +4.00‘XS'XS)‘(Y(3)/R/TS/(1.DO+XS))"2 DHSSDXW‘DH(TS)'(1.DO-R2‘R2)“2 DXS=(Y(N-N2)/Y(N)I(1.DO+Y(N-N2))+ +DXW‘(1.DO-R2'R2)“2/KW'(1.DO-XS/(1.DO+XS))/AS‘R/Y(3))‘TS XS1=DXS/(1.DO-st) DXSsXST-XWI XS1-DMAX1(XS1.O.DO) XS1=DMIN1(XST.1.DD) TSA=TS1-DHS/AS/(HO+HI) DTS=TSA-TSI IF(GZ.GT.2.DO)WRITE (44.‘).TSI.DTS.XS.DXS.' L5=’.L5.' L7=’.L7. +' LB='.LB IF(NCOUNT.LT.1000)GO TO 300 G2BGAMMA+1.DO NCOUNT=O C ERROR BAILOUT. IF OWN INTEGRATOR SUBROUTINE. CASCADING RETURNS. C IF LIBRARY INTEGRATOR. CALL TO RECOVERY ENTRY POINT IN MAIN C SUBROUTINE. 310 300 134 IF (LCOUNT.LT.10)GO T0 310 IF1=1+IF1 IF(IF2.NE.O)RETURN WRITE (44.‘).'BOUNDARV CONDITIONS HUNG UP' CALL EGRESS LCOUNT=LCOUNT+1 NCOUNT=NCOUNT+1 IF(DABS(DTS).LT.ERR'TS.AND.DABS(DXS).LE.ERRX) GO TO 151 IF(L7.NE.1)GO TO 250 IF(L5.NE.2)GD TO 134 IF(DABS(DXS).GT.ERRX)GO To 131 IF(JS.GT.2)GO TO 134 IF(Gz.GE.2.D0)WRITE (44.‘).'J5='.J5.TS.XS xsx(J5)=xs TSX(J5)=TSI XTN(1.2‘J5)=XS XTN(2.2‘J5)=TSI FXT(Z‘JS-T)=(DTS/Y(N))“2 IF(V(N-N2).NE.0.DO)FXT(2‘JS-1)=FXT(2*J5-1)+(DXS/Y(N-N2))“2 J5=J5+1 IF(J6.GE.3)GO TO 180 CONTINUE IF(L5.EQ.2)LT=1 L5=1 IF(LT.EQ.1)GO TO 124 IF(LT.EQ.2)GO T0 122 IF(DTS‘DTSO.LT.0.DO)GO T0 122 IF(DTS‘DT52.LT.0.DO)GO T0 122 IF(DABS(DTS).LT.DABS(DTSO))GO TO 122 IF(DABS(DTS).LT.DABS(DT52))GD TO 122 C C IF TEMPERATURE BOUNDARY CONDITION GENERATES MINIMUM WITHOUT CROSSING C ZERO. FITS QUADRATIC CURVE AND FINDS MINIMUM FOR NEW VALUE C AO=((DTSZ-DTS)/(T52-TSI)‘(DTSO-DTS)/(TSO-TSI))/(T52TTSO) BQ=(DTSO-DTS)/(TSO-TSI)‘AQ‘(TSO+TSI) 122 131 132 121 124 125 133 123 C 202 TS=—BQ/AQ/2.DO IF(DABS(TS-TSI).LT.ERR.0R.DABS(TS-TSO).LT.ERR.0R.DABS(TS-T52).LT. +ERR)GO TO 131 Go To 133 IF(TSI.EQ.TSO) GO TO 131 IF(DTS.EQ.DTSO) GO TO 131 IF(DABS(DXS).LT.ERRX)GO TO 131 TS=TSI+DTS-(TSO-TSI)/(DTS-DTSO) IF(G2.GE.2.OO)WRITE (44,').TS.XS IF(TS.LT.TSMIN.ANO.(TSI-TSMIN)‘(TSO-TSMIN).E0.0.00)G0 T0 131 IF(L5.EQ.1.AND.DABS(DTS).GT.ERR'TS)GO T0 133 IF(J6.GT.2)GO To 131 IF(G2.GE.2.DO)WRITE (44.¢).'JB-'.JB.TS.XS XSN(JB)-Xs TSN(J6)=TSI XTN(1.20J6-1)=xs XTN(2.2‘J6-1)=TSI FXT(Z‘J6-1)=(DTS/Y(N))“2 IF(Y(N-N2).NE.0.DO)FXT(2‘J6-1)=FXT(2‘J6-1)+(DXS/Y(N-N2))0‘2 DXN(J6)=DXS J6=J6+1 IF(J6.GE.3)GO T0 190 L5=2 TS=TSI IF(LX.EQ.1)GO T0 121 IF(DXs.Eo.DXSO)GO T0 132 IF(DABS(DXS-DXO).LT.DABS(DXS‘(XWO-XWI)))THEN xs=O.DO IF(DXS‘(XWO-XWI)‘(DXS-DXSO).GT.0.DO)XS=1.DO ELSE XS=DMAX1(0.DO.XWI+DXS'(XWO-XWI)/(DXS-DXSO)) ENDIF IF(XS.E0.0.DO.AND.XWO'XWI.EQ.0.00)GO TO 132 XS=DMIN1(1.DO.XS) LX=2 IF(XWO-XS)121.132.121 XS=X$1 LT=1 L5=1 GO TO 125 IF(LX.EQ.1)XS=X51 GO TO 125 LT=2 TS=TSA IF(LX.NE.1)GO T0 123 XWO=XWI DXSO=DXS IF(LX.NE.1)GO To 150 LX=2 LT=3 Tsz=Tso DT52=DTSO TSO=TSI DTSO=DTS GO TO 150 C PROJECTS TO NEW START POINT BY TAKING ONE LINEARIZED ZERO LINE AND C LINEAR EXTRAPOLATION OF OTHER FUNCTION ALONG LINE C 180 181 CONTINUE IF(G2.GE.2.DO)WRITE (44.‘).'PROJECTING’.TSX.XSX.TSN.XSN IF(XSN(2).EQ.XSN(1))GO TO 170 IF(DXN(1).NE.DXN(2))XS=XSN(1)+DXN(1)‘(XSN(2)-XSN(1))/ *(DXN(1)-DXN(2)) TS=TSN(2)+(TSN(2)-TSN(1))‘(XS-XSN(2))/(XSN(2)-XSN(1)) IF(GZ.GE.2.DO)WRITE (44.‘).'PROJECTION'.TS.XS 203 IF(XS.LT.O.DO.AND.TS.LT.TSMIN)GO To 251 XS=DMAX1(XS.O.DO) XS=DMIN1(XS.1.DO) L5=1 IF(XS.EQ.O.DO)L5=2 GO To 170 C c FINDS MINIMUM OF SUM OF SQUARES OF EQUATIONS BY SIMPLEX METHOD c . 250 FXT(LAST)=(DTS/Y(N))“2 IF(Y(N-N2).NE.0.DO)FXT(LAST)=FXT(LAST)+(DXS/Y(N-N2))‘*2 IF(LB.NE.1)GO TO 251 LAST=MINO(4.LAST+1) XSIXTN(1.LAST) TS=XTN(2.LAST) IF(LAST.LE.3)GO To 150 251 L7=2 LB=2 IMIN=1 IMAX=1 DO 252 L=2.3 IF(FXT(L).GT.FXT(IMAX))IMAX=L IF(FXT(L).LT.FXT(IMIN))IMIN=L 252 CONTINUE IF(IMAX.NE.LAST.AND.IMAX.NE.IMIN)GO TO 260 OD 525 L=1.3 Do 525 J=1.2 525 XTN(J.L)=(XTN(J.L)+XTN(J.IMIN))/2.DO LAST=1 LB=1 GO TO 290 260 CONTINUE IF(Gz.GE.2.DO)WRITE (44.‘).XTN.FXT DO 261 L1=1.2 SUM(L1)=O.DO Do 262 L2=1.3 262 SUM(L1)=SUM(L1)+XTN(L1.Lz) 261 CONTINUE LAST=IMAX Do 263 L3=1.2 263 XTN(L3.IMAX)=SUM(L3)-2.DO‘XTN(L3.IMAX) IF(DABS(SUM(1)-3.DO‘XTN(1.IMAX)).LT.ERR.AND. +0ABS(SUM(2)-3.DO‘XTN(2.IMAX)).LT.ERR‘Y(N))GO T0 291 XTN(2.IMAX)=DMIN1(XTN(2.IMAX).Y(N)+TDOOO.DO) XTN(2.IMAX)=DMAX1(XTN(2.IMAX).TSMIN) XTN(1.IMAX)=DMIN1(XTN(1.IMAX).1.DO) XTN(1.IMAX)=DMAX1(XTN(1.IMAX).0.D0) XS=XTN(1.LAST) TS=XTN(2.LAST) 280 IF(DABS(XTN(2.1)-XTN(2.2)).GT.ERR‘Y(N).OR. +DABS(XTN(2.1)-XTN(2.3)).GT.ERR‘Y(N).OR. +DABS(XTN(1.1)-XTN(1.2)).GT.ERR.OR.DABS(XTN(1,1)-XTN(1.3)).GT.ERR) +GO TO 150 281 XS=XTN(1.IMIN) TS=XTN(2.IMIN) GO TO 170 151 QI=(Y(N)-TS)‘HI‘AS QO=QI-DHS c c MATERIAL BALANCE BASED ON AXIAL FLOW AND GAS-PHASE REACION ONLY c DXR=AR(Z)/F0‘(K(2.Y(N))‘2.DO‘Y(N-N2)+K(1.Y(N))‘(1.DO-Y(N-N2)))‘ +(KEO(V(N))‘(1.00-Y(N-N2))/Y(3)‘(1.DO+Y(N-N2))‘4.DO‘Y(N-N2)"2)‘ +(Y(3)/R/Y(N)/(1.DO+Y(N-N2)))"3/(1.DD-RZ‘R2) C c MATERIAL BALANCE INCLUDING RADIAL FLUX AND SURFACE REACTION 204 C DY(N-N2)=DXR+(DXW+2.DO‘PI‘(RNH2)/(1.DO-RZ‘R2))/FO c C ENERGY BALANCE EQUATION INCLUDING HEAT OF REACTION. KINETIC ENERGY c AND THERMAL ENERGY TERMS c DY(N)=(((EK(Y.Z)/(1.D0+Y(N-N2)))‘DY(N-N2)+DH(Y(N))‘DXR-EK(Y.Z)‘ +0AR(Z)/AR(Z))+(-R02‘2.DO‘PITQI)IF0/(1-R2‘R2)“2)/(EK(Y.Z)/Y(N)+ +(1.00-Y(N-N2))'CP(1.Y(N))+2.DO‘Y(N-N2)‘CP(2.Y(N)))‘(-1.00) DY(1)-AVG(DYIE).N2) DY(2)=AVG(DY(N3).N2) C c MECHANICAL ENERGY BALANCE FOR NON-VISCOUS FLOW c DY(3)=-EK(V.Z)'(DY(2)/V(2)+DV(1)/(1.DO+V(1))-DAR(Z)/AR(Z))' +v(3)/R/V(2)/(1.DO+Y(1)) c c ENERGY BALANCE ON COOLANT. INCLUDING RADIATIVE AND CONVECTIVE LOSSES c TO ENVIRONMENT c DY(4)=(OO+1.6400*(300.00-Y(4)))/TAIR/3.500/R IF(GAMMA.NE.O.DO)WRITE (44.').ts.xs RETURN END REAL FUNCTION AVG‘B(Y.N) c c COMPUTATION OF BULK (MASS-) AVERAGES BASED ON PARABOLIC FLOW c PROFILE c IMPLICIT REAL'B (A-H.K.O-2) DIMENSION Y(N) AVG=Y(1) IF(N.EQ.1)RETURN A1=N AVG=O.DO A=O.DO DO 10 I=1.N B=A A=DMIN1(1.DO/A1¢I.1.DOO) 1O AVG=AVG+Y(I)‘(A‘A‘(2.DO-A°A)-B‘B‘(2.00-8‘8)) RETURN END REAL FUNCTION VS'B(Y) C C COMPUTATION OF SPEED OF SOUND BASED ON IDEAL GAS C IMPLICIT REALFB (A-H.K.O-Z) DIMENSION Y(2) COMMON/BRAVO/AMO.FO.R VS=DSQRT((1.DO+Y(1))tRFY(2)/AMO/(1.DO-Rt(1.DO+Y(1))/ +CPM(Y(1).V(2)))) RETURN END REAL FUNCTION EK‘B(Y.Z) IMPLICIT REAL‘B (A-H.K.O-Z) DIMENSION Y(A) COMMON/BRAVO/AMO.FO.R C C KINETIC ENERGY OF FLOW As FUNCTION OF X. T. P. AND 2 C EK=1.DO/(1.DO/(AMO‘UB(Y.Z)"2)-1.DO/R/Y(2)/(1.DO+V(1))) RETURN END REAL FUNCTION CP‘B(I.T) C C HEAT CAPACITY COMPUTATION FOR PURE GASES BASED ON QUADRATIC 205 C FIT DATA 000 0000000 100 101 102 103 C IMPLICIT REAL‘B (A-H.K.O-Z) COMMON/CHARLIE/DHO.DGO.TO.C(3.3) CP-C(1.I)*T‘(C(2.I)+T‘C(3.I)) RETURN END REAL FUNCTION CPM‘B(X.T) HEAT CAPACITY BASED ON IDEAL MIXTURES IMPLICIT REAL‘B (A-H.K.O-Z) CPMI(1.DO-X)‘CP(1.T)*2.DO‘X‘CP(2.T) RETURN END REAL FUNCTION AR'8(Z) CROSS'SECTIONAL AREA OF TUBE AS FUNCTION OF LENGTH. STRAIGHT TUBE OF AREA AK (USER INPUT. DEFAULT 3.8E-4 50. M.) FOR .025 M. THEN STRAIGHT CONICAL CONVERGENCE TO THROAT OF AREA A2 (USER INPUT. SAME DEFAULT). AT .05 M. THEN STRAIGHT CONICAL DIVERGING SECTION OF .25 M LENGTH TO AK. IMPLICIT REAL‘B (A-H.K.O-Z) COMMON/ALFA/A2.AK DATA 21.22.23.24/0.02SDO..OSDD..OSDO..3000/ AR=AK IF(AK.EQ.A2)RETURN IF((Z.LE.ZT).OR.(Z.GE.24))RETURN IF(Z.LT.22)GO T0 101 IF(Z.GT.23)GO T0 102 AR=A2 RETURN ZA=Z1 ZB=22 GO TO 103 ZA=24 ZB=Z3 AR=(DSQRT(AK)+(DSQRT(AK)-DSQRT(A2))‘(Z-ZA)/(ZA-ZB))“2 RETURN END REAL FUNCTION DAR¢B(Z) C COMPUTES DERIVATIVE OF AREA WITH LENGTH BY FINITE DIFFERENCE METHOD C STEP SIZE 1.E-6 M C C IMPLICIT REAL‘B (A-H.K.O*Z) DY=1.D'06 DAR=(AR(Z+DY)-AR(Z))/DY RETURN END REAL FUNCTION UB‘8(Y.Z) C COMPUTES BULK VELOCITY BASED ON AVERAGE CONDITIONS AND IDEAL GAS LAW C IMPLICIT REAL'B (A-H.K.O-Z) COMMON/BRAVO/AMO.FO.R DIMENSION Y(4) UB=F0‘(1.DO+Y(1))‘R‘Y(2)lV(3)/AR(Z) RETURN END 206 Program 5 Runge-Kutta-4 Integrator Stiff Equations Control, Cascade Error Return, Dependent Variable Control 207 SUBROUTINE RK4(FCN.X,XF.V.W.N.TOL.IF2.IT.zPR) c C INTEGRATION OF DIFFERENTIAL EQUATIONS BY RUNGE-KUTTA 4 METHOD C IMPLICIT REAL‘B (A-H.O-Z) DIMENSION Y(N).W(N.9) DXO=XF-x DXL-1.Do N1=(N-7)/4 N2=N1+6 N3=2'N1+5 CALL FCN(N.X,Y.W(1.B)) IF(IF2.NE.O)RETURN KOUNT=—1 100 KOUNT=KOUNT+1 XC=.IDO XL=O.DO DO 5 I=1.N W(I.2)=W(I.8) 5 W(I.1)=Y(I) GO TO 130 101 CONTINUE CALL FCN(N.X+XL'DXO.W(1.1),W(1.2).DXL‘DXO) IF(IF2.NE.O)RETURN 130 CONTINUE DO 120 K=N2,N3 105 IF(DABS(W(K.2)‘DXL‘DXO).LE.0.100‘DABS(W(K.1)))GO To 114 DXL=DXL/2.00 GO To 105 114 CONTINUE 120 CONTINUE IF(DXL.GT.2.DO“(-20))GO T0 110 WRITE (44.t).' INTEGRATION DIVERGED' IF2=1 RETURN 110 XL=XL+DXL/2.DO X1=X+XL‘DXO DO 10 I=1.N 10 W(I.6)=W(I.1)+DXL‘DXO‘W(I.2)/2.00 CALL FCN(N.X1.W(1.6).W(1.3)) IF(IF2.NE.O)RETURN DO 20 I=1.N 20 W(I.7)=W(I.1)+DXL‘DXO‘W(I.3)/2.DO CALL FCN(N.X1.W(1.7).W(1.4)) IF(IF2.NE.O)RETURN DO 30 I=1.N 30 W(I.7)=W(I.1)+DXL‘DXO‘W(I.4) XL=XL+DXL/2.00 CALL FCN(N.X+XL‘DXO.W(1.7).W(1.5)) IF(IF2.NE.O)RETURN DO 40 I=1.N W(I.1)=W(I.1)+(W(I.2)+2.DO‘(W(I.3)+W(I.4))+W(I.S))‘DXL‘DXO/ +6.00 40 CONTINUE IF (XL-XC) 1001.1000.1000 1000 XC=XC+.100 1001 CONTINUE XIT=0.DO IF(IT.E0.0.)GO T0 901 IF(Y(IT).NE.ZPR)XIT=(W(IT.1)-Y(IT))/(ZPR-Y(IT)) 901 IF(XL.LT.1.00.AND.XIT.LT.1.DO)GO T0 101 DXL=DXL/2.Do IF(KOUNT.E0.0)GO T0 102 ERR=0.DD DO 50 I=1.N IF(W(I.9).NE.Y(I).AND.I.NE.IT)ERR=ERR+ 208 +((W(1.9)-W(I.1))/(W(I.9)-V(I)))"2 50 CONTINUE IF(ERR.LE.TOL)GO T0 103 102 00 60 I=1.N 60 W(I.9)=W(I.1) GO TO 100 103 X=X+XL¢DXO DO 70 I=1.N 70 Y(I)=W(I.9) RETURN END APPENDIX C APPENDIX C: Relations For Physical Properties, Transport Properties, and Rate Constants Physical Properties The reacting gas was assumed to be an ideal mixture of ideal gases. Heat Capacity (CP, CP): For atomic hydrogen - C - 5/2 R - 20786 J / kg-mol K P For molecular hydrogen - CP - 26880 + 4.347 T - 0.0002645 T2 J / kg-mol K (T in K) For vibrationally equilibrated hydrogen. Correlation taken from Balzhiser, et a1. [1972]. For vibrationally excited hydrogen - c T - 7/2 R - 29,100 J / kg-mol K P h V / k T CPV(T) "' (EV(T) / T) e / R EV(T) - NAv h u / (eh V / k T - 1) using h u / k - 5990 K based on the ground-state vibrational frequency for molecular hydrogen we - u / c - 4395.24 cm.1 (Herzberg [1950]). For mixture, based on mole fractions: CPM - [(1 - X) CPH2 + 2 X CPH] / (1 + X) 209 210 Heat of Reaction (DH DH): R) For vibrationally equilibrated system, computed by integrating the ACP from 0 K to the temperature. DHR(O) - 431.41 MJ / kg-mol. Extra- polated from Balzhiser, et a1. [1972] using CP correlation. In vibrationally excited case, use DHR(T,TV) - DHR(O.,0.) + 3/2 R T - E (TV) with DHR(0.,0.) - 431.885 MJ / kg-mol, extrapolated from Balzhiser, g; a1. (Ibid.) value at T - Tv - 298.15 K. Eguilibrium Constant (K, KEQ): Computed using standard relationships; In K - DCO / R T 0 0 d 1n K / d T - DHR / R T2 and above relationships for DHR, based on equilibrium vibrational excitation. DG - 406.48 MJ / kg-mol, at T 0 - 298.15 K, taken from 0 Balzhiser, et al., (Ibid.). 211 Transport Properties Thermal Conductivity: A - KTH - 1.843 x 10'3 TO'8 W / m K (T in K) Fitted to Perry and Chilton [1973] (p. 3-215) corrected for apparent two order of magnitude error in table, discovered by Prandtl number calculations. Assumed independent of atom concentration from Prandtl number argument. Under ideal gas conditions, the Prandtl number, the ratio of the product of the constant-pressure heat capacity and the viscosity to the thermal conductivity is roughly the same for all simple gases, around 0.7, from the kinetic theory, supported by experiment. At low temperatures, this heat capacity for atomic hydrogen is 5/2 R, and that for molecular hydrogen is 7/2 R, on the same bases. When corrected for the molecular weights, on the mass basis appropriate for Prandtl number calculations, the heat capacity at constant pressure for atomic hydrogen is 10/7 that of the molecular gas. Similarly, kinetic theory predicts that the viscosity, to a first, very rough, approximation, should be proportional to the square root of the molecular weight. Browning and Fox [1964] verified this experimentally, measuring a viscosity for atomic hydrogen 0.7 times that of the molecular gas. Hence, the product of the heat capacity and the viscosity is roughly the same for molecular and atomic hydrogen, and, to give the same Prandtl number, the thermal conductivity should also be roughly the same, and independent of composition, in a mixture of atomic and molecular hydrogen. such as a partially dissociated hydrogen gas. This 212 is probably only rigorous at low temperatures, though it was used for all temperatures, and the tabulated values extrapolated freely. -3 T1.72 Diffusion Coefficient: D P - DHZHP - 1.264 x 10 Pa m2/s H-H2 Taken from Lede and Villermaux [1976], extrapolated freely beyond the range reported, 270-320 K. Fits well high temperature dependence of Sancier and Wise [1969]. Assumed constant over all dissociations. Viscosity: p - VISC (For one-dimensional model, friction factor calculation) For mixture, computed from simple mixture heat capacity and conductivity relations (vide supra), assuming Prandtl number of 0.7. 213 Rate Constants l. Recombination Reaction Rates a. Gas Phase Reaction 1. Hydrogen Molecule as third body k _ 2.014 x 1011 T-0.7497 e-20.43/T m6 / kg-mol2 / s (T in K) Fit to curve in Shui and Appleton [1971] ii. Hydrogen Atom as third body k _ 8.816 x 1011 T-o.7219 9-92.09/T m6 / kg-mol2 / s (T in K) Fit to curve in Shui [1973] b. Surface Reaction k - 3.855 x 1018 T'3 e-2900/T m4 / kg-mol / s (T in K) Fit to data for quartz in Wood and Wise [1962]. Assumes second-order surface reaction, and adjusts low-temperature values to give pseudo-second order rate constant. Neglects high-temperature 0. rate constant. Extrapolated freely, assuming temperature reported is surface temperature. 214 2. Collisional Relaxation Rates a. Hydrogen molecule as collision partner .100/171/3 kVT - 25300 e P s'1 (T in K, P in Pa) From Kiefer and Lutz [1966] extrapolated freely to lower temperatures, after Heidner and Kaspar [1972]. b. Hydrogen atom collision partner / exchange relaxation kVT - 6.6 x 10 From Cohen [1978]. 9 e’3800/T P / T s'1 (T in K, P in Pa) BIBLIOGRAPHY 10. 11. l2. 13. REFERENCES Anderson, P., The Makeshift Rocket. Ace Books, N. Y. [1962]. Asmussen, J., Whitehair, S., Nakanishi, 8., "Experiments with a Microwave Electrothermal Thruster Concept", AIAA\JSASS\DGLR 17th International Electric Propulsion Conference, IEPC 84-74 [1984]. Balzhiser, R. E., Samuels, M. R., Eliassen. J. D., Chemical Engineering Thermodypamics, Prentice-Hall, Englewood Cliffs, N. J. [1972]. Bennett, J. E., Blackmore, D. R., "The Measurement of the Rate of Recombination of Hydrogen Atoms at Room Temperature by Means of E. S. R. Spectroscopy", Proc. ROY. Soc. (London), Ser. A, 305, 553 [1968]. Bhatia, Q., S., Hlavacec, V., "Integration of Difficult Initial and Boundary Value Problems by an Arc-Length Strategy", Chem. Eng. Comm.. 22. 287 [1983] Bird, R. B.. Stewart, W. E., Lightfoot, E. N., Transport Phenomena, Wiley & Sons, N. Y. [1960]. Browning, R., Fox, J. W., "Coefficient of Viscosity of Atomic Hydrogen and the Coefficient of Mutual Diffusion for Atomic and Molecular Hydrogen", ProcI Roy, SocI (London), Ser. A., 278(1373), 274 [1964]. Cash, J. R., "A Class of Implicit Runge-Kutta Methods for the Numerical Integration of Stlff Ordinary Differential Equations", J; Assoc, Computing Machinery, 22 (4), 504-511 [1975]. Chapman, R., Filpus, J., Morin, T., Snellenberger, R., Asmussen, J., Hawley, M., Kerber, R., "Microwave-Plasma Generation of Hydrogen Atoms for Rocket Propulsion", J. Spaceprgfp and Rockets, 19, 579 [1982]. Cohen, N., A Review of Rate Coefficients in the H2;§2 Chemical Laser System - fippplement (1977), TR-0078(3603)-4, The Aerospace Corp., El Segundo, California [1978]. Coste, J., Rudd, D., Amundson, N. R., "Taylor Diffusion in Tubular Reactors", Can. J. Chem. Eng., 39, 149 [1961]. Filpus. J., Hawley, M., "The Energetics of Hydrogen Atom Recombination: Theory, Experiments, and Modelling", AIAA/JSASS/DGLR 17th Internatonal Electric Propulsion Conference, IEPC 84-77 [1984]. Golden, S., "Free-Radical Stabilization in Condensed Phases", Jy_ Chem, Phys., 29, 61 [1958]. 215 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 216 Greco, R. V., Byke, R. M., "Resistojet Biowaste Utilization-- Evaluation and System Selection", J, Spacecraft 5nd Rockets, 6(1), 37 [1969]. Halbach, C. R., Yoshida, R. Y., "Development of a Biowaste Resistojet". J. Spacecraft and Rockets, 8(3), 273 [1971]. Hawkins, C., Chapman, R., private communication [1982]. Hawkins, C., Nakanishi, 8., Free Radical Propulsion Concept, NASA Technical Memorandum 81770 [1981]. Heidner, R. F., III, Kaspar, J. V. V., "An Experimental Rate Constant for H + H2(u" - l) » H + H2(v" - 0)", Chem. Phy§.,L§tt., 15(2), 179 [1972]. Herzberg, G., Spectr§_of Diatomic Molecules, 2nd ed., Von Nostrand, N. Y., [1950]. Kallis, J. M., Goodman, M., Halbach, C. R., "Viscous Effects on Biowaste Resistojet Nozzle Performance", J, Spacecrafp and Rockets, 9(12), 869 [1972]. Kiefer, J. H., Lutz, R. W., "Vibrational Relaxation of Hydrogen", Jp_ Chem, Phys,, 44, 668 [1966]. Kingsbury, D., "Atomic Rockets", Analog, 95(12). 38 (Dec.)[1975]. Lede, J., Villermaux, J., "Measurement Of Atom Diffusivity by a Method Independent of the Wall Effects", Chem. Phys. Letp,, 43(2), 283 [1976]. Levitskii, A. A., Polak, L. 8., "Study of the Recombination Reaction H + H + H » H2 + H by the Classical Trajectories Method", Khim. Vys. Energii, 12, 291 [1973]. Mallavarpu, R., Asmussen, J., Hawley, M. C., "Behavior of a Microwave Cavity Discharge over a Wide Range of Pressures and Flow Rates", IEEE Ipans. Plasma Science, PS-6(4), 341 [1978]. Mearns, A. M., Ekinci, E., "Hydrogen Dissociation in a Microwave Discharge", J, Micpowave Power, 12, 156 [1977]. Mertz, S. F., Asmussen, J., Hawley, M. C., An Experimental Study of C0 and Hzin a Continuous Flow Microwave Discharge Reactor". IEEE Ipans, Plasma §gience, PS-2(4), 297 [1974]. Mertz, S. F., Asmussen, J., Hawley, M. C., "A Kinetic Model for the Reactions of C0 and H2 to CHhand 02H2 in a Flow", IEEE Trans. Plasma Science, PS-4(1), 11 [1976]. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 217 Morin, T., Chapman, R., Filpus, J., Hawley, M., Kerber, R., Asmussen, J., Nakanishi, 8., "Measurements of Energy Distribution and Thrust for Microwave Plasma Coupling of Electrical Energy to Hydrogen for Propulsion", AIAA/JSASS/DGLR 16th International Electric Propulsion Conference [1982]. O'Neill, G. K., The High Frontier, [1977]. Palmer, H. B.. "Burning Rate of an H-Atom Propellant", ARS J., 29, 365 [1959]. Perry, R. H., Chilton, C. H., eds., Chemical Engineers' Handbook, 5th ed., McGraw-Hill, N. Y. [1973]. Sancier, K. M., Wise, H., "Diffusion and Heterogeneous Reaction. XI. Diffusion Coefficient Measurements for a Gas Mixture of Atomic and Molecular Hydrogen". J, Chem, Phys,, 51, 1434 [1969]. Shapiro, A. H., e D namics nd ermod amics of Com ressible Fluid Flow, Vol. 1, Ronald Press. N. Y., 237 [1953]. Shaw, T. M., "Dissociation of Hydrogen in a Microwave Discharge", ;y_ Shem. Phys., 30, 1366 [1959]. Shui, V. H., "Thermal Dissociation and Recombination of Hydrogen According to the Reactions H2 + H «a H + H + H", 4. Chem. Phys,, 58, 4868 [1973]. Shui, V. H., Appleton, J. P., "Gas-Phase Recombination of Hydrogen. A Comparison between Theory and Experiment", J. Chem Phys., 55, 3126 [1971]. Silvera. I. F., Walraven, J., ”The Stabilization of Atomic Hydrogen", Scientific American, 246(1), 66 (Jan.) [1982]. Snellenberger, R. W., "Hydrogen Atom Generation and Energy Balance - Literature Review, Modeling, and Experimental Apparatus Design", M. S. Thesis, Michigan State University [1980]. Stewart, P. A. E., "Liquid Hydrogen as a Working Fluid in Advanced Propulsion Systems", JI B. I. 8., 18, 225 [1962]. Sutton, G. P., Ross, D. M., Rocket Propuision Elements, 4th ed., Wiley, N. Y. [1976]. Trainor, D. W., Ham, D. 0., Kaufman, F., "Gas Phase Recombination of Hydrogen and Deuterium Atoms", J, Chem. Phys., 58, 4599 [1973]. Villermaux, J., "Etude de L'Hydrogene Actif. I. Analyse Thermocinetique de la Recombinaison de L'Hydrogene Atomique", ;y_ Chim. £hys., 61, 1023 [1964a]. 44. 45. 46. 47. 48. 218 Villermaux, J., "Etude de L'Hydrogene Actif. II. Etude de la Relaxation de Recombinaison Homogene de L'Hydrogene Atomique", J. Chim. Phys., 61, 1033 [1964b]. Whitehair, S.. Asmussen, J., Nakanishi, S., "Demonstration of a New Electrothermal Thruster Concept", Appi. Epya. Letp., 44(10), 1014- 1016 [1984] Whitehair, S., Asmussen, J., Nakanishi, 8.. "Microwave Electrothermal Thruster Performance in Helium Gas", accepted by l; Propulsion and Power [1986]. Whitlock, P. A., Muckerman, J. T., Roberts, R. E., "Classical Mechanics of Recombination via the Resonance Complex Mechanism: H + H + M 4 H + M for M - H H He, and Ar", J. Cham. Phys., 60, 3658 2 i 2 I [1974]. Wood, B. J., Wise, H., "The Kinetics of Hydrogen Atom Recombination on Pyrex Glass and Fused Quartz", J. Shem. £bys., 66, 1049 [1962]. GENERAL BIBLIOGRAPHY Amdur, I.. Robinson, A. L.. "The Recombination of Hydrogen Atoms. I ", J. A. C. 8., 55, 1395 [1933]. Amdur, I.. ”The Recombination of Hydrogen Atoms. 11. Relative Recombination Rates of Atomic Hydrogen and Atomic Deuterium", J. A. C. S., 57. 856 [1935]. Amdur, I.. "The Recombination of Hydrogen Atoms. 111.", J. A. C. 5., 60, 2347 [1938]. Audibert, M. M., Joffrin, C., Ducuing, J., "Vibrational Relaxation of H2 in the Range 500-40°K", Chem. Phys. Lett., 25(2). 158 [1974]. Audibert, M. M. , Vilaseca, R. Lukasik, J. , Ducuing, J., "Vibrational Relaxation of Ortho and Para- -H2 in the Range 500- 40° K", Chem. Phys. Lett., 31(2) 232 [1975]. Bancroft, J. L.. Holeas, J. M., Ramsay, D. A., "The Absorption Spectra of HNO and DNO". San. J. Phys., 40, 322 [1962]. Brown, S. C., Basic Data of Plasma Physics, Wiley & Sons, N. Y. [1959]. Clyne, M. A. A., Stedman, D. H., "Reactions of atomic H with HCl 2 and N001", Trans. Faraday Soc., 62, 2164 [1966]. Poole, H. G., "Atomic Hydrogen. I. Calorimetry of Hydrogen Atoms", Epoc. Roy. Soc. (Lppdgp), Ser. A., 162, 404 [19378]. 10. 11. 12. 219 Poole, H. 6., ”Atomic Hydrogen. II. Surface Effects in the Discharge Tube", Eppp. Roy. Soc. (London), Ser. A., 162, 415 [1937b]. Poole. H. G., "Atomic Hydrogen. III. The Energy Efficiency of Atom Production in a Glow Discharge". oc Ro oc ndon , Ser. A., 162, 424 [1937c]. Sutherland, G. 8., "Recent Advances in Space Propulsion", ARS J., 29, 698 [1959]. HICH G 1111111[111111115