UM 320057-3-P COMBUSTION DYNAMICS AS RELATED TO AIR POLLUTION by J. A. Nicholls M. Sichel J. Ay D. E. Geister H. C. Liang Gas Dynamics Laboratories Department of Aerospace Engineering The University of Michigan prepared for American Gas Association A. G.A. Project BR-26-1 L.A. Sarkes, Project Monitor September 1973

TABLE OF CONTENTS Page FOREWORD SUMMARY LIST OF FIGURES iv LIST OF TABLES v I. INTRODUCTION 1 II. ANALYTICAL STUDIES 3 A. Mechanism of NO Generation in Exhaust Systems 3 B. Equilibrium Theories of NO Formation 6 C. NO Formation in a Streamtube With Variable Temperature 13 D. The Detailed Kinetic Model 16 E. Discussion of Analytical Results 22 III. EXPERIMENTAL ARRANGEMENT AND RESULTS 124 A. Facility 24 B. Results 27 REFERENCES 29 APPENDIX A 32

FOREWORD This report describes the progress made in the past year on American Gas Association Project BR-26-1. This study is being conducted under the direction of Prof. J.A. Nicholls of the University of Michigan while Mr. Louis A. Sarkes of A. G.A. is the project monitor. Professor M. Sichel is leading the analytical portions while Mr. D.E. Geister, Research Engineer, assists on the experimental aspects. Mr. Ay is conducting the analytical work as part of his Ph.D. thesis. Mr. Liang, graduate student, is doing the experimental work.

SUMMARY The objective of this study is to increase the fundamental knowledge of the combustion of natural gas and to reduce the pollutants in the combustion products that are exhausted to the atmosphere. There are two parts to the problem; first the generation of pollutants in the main combustion zone and secondly, the change in the pollutant concentrations as the combustion products flow through the exhaust system. In assessing these aspects various chemical kinetic schemes for the formation of NO have been considered and characteristic formation times predicted. Flows in the exhaust system where the temperature drops due to heat loss have been evaluated and "freeze-out" concentrations of NO estimated. Further, the numerical integration of the stiff differential equations describing the constant pressure combustion of H2-02 has been successfully carried out. The method will now be applied to the methane-air case. An experimental facility simulating a simple exhaust system has been set up and operated. Experiments have been run with methane-air. Satisfactory temperature profiles and qualitative species concentrations determined by a mass spectrograph have been obtained. iii

LIST OF FIGURES Figure 1. Flow in the Exhaust System Figure 2. X versus t Profiles for Various Simplified Equilibrium l'oels of CH4-Air Mixture at $ = 1.0, p = 1 atm and Tc = 22270K. Figure 3. XNO versus t. Profiles for CH4-Air Mixture of 4 = 0. 6 and 1.6 (2-Reaction and 3-Reaction Models) p = 1 atm. Figure 4. Characteristic Time, TNO of NO formation versus EquivalenceRatio, 4, Profiles for CH4-Air Mixture at 1 atm and To = 298 and 10000K. Figure 5. Characteristic Time, TNO of NO Formation versus Equivalence Ratio,,, for CH4-02-N2 Mixture at p = 1 atm, T Reactants = 298. 150K and 10000K. Figure 6. XNO versus IdT/dtl Curves at Various x/vc of CH4-Air Mixture at= 1. O, p = atm. Figure 7. T versus x Curve (T = (T/TC, x = (x/L), L = (Tcvc/a), and dT/dt = - a = constant). Figure 8. Mole Fraction and Temperature Profiles of H2-02 Mixture at 4 = 1.0 and 1 atm. Figure 9. Combustion Chamber and Exhaust Stack. Figure 10. Water Cooled-Gas Probe. Figure 11. CH4-Air Axial Temperature Profile, 4 = 1.0. Figure 12. CH4-Air Axial Temperature Profile; Stoichiometric, Lean, Rich. Figure 13. CH4-Air, Temperature vs Equivalence Ratio. iv

LIST OF TABLES Table 1. Reaction Constants of Nitric Oxide Mechanism Table 2. Chemical Model Numerical System of H2-02 Mixture Table 3. Forward Rate Constants of H2-02 System Table 4. Equilibrium Constants of H2-02 System Table 5. Coefficients for Evaluation of Molar Enthalpies of Species

I. INTRODUCTION Natural gas is one of the cleaner fuels when it comes to the production of pollutants. However, the vast quantities that are burned and exhausted to the atmosphere in our technological society make it important to minimize the resulting air pollution. The objective of this study, then, is to increase the fundamental knowledge of the combustion of natural gas and to reduce the pollutants in the combustion products that are exhausted to the atmosphere. There are two main factors that enter; first there are the pollutants generated in the primary combustion process and, second, there is the subsequent generation and/ or disappearance of pollutants in the exhaust system. In order to assess the primary combustion zone, various chemical kinetic schemes are being for the production of NO in methane-air are being analyzed in order to predict the NO levels possible, the times required to produce these levels, and hence possible means of control preventing the higher concentrations. The exhaust phase is also being treated analytically with the aim of identifying the importance of heat loss to the walls, residence time (flow velocity), duct size and material, and nature of the flow. An experimental facility has been developed -to study this latter phase and to integrate the results with the theory. 1

The analytical work will be described in detail in the next section followed by a description of the experimental work. 2

II. ANALYTICAL STUDIES A. MECHANISM OF NO GENERATION IN EXHAUST SYSTEMS As indicated in the Introduction, the main purpose of the present investigation is to determine the influence of exhaust system parameters upon pollutant generation. Development of a reasonably accurate analytical model which will predict the influence of exhaust system parameters upon pollutant generation is the central objective of the analysis. This model, which must be substantiated by experiment, should supplement and guide the experimental portion of the present study; suggest and assess exhaust system modifications that will reduce pollutant emission; and permit extrapolation of the experimental results to exhaust systems other than the specific one being studied here. A simple exhaust system corresponding to that being studied experimentally is shown in Fig. 1. The fuel-air mixture enters the burner with velocity vo and temperature To. vc and T are the velocity and temperature downstream of the flame with T the adiabatic flame temperature if combustion is complete. Fluid particles in different streamtubes are subject to different residence times and temperature histories even when the exhaust system is a simple tube because of the turbulent boundary layer and the temperature variations across the duct. As indicated by Nicholls et al., the first step in dealing with complex systems such as that shown in Fig. 1 is to analyze the processes which 3

occur in a single streamtube. The streamtube analysis of NO generation in the products of methane-air combustion has been the main subject of research during the current contract year. Experimental evidence indicates that the key nitric oxide formation reactions are the following two, which were first proposed by Zeldovitch et al.3 kif O + N2 NO +N ( klb k2f N + o2 I NO + O (In k2b while a third reaction k3f N + OH- NO + H (II) k3b becomes important for rich mixtures, and low temperatures. The forward and reverse rate constants above used by Bowman2 and several others are tabulated by Nicholls et al. Reactions (I), (It, and (TI) are slower than the main combustion reactions so that the bulk of NO formation occurs downstream of the flame front. Zeldovitch et al. were the first to suggest that the NO reactions are thus almost completely uncoupled from the kinetics of the hydrocarbon combustion. This idea 4

has been the basis for a number of simplified models of NO chemistry (e.g. Refs. 4-8) and for various techniques for reducing NO generation by manipulation of the post flame region. In systems such as the internal combustion engine the post flame region, which lies within the cylinder is not readily accessible. However, in commercial burners 9 10 using natural gas such as discussed by Kremer, and Sarofim et al. among others, it is possible to reduce NO generation by appropriate modification of the flow downstream of the burner. It is of course this fact which provides the impetus for the present study. Kinetic schemes for the formation of NO, based mainly upon reactions (1), (If), and (III), have been proposed by a number of investigators in recent years. Approximations have been introduced to obtain solutions for the concentration of NO in simple analytical form. Some of the more common kinetic schemes have been compared in Section B below by using them to compute the post flame formation of NO in an isothermal streamtube. In a real exhaust system the flow downstream of the flame will not be isothermal, rather, the temperature will decrease at a rate which depends on the details of the flow through the exhaust system. In order to investigate the effect of post flame cooling the Zeldovitch mechanism has been used to compute the rate of NO formation in a streamtube with decreasing temperature as described in Section C below.

A common feature of the post flame theories of NO generation is that the concentrations of some of the species which take part in reactions (1), (II), and (II1), are constant at the post flame equilibrium value. Bowman 1 therefore refers to these as "equilibrium theories". The concentration and rate of formation of NO obtained from these equilibrium theories are always below experimentally observed values (Ref. 2, 10, 11). Bowman and Sarofim and Pohl12 suggest that these discrepancies arise because the concentrations of 0, 02, and OH may be orders of magnitude above the post flame equilibrium values. Once this factor is taken into account, reactions (), (IU), and (IIi), appear to yield reasonable results. An analysis of this effect requires detailed computation of flame kinetics. The numerical solution of the kinetic equations is one of the key difficulties in such an analysis. This step has been completed using a relatively new computational scheme, STIFF3, as described in Section D below. B. EQUILIBRIUM THEORIES OF NO FORMATION As indicated above, reactions I, II, and III, are the basis of a number of approximate kinetic schemes for the post flame generation of NO. Some of the more common schemes are presented below and have been used to compute the variation of NO downstream of the flame, where the flow is assumed to be adiabatic and at atmospheric pressure. Since the mass fraction of NO is always very small, the NO reactions 6

have a negligible effect on the enthalpy of the fluid so that the flow will be isothermal. The rate constants which have been used for reactions I, II, and III, are tabulated in Table 1. Each of the reaction schemes considered is described in detail below: 1. Newhall Model Reactions (1) and (II). Concentration of O, N2, N, 02 at the post flame equilibrium value. Only NO concentration is variable. The analytic solution for XNO, the ratio of mole fraction, XNO of NO TheNXNO O o to that of equilibrium value, XNOe, at adiabatic flame temperature, Tc, is then XNO = - exp RT (klb XN +k2b O t (1) C T=T where the subscript e refers to the post flame equilibrium value. 2. Marteney Model4 Reactions (1) and (II) plus the reaction k4f NO+N2 kf N+O+N2 (IV) 4b Only forward reactions considered. Concentrations other than NO and N at the post flame equilibrium value. Analytic solution for XNO is then 7

2P klf Xe XN e RT XNO T T=Tc klf Oe X eN2 /kX f1i2f f (1 - exp 2f NO) t (2) RT T=T for k X2f Xe >> k4f XN and negligibly small initial value of XNo and XN. 3. Zeldovich Model3 Reactions (1) and (II) Steady state approximation for 0 and N, i. e. dXN dXo dt - dt Concentration of 0, 02 and N2 at post flame equilibrium value. An implicit analytical solution for 5NO is then

2P k2b XO e R tan Tb X /(XNo T=T f N T-T C -ln(1 xNOe ) (I 4. Zeldovich Model3 With k2f 0 - >> kb NO T=T C The analytic solution becomes,V2Pk k X e=. e\ Xon tanh I2P k1 f k2b NOe Xt (4) NNO Tkf e RT kaf X T=T C 5. Zeldovich Two-Reaction Model Steady state approximation relaxed Concentration other than NO and N at the post flame equilibrium value Numerical integration of equations using STIFF3 algorithm2

6. Three-Reaction Model Reactions (1), (It), and (II1) Concentrations other than NO and N at the post flame equilibrium value Numerical integration STIFF3. The last two models, the Zeldovich two-reaction model and the threereaction model, are the simplified kinetic models employed in the studies presented below. The former model has the fewest approximations in comparison with all other two-reaction models. Model No. 6 introduces reaction (II) with, however, the same assumptions as model No. 5. The variation of XNO, the mole fraction of NO, with time for all the models mentioned above is shown in Fig. 2 for the region downstream of a stoichiometric methane-air flame (equivalence ratio q = 1.0) at atmospheric pressure. From Fig. 2 it can be seen that the Zeldovich, Newhall, and Zeldovich two-reaction models (1, 3, and 5) yield similar results. A detailed study of the two-reaction models shows that the discrepancy in initial rates of NO formation is basically due to the different assumptions regarding the concentration of the species involved. The N concentration of the Zeldovich model with the steady state approximation (model No. 3), for N, i.e., d XN/dt = 0, decreases from its initial value asymptotically to the post flame equilibrium value. The concentration of N in the Zeldovich two-reaction model without any restriction 10

on N rises rapidly from zero to the steady state value and then relaxes slowly to the post flame equilibrium value. Both of the above two models give a higher initial rate of NO formation than the Newhall model, which assumes a constant post flame equilibrium concentration for N. The Marteney model (2) can give reasonable results only for very short times since reverse reactions are neglected. The three reaction model (6) and the approximate Zeldovich model (4) give similar results but this similarity is probably coincidental. Further studies verify Bowman's2' remark that reaction (IITI plays an important role in NO formation only behind rich flames. This can be seen from Fig. 3 which compares the Zeldovich two-reaction (5), and three-reaction (6) models for equivalence ratios of 0. 6 and 1. 6. For the lean flame, J = 0. 6, the results of the two models are indistinguishable but for the very rich 1 = 1.6 flame the rates of NO formation differs by two orders of magnitude. For these two extreme cases the rate of NO formation is in any case very small. In the range 0. 8 < p < 1.0, reaction (III) while important will not have a major effect on the rate of NO production. A characteristic NO formation time, r, can be readily obtained from the rate equation of the Zeldovich model (3). From the rate equation in the form 11

d XNO 1 -XNO dt T(5)'NO it is possible to identify the characteristic reaction time k2f (Xo + klb (XNNO 2 02 T=T ccT=T T=-T c c The time TNO provides a measure of the time required for NO formation behind the flame. The calculated -NO versus equivalence ratio, 4, curves for an atmospheric methane-air flame at two different initial reactant temperatures are shown in Fig. 4. As is to be expected TNO reaches a minimum value at k = 1.0 where the adiabatic flame temperature, Tc, reaches its maximum. The minimum value of'NO at k = 1.0 can be reduced by an order of magnitude if the reactants are preheated from 298~K to 1000 K as shown in Fig. 4. The same TNO curves for CH4-N2-o2 mixtures with the ratio of N2 to 02 by volume as a parameter are shown in Fig. 5. For Oxygenrich flames with N2/02 = 1/4, and 1.0, the minimum values of TNO are three orders of magnitude lower than for air. The minima of the TNO curves are very flat for the oxygen-rich and preheated flames as compared to the methane air flame without preheating. 12

C. NO FORMATION IN A STREAMTUBE WITH VARIABLE TEMPERATURE In order to investigate the effect of post flame cooling a streamtube with a temperature variation such that the Lagrangian time derivative of the temperature is constant was considered. Thus, the temperature of each fluid particle decreases at a constant rate so that dT dT dt =v dx = a = const (7) As shown by Nicholls et al.l the molecular weight and pressure downstream of the flame will be approximately constant. Then the continuity equation becomes Pc c= pv (8) so that PC T v -v = v (9) p c T C and from (1) and (3) the spatial distribution of temperature, T(x), is then 1/2 T: (1_-2xL) (10) T;\ L c where L is a characteristic cooling length given by v T c C L= c c (11) a 13

Here v and p are the velocity and density downstream of the flame, and L is of the order of duct length required for the gas to cool from a temperature of Tc to zero with dT/dt = - a. The temperature distribution corresponding to (4) is shown in Fig. 7. The Zeldovich Model (No. 5) was used to compute the rate of NO formation; however, with the rate constants evaluated at the local value of the temperature T. The resultant rate equations are then dXNO P e e dt RT (kif XO XN2 2f N 02XO lb N NO -k 2bX eXN (12) dXNP k XX k X X dt =RT (klfXO XN2 2f XN 02 lbNOXN +k2bXO e XNO) (13) 2b OXN The results of these calculations for f = 1.0 are given in Fig. 6 which shows the variation of XNO the ratio of NO concentration to the NO equilibrium concentration (XNOe) T=T at the adiabatic flame temperac ture, with the cooling rate, - a = dT/dt, and a residence time x/vc. For each value of the cooling rate XNO approaches a limiting "freeze out" value with increasing distance x along the stream tube. With increasing IdT/dtl 1 the maximum or "freeze out" value of XNO decreases. Such 14

chemical freezing also occurs in high temperature nozzles as discussed by Bray 8' 9,10 Hall and Russo and Cheng and Lee' among others. The nature of such reactive flows depends on the relative values of the characteristic flow time Tf and a characteristic chemical reaction time T. Detailed study of the numerical results shows that the "freezing" of NO occurs in the temperature range Tc > T > (Tc - 200) OK. As can be seen from Fig. 7, the temperature decreases almost linearly with distance in this range, and this result has been used as the basis of an analytical study of NO freezing. The details of this analysis will be presented in a later report; however, preliminary results indicate that, X*, the freeze out value of XNO and the cooling rate dT/dt satisfy a relation of the form (X*)n dT= f( (14) dt = f(p) where f(q) is a function which depends upon the fuel under consideration and the equivalence ratio, $, of the combustible mixture. The exponent n lies in the range 1 < n < 2. For a stoichiometric CH -air mixture at atmospheric pressure, n = 1. 6 and f(A) = 65. The value of X* calculated from (14) is compared to the numerically determined value in Fig. 6 and the two values are found to be in excellent a greement in the range 2 x 102 < I dT/dt I < 2 x 103 (~K/sec). Relations of the form (14) should be extremely 15

useful in dealing with complex exhaust systems in which the temperaturetime variation is different for each streamtube. D. THE DETAILED KINETIC MODEL The concentrations and rates of formation of NO calculated using the simplified streamtube models discussed in Section B, are found to be below those observed experimentally. Another related obser vation is that the experimental variation of NO when extrapolated backward leads to a non-zero value of NO concentration at the flame, sometimes referred to' 19 as "quick NO". Bowman ascribed this effect to the non-equilibrium concentrations of the species, 0, 02 and OH, immediately behind the flame front. However, once the proper nonequilibrium values are used as initial conditions, the Zeldovich mechanism appears to properly predict the formation of NO 2 1 As indicated in Section A, determination of the nonequilibrium values of 0, O22 and OH behind the flame requires detailed consideration of flame kinetics. The solution of the kinetic equations is a key source of difficulty in such an analysis. As a first step in the kinetic calculation a relatively new computational scheme called STIFF3 has been successfully applied to the calculation of H2-02 kinetics. The detailed evolution of the H2-02 reaction at atmospheric pressure has been determined and the results for a stoichiometric mixture with an initial temperature of 1000~K are shown in Fig. 9. The 16

concentration ofthe species are shown as a function of time; however, the time scale is readily converted to a distance scale if the mass flow density of the flow is specified. The influence of heat conduction and diffusion, which plays an important role in flame structure has been omitted in the present calculations, which are, essentially, a test of the STIFF3 algorithm. Eight chemical reactions involving the six species H2, 02, H, O, 2. OH, and H20, as listed in Table 2, were used in these calculations. The forward rate constants of the reactions are listed in Table 3. The reverse rate constants were evaluated from the equilibrium constants of the reactions taken from the JANAF Thermodynamic Tables, and listed in Table 4. The enthalpy of each species is calculated by the polynomial equation H. 2 3 4 A6'-A +A A1 T T T _(15) RT 1 2 + A3 3 +A4 + A5- + - (15) These coefficients are the same as those evaluated and tabulated in the NASA Program (SP-273) by McBride and Gordon. The enthalpies calculated from the above equation have the same values as those 20 listed in the JANAF Thermodynamic Tables. The coefficients in (15) are tabulated in Table 5. The six species equations written in terms of the logarithm of (Yi/MI) are listed below. Use of the logarithmic form permitted 17

larger step sizes in the numerical integration. d H X dt PZ (16) dt p Z dZo co H2o d Z00 WO (18) dt p Z0O dZOH OH0 2= 2 (19) dt p ZO d Z OH0OH H~, e (20) dt p Zo 2 dZH wH dt pZ H2 where Z. = In Zi = In (Yi/Mi) (22) Y. = mass fraction of species i M. = molecular weight of species i o. = the net chemical rate of formation of species i, 1 and 18

wH = - Y1 + Y2 - Y3 - 2Y5 - Y7 (23) ~ = Y1 + Y2 - Y3 - 2Y5 - Y7 (24) - Y3 - Y4 + Y6 (25) WOH = Y1 + Y2 + Y3 + 2Y4 -Y6 (26) wO = - Y1 + Y8 (27) WH = - Y2 + Y3 + Y4 (28) 2 and Y1= p2 (klf ZH Zo klb ZOH Z) (29) 2 Y2 = p (k2f ZO Z H - k2b ZOH ZH) (30) Y3 = p (k3f ZHZHO -k3b Z ZH)(31) 3f H H~0 2 (31) 2 2 Y4 = p (k4f ZO ZH20 k4b ZOH (32) 4S —-H 2f*H kb O 3=p Z H k5b 22 ) (33) Y5= Z k,33 Y (k6f O ZOH k6b ) (34) M7'Z7f ZH k7b pOH (35) 3 (k-k (36) M8 0 8bOp

The following two equations for the conservation of hydrogen and oxygen atoms were used in addition to the species equations. H Conservation Z + 2Z + ZO + 2ZH =constant (37) O Conservation O +Z +Z +2Z2 =constant (38) 0 H20 OH 0 The energy equation was written in differential form as indicated below: (HH' O O H20 H20 OH OH +H o +HH d2 2 2 2 dT ~~~~d~~~~t~~~~~~~~ = ~(39) dt / d HH dH0 dH20 Z dT ZO dT H2O dT d dH dH d HOH 02 dH2 oH dT o2 dT where H. = molar enthalpy of species i and is determined using (15). The above set of differential equations must be integrated numerically. However there is a large variation in the time constant associated with the different chemical reactions so that the above system is what is referred to as a stiff system of equations. The usual methods of 20

numerically integrating ordinary differential equations, e. g. Runge Kutta, Adams-Moulton, and Hamming's methods, then require inordinately small step sizes for a stable solution. Special methods have been developed for dealing with such equations including the algorithm22 STIFF3 used in the present study. The method uses an implicit technique to solve the equations which are first locally linearized. With this technique the step size is no longer controlled by the largest negative eigenvalue of the system. A detailed description of the program used including a computer printout is presented in Appendix A. The numerical results for the reaction of a stoichiometric H2-02 mixture at atmospheric pressure and an initial temperature of 10000K is shown in Fig. 9. The reaction can be divided into an induction region, a main reaction region, and a recombination region. In the induction region the mole fractions of O, OH, H, and H20 increase rapidly from zero with a negligible change in the temperature, and the mole fractions of the reactants H2 and 02. In the main reaction region the major portion of the reactants is consumed, the bulk of H2O0 is formed, the temperature rises rapidly, and the maximum overshoot in O and H occurs. In the recombination region the free radicals recombine and the mole fractions and temperature approach the equilibrium values. A check on the non-equilibrium calculations is provided by the fact that the final values of the mole fractions and the temperatures agree 21

with the equilibrium values computed for the same mixture with the McBride-Gordon program. The STIFF3 algorithm is very efficient in that only 2 minutes of central processor time were required to complete the calculation described above. E. DISCUSSION OF ANALYTICAL RESULTS Comparison of the various "Equilibrium Models" for NO formation based on the Zeldovich mechanism, shows that they yield comparable results with the exception of the Marteney model which neglects reverse reactions. Relaxation of the steady state approximation for N results in little change in the final results. In fact, the concentration of N was found to rapidly approach the steady state value. It should be possible to develop analytical criteria for the range of validity of the steady state approximation. The calculations presented here confirm the conclusion of Bowman 2 1 that the reaction III must be added to the Zeldovich reactions I and II for rich mixtures. The characteristic reaction times shown in Figs. 4 and 5 provide an excellent indication of rapidity of NO formation behind CH4-flames under a variety of conditions. If the velocity of the post flame gases is known, this time is readily converted to an NO formation distance. From these figures it is clear that NO forms most rapidly near the stoichiometric mixture ratio. NO formation is particularly rapid when the reactants are preheated and the amount of 02 relative to N2 exceeds that for air. 22

With cooling behind the flame front the NO formation reactions freeze below a certain temperature level so that the final NO mole fraction is less than the equilibrium value at the adiabatic flame temperature. Preliminary analysis indicates that simple relations can be developed between the mole fraction of NO at freeze out and the temperature gradient behind the flame front. Such relations will be very useful in estimating NO formation in complex systems and in evaluating the effectiveness of post flame cooling devices for the reduction of NO. As a first step in the evaluation of non equilibrium concentrations behind the flame front the algorithm STIFF3 has been successfully applied to the calculation of H2-02 kinetics. Succeeding steps in this study will be to treat the more complicated kinetics of CH4 combustion followed by the inclusion of heat and mass transfer within the flame front. It is hoped that these studies will lead to simple criteria for estimating non equilibrium concentrations behind the flame. 23

III. EXPERIMENTAL ARRANGEMENT AND RESULTS A. Facility The experimental setup that has evolved for studying the products of combustion in the exhaust system for a methane-air mixture is shown in Fig. (9). The combustion and exhaust systems are oriented vertically. The combustion chamber is a commercially available wall burner capable of operation over an extremely large range of fuel to air ratios. The burner consists of a gas-air mixing tube, injection nozzle, and a ceramic burner block with a mounting flange. Combustion occurs at the burner which is located at the base of the burner block. A pilot is installed in the burner to assure and maintain combustion at all times. The exhaust system simulates a common commercial chimney and consists of an 8 ft length of 4 in. steel tube. Flow sampling ports have been provided at one foot intervals on the tube wall. The flow control system for the burner and pilot has also been installed on the test fixture. This system consists of a number of metering valves, rotameters, pressure gauges, and turbine flowmeters, and allows individual control and monitoring of both fuel and air flows to the burner and pilot. Each unit is independently controlled and monitored. In this facility methane is supplied by a pressurized cylinder and filtered air from the laboratory's high pressure air supply. These 24

gases enter the stand separately through quick-disconnect fittings. Rotameters provide for measurement of the fuel and air flow rates for both the main burner and the pilot. Inlet pressures of the CH4 and air at the burner assembly are measured by inclined manometers and have been found to be the same as the pressures of the supplied fuel and air. The pilot mixture is first ignited by an electrical spark; it then is used to ignite the main burner mixture. The pilot may burn continuously during operation of the main burner or it may be turned off and re-ignited at anytime. The exhaust duct is mounted at the exit of the burner block. The stand is located so that the combustion gases leaving the exhaust duct are drawn through louvers in the ceiling by an exhaust fan and released to the atmosphere at the roof of the building. Measurements are made on the combustion products as they pass through and along the exhaust system. Probe access holes are located at one foot intervals along the flow direction. Experiments were conducted wherein the chemical composition of the combustion product gases was indicated by use of a mass spectrometer. Temperature is measured using a pt/pt + 10% Rh thermocouple. The preliminary experiments indicated that unduly large heat loss was being experienced. In order to rectify this, two major changes were incorporated. First, the steel pipe was replaced by a 25

3 in. I.D. (4 1/2 in. O.D.) ceramic tube. The low thermal conductivity of this tube lowers the heat losses. Further, the smaller diameter leads to higher flow velocities with consequent lower residence times and heat loss. The second change consisted of modification of the burner in order that better flame holding and a more uniform temperature profile would be obtained. Two ceramic slabs with many axial holes in them were mounted in the burner block as shown in Fig. ( 9 ). One is 3 in. above the exit of the burner and the second is 3 in. higher. Alumina pebbles (1/ 3 in. diameter) were put between the two slabs and also on top of the second one. Finally, a 3 in. thick blanket of ceramic fiber was wrapped around the burner block and the lower part of the ceramic tube to further reduce heat losses. In order to check the exhaust temperature measured by the pt/pt + 10%1 Rh thermocouple and compare it with the theoretically computed exhaust temperature for the gases flowing in the exhaust duct, the Extreme-Temperature Gas Sampling probe (Greyrad Probe) was used to deduce primary flow characteristics such as the total enthalpy and chemical composition along the exhaust duct. This calorimetric probe consists of a straight double-jacketed copper sampling probe, as shown in Fig. (10). The gas sample passes through the center tube, with cooling water flowing around the double jacket to provide both positive tip cooling and facility for calorimetric measurement. Thermocouple 26

emplacements are located at the cooling water inlet and exit ports and also at the gas sample exit port. All required auxiliary equipment is also available, and consists of a portable, self-contained 800 psi coolant system, incorporating the necessary calorimetric instrumentation and a complete, portable gas-sample processing and analysis system. B. Results Axial and radial temperature measurements are made by inserting a probe, connected to a thermocouple, through the sampling ports along the exhaust duct and to a specific radial position. Figure (11) shows that the temperature level obtained from the ceramic exhaust pipe is substantially higher than with the steel pipe. Extrapolation of these results back to the burner indicates that temperatures close to the adiabatic flame temperature are being realized. Some other results showing the variation temperature with distance along the ceramic exhaust duct are shown in Fig. (12) for stoichiometric, lean and rich mixtures. The rapid change in temperature is to be noted and the question of concern here is, of course, whether the chemical kinetics can keep up with this temperature gradient. The emphasis is now being placed on this aspect. Figure (13) shows the temperature at 3 axial stations (1, 2, and 3 ft) in the duct for runs of different fuel-air ratio. As can be seen, 27

the temperature is rather insensitive to mixture ratio at the more downstream locations. 28

RE FERENCES 1. Nicholls, J.A.,Sichel, M., Mason, C.J., Garrison, R.P., Geister, D.E., Liang, H.C., and Paul, J., "Combustion Dynamics as Related to Air Pollution," Report UM 320057-2-P, Sept. 1972, prepared under A.G.A. Project BR-26-1. 2. Bowman, C.T., "Kinetics of Nitric Oxide Formation in Combustion Processes, " XIV Symposium (International) on Combustion, Pennsylvania State University, Aug. 1972. 3. Zeldovich, Ya. B., Sadovnikov, P. Ya., and Frank-Kamenetski, D.A., "Oxidation of Nitrogen in Combustion, " Academy of Sciences of USSR, Inst. of Chemical Physics, Moscow-Leningrad, (trans. by M. Shelef), 1947. 4. Marteney, P. J., "Analytical Study of the Kinetics of Formation of Nitrogen Oxide in Hydrocarbon-Air Combustion, " Comb. Sci. and Tech., 1, 461, 1970. 5. Newhall, H. K. and Shahed, S. M., "Kinetics of Nitric Oxide Formation in High Pressure Flames, " XIII Symposium (International) on Combustion, Combustion Inst., 1971. 6. Lavoie, G.A., Heywood, J.B., and Keck, J.C., "Experimental and Theoretical Study of Nitric Oxide Formation in Internal Combustion Engines, " Comb. Sci. and Tech., 1, 313, 1970. 7. Heywood, J. B., Fay, J.A., and Linden, L.H., "JetAircraft Pollutant Production and Dispersion, " AIAA Paper No. 70-115, 8th Aerospace Sci. Meeting, 1970. 8. Altenkirch, R.A., Shahed, S.M., and Sawyer, R. F., "Nitric Oxide Formation Around Droplets Burning at Elevated Pressures," Comb. Sci. and Tech., 5, 137, 1972. 9. Kremer, H. and Schafer, G., "Verbrennungsverlauf von Methan/ Luft-Gemischen in einer Brennkammer mit heissen Wanden, " Chemi-Ing-Techn., 43, Jan., 1971, Nr. 7. 10. Sarofim, A. F., Williams, G. C., Lambert, N., and Padia, H., "Control of Emission of Nitric Oxide and Carbon Monoxide from Small Scale Combustors," Proc. 2nd Conf. on Natural Gas Research and Technology, June 1972. 29

11. Bowman, C.T., "Investigation of Nitric Oxide Formation Kinetics in Combustion Processes: The Hydrogen-OxygenNitrogen Reaction, " Comb. Sci. and Tech., 3, 37, 1971. 12. Sarofim, A. F. and Pohl, J., "Kinetics of Nitric Oxide Formation in a Premixed Laminar Methane Flame, " XIV Symposium (International) on Combustion, Pennsylvania State University, Aug. 1972. 13. Bray, K.N. C., "Atomic Recombination in a Hypersonic Wind Tunnel Nozzle," J. of Fluid Mech., 6, 132, 1959. 14. Bray, K.N. C., "Chemical and Vibrational Non-equilibrium in Nozzle Flow," Transaction of 9th Symposium on Combustion (1963). 15. Bray, K.N. C., "Chemical and Vibrational Non-equilibrium in Nozzle Flow," Non-Equilibrium Flow, Part II, Edited by P. P. Wegener, Marcel Dekker, New York. 16. Hall, J. G. and Russo, A. L., "Studies of Chemical Nonequilibrium in an Expanding Gas Stream, " Rept. AD 1118-A-6, AFOSR TN 59-1090, Cornell Aeronautical Lab., Buffalo, N.Y., 1959. 17. Cheng,H.K. and Lee, R.S., "Freezing of Dissociation in Supersonic Nozzle Flows, " AIAA J., Paper No. 66-1, Jan. 1966. 18. Cheng, H.K. and Lee, R.S., "Freezing of Dissociation and Recombination in Supersonic Nozzle Flows, " Part I: General Discussion and Special Analysis. Part II: More General Analyses and Numerical Example, AIAA J., Vol. 6, No. 5, May 1968. 19. Fenimore, C.P., "Formation of Nitric Oxide in Premixed Hydrocarbon Flames, "XII Symposium on Comb., Combustion Institute, 1971. 20. JANAF Thermodynamical Tables, The Thermal Research Lab., Dow Chemical Co. 21. Gordon, S. and McBride, S., "Computer Program for Calculation of Complex Chemical Equilibrium Composition, Rocket Performance Incident and Reflected Shock and C-J Detonations, " NASA SP-273, Lewis Research Center, Cleveland, Ohio, 1971. 30

22. Bowen, S.W., AMPLCT, A Numerical Integration Routine for System of Stiff Differential Equations. University of Michigan Office of Research Adminstration, Report 033390-2-7 June 1971. 23. D'Souza, M.V. and Karim, G.A., "An Analytical Study of Methane Oxidation in a Steady-Flow Reactor, " Combustion Science and Technology, Vol. 3, pp. 83-89, 1970.

APPENDIX A Program Used for the Detailed Kinetic Model of H2-02 Mixture The program consists of a main program,two subroutines and block data. A description of the function of each of these is given below: (1) Main Program The program performs the following functions (a) Read in and check the input data. (b) Call Subroutine STIFF3 to initiate the integration. (c) Further control the step size in order to obtain a convergent and stable solution in the course of integration through induction period, main reaction, and radical recombination regions. (d) Print out the calculated results. (2) Subroutine STIFF3 This Subroutine locally linearizes the system of differential equations and converts them into a set of algebraic equations. The algebraic equations in matrix form are triangularized by Gaussian elimination procedures and solved by the method of linear algebra. An implicit technique was used to obtain an accurate and stable solution for all negative eigenvalues of the 32

system of linear, algebraic equations. The Jacobians of the right hand sides of the differential equations are evaluated numerically with the subroutine DER. (3) Subroutine DER This subroutine sets up expressions for the evaluation of the following: (a) The rate and equilibrium constants for all chemical reactions. (b) The species enthalpy, Hi, and Cpi, i.e. dHi/dT. (c) Thermodynamic properties: density, molecular weight and total enthalpy of the gas mixture. (d) The time derivatives of all dependent variables, In Zi, and In (1/T). (4) Block Data Block data contains (a) Molecular weight of all species. (b) Constants of Arrhenius equation for the evaluation of rate and equilibrium constants. (c) Coefficients of polynomials for the computation of the species enthalpy. (5) The Print Out of the Complete Program: 33

I C MAIN PHIOGilAN OF KINETIC MODEL > C iI2 —)2 lIXXTUFhL 3 > P=}iESSUEE (ATlNOSPHEI E) 4 C T-TLPEliATUILE (DLEGREE K) >'CT G TICiEI(,TCITL;A2=l PAInAMETERS FOR CONTROLLING STEP SIEZE 6 X=TIME (SL0ND) 7 Z(I)=N0. 0F MOLLS OF SPECIES I PER UNIT WEIGHT OF MIXT. > 6 SUMh=TUTAL MOILES OF ATOMIC H PER UNIT WEIGHT OF MIXT. 9 C SUbNO=TOTAL MOLES -F ATOMIC 0 PER UNIT WEIGHT OF MIXT. i 0 WI'NI=NO'LEGULAiR WIEI GHT OF MIXTURE(GM/MOLE) 1 C RHC=DENSITY OF MIXTURE(GM/CM3) 112 M (I)=M OLECULAR WEIGHiT 0F ITH SPECIES 1 13 C Y( I )=DLOG(Z(I) ) I=1,4 14 C Y(1,5)=-DLOG(T) 15 C I1 IMPLICIT REAL*8(A-H,0-Z) 17 REAL*8 Y(8,5),DY(5),M(6)JZ(6),A(6,2,6),RC(3o2,8), 1 1 U (2, 5), DELY(5), MWM 19 INTEGER KFLAG O0 DIMENSION XM(6),L(5),PW(5,5),QMAT(5,5) 21 COMMON PSUM0,SUMH, HMDELY 22 COMMON / BLDATA / MA,RC 23 EXTERNAL DER 24 C 25 C iREAD AND CHECK TiE INPUT DATA,, C 27 1 READ(5,100)P, TIN IT,XINIT,XF INAL,H,HMAX1,HMAX2,YMIN, 2 1 EPS1,L PS2, iRTNJ0,TCHEK1, (XN( I ), I= i, 6) > 9 100 FORMAT('',4D18. 10,/4L;18. 10/4D18. 10, O0 1 /4D16. 10,/2iL18. 10) 31 WRITE(6,101)P,TINITXINITXFINAL, II'MAX IHMAX2, 32 1 YlINI8JEPS1 EPS2,PRTNOT:CHEK1 33 101 FORMAT(//,'PRINJT P, TINIT,XINITXFINAL',/4D13-5, > 34 1 //,'PRINT HHMAX1,HMAX2,YrlIIN',/4DI13.5 3:5 2 //,'PRINT EPS1,EPS2, PRTN0,TCHEII' /4L13.5) 37 t CALCULATE MWM, RHOAND Y ( 1,I),I=1, 5 3& 8 39 X=XIN T 40 XCIiEK=0.0 41 HMAX=HMAX1 /42 EPS=EPSI 45 TCHEK2=TIN IT+0. DO 44 N:WM=0. 0 > 4 5 $ U=0 ~ 0 4ow:.:L 2 I=1,6 47;2 It.'di4=Ml.T+XN~ ( I ) *MXi ( I ) > d LG 3 I=1,6 9 3 I ( I )=>;M I C/-jM,.9 3 Z()'=XM(I)/M!WM ~ 5C 4 I=l,4 51 DY( I )=0-0 > a Y( 1, I )=DLG (Z( I ) ) 54 U( 2, I)=Y( 1, I) 55 DELY( I ) =0. 5D-2,6 SUMtSA=SUMA+Z ( I ) 57 4 CONTINUE LC T=TINJIT SL Y( 1,5 )-DLOG(T) 60 U(,5)=Y(1,5) 61 U(2,5)Y(1,5) 34

>.'' b LLdL (LbY)=O.bu-2_ ~> 3 DY(5)=0.O > 64 SUM O=SUMt,-Z ( I ) +2. D0*Z (5) 0o SUNHti=SUirA-Z ( 2) +Z ( 3 ) + 2. DO*Z ( 6 ) 06to H0i;3=,; =P*NI/ ( b62. O05*T ) 07 CALL DEtli(X,Y, bY) oo WiWITL (6, 102)X, T, 9 IUiti,H, (XM(I), I=1, 6 ), It;i > 6"'? 10 2 FORMlAT ( //bX,'T I M', 5, T 1IPET'ATUL',6X, 6,'Mt', OX, 70> 1'IriC', 10X,'1i',/5X,'XM( 1 )' 7X,'XM(2)', 8X 71 2'XM(3)',bX,'XM(4)' 7Xv'XM(5)'/5X,'XM(6)',8X' iM;'. //, 5D1 3 5 J /5D1 3 5. /2D 1 3. 5) 7> /3 1 1 CALL STIFF3( 5,X Y, DY, DER, H, EPS HMAX,YMINQMAT, > 7, 1 L, PW, KFLAG) > 75 DO 41 1=1,5 76 U(2, I)=U(I 1I) > 77 U(1.,i)=Y( 1,I) > 78 DELY(I )=DABS( (U( 1, I )-U (2, I ) )/(U( 1, I )+U(2, I) ) ) > 79 I (ELY ( I ).LE. D- 10) DELY(I)=I. D-10 > 60 41 C0NTINUL > 61 XCHEK=XCHEK+1 DO > k82 IF(KFLAG.NE.*) G0 T0 21 > b;3 IF(XCHEK.LT.PRTN0) G0 T0 11 84 XCHEK=0. 0 > 85 SYM=O.0 > 86 DO 13 I=1,4 > 87 Z(I)=DEXP(Y(1,I)) > 88 SYI-=SYE1+Z ( I) > 89 13 CONTINUE 90 Z( 5 ) = ( SUM - SY+Z( 1) ) / 2. O 91> i Z (6 ) = ( SUI31I{- SYM1+Z ( 2 ) -Z ( 3) ) /2. DO > 92 M;tt<= 1. DO/( SY~l.+Z (5 ) +Z (6) ) > 93 F=IiG=P*,MWM/(82.05*DlEXP(-Y(1, 5))) > c94 DO 14 1=1,6 95 14 X1 II ) =Z(I)*MWM > 96 CALL DER(X.Y, DY) > 97 T=DEXP(-Y(1,5)) > 98 IF(T.LT.TCHEK2) GO T0 31 > 99 IF(T.GT.TCHEK1) G0 TO 30 > 1(00 EPS=EPS2 > 101 HMAX=HMAX 2 > 102 G0 TO 31 > 103 30 EPS=2.DO*EPS > 104 HMAX= 2 DO*HiRAX > 105 IF(!M1AX.GT. tiMiAX1) HiiAX=HMAX1 > 106 IF(EPS.GT.EPS1) EPS=EPS1 > 107 31:RFITE(6.103)X, T.bW.,F.,H,(XM(I),.I=1,6),HI. > 108 103 FORMAT(/5SDI3.5,/5D13.5,/2D13.5,/) > 109 IF(X-LT*XFINAL) GO TO 11 > ii10 GO TO 1 > 111 21 WRITE(6. 121) > 1 C121 FORFiAT ( / / KFLAG= I3) > 113S G0 T3 1 > 114 END;2JD 0F FILL. 35

SUbROUT INE ST IFf 3 (N JT Y. DVVDEF(r, HEP S'tAX, I YNINQMAT.,L,PW, KFLAG) 3 C 4 C N=NUMLER OF AUTOMOTRPHICI DIFFERNTIAL EOUATIO4NS 5 C T= INDEPENDENT VARIABLE 6 C Y= 8 BY N AtiiAY 0F WHIChi NLY FIRiST ROW ( 1,1) T0 ( 1,N) 7 C IS EXTERNALLY USED. FIRST R0OW CONTAINS DEPEtJDENT 6 C VACIAiAbLE 9 C DY=DERIVATIVE VECTORj LENGTH N, SUPPLIED FROM DER l O C DE = SUb i,0li UTINE HAVING AR GUMENTS(T. YJ DY) THAT EVALUATES 1 1 C VECTOR 12 C H=INITIAL AND CURRENT STEPSIZE 13 C EPS=MAXINUM ALLOWED RELATIVE CHANGE FOR ANY DELTA Y/Y 14 C Yl'IN= MINIMUM VALUE OF ANY DEPENDENT VARIABLE CONSIDERED 15 IN FINDING DELTA Y/Y 16 C IFLAG =1 NORMAL RETURN 17 C =-4 SINGULAR MATRIX 18 C QIAT N BY N ARRAY 19 C L= INTEGER WORK VECTOR, DIMENSION N 20 C PW= JACOBIAN MATRIXoN BY N 21 C PWL.,Y DY AND QMAT MUST BE EXPLICITLY DIMENSIONED 22 C FROM THE CALLING PR0GRAM 23 C DELY( I )=DABS ( U( 1 I )-U(2, I) )/(U ( 1, )+U(2, ) )) 24 IMPLICIT REAL*8(A-H, -Z) 25 REAL*8 Y(8,5) PW(5 5),DY(5),DELY(5) 26 DIMENSION QMlAT(5,5),L(5) 27 COMMON P. SUNO, SUMH, HM, DELY 26 EXTERNAL DER 29 vKFLAG=1 30 CALL DER(T,Y, DY) 31 DO 1 JEQ=I,N 32 Y(2,JZQ)=Y(1,JEQ) 33 1 Y(5,JEQ)-DY(JEQ) 34 DO 2 JTERM=IN 35 DC 9 JT=1,N. 36 9 Y(1,JT)=Y(2oJT) 37 Y( 1 JTERM)( I. DO+DELY(JTERM) )*Y(2,JTERM) + 1 ~ D-30 38 CALL DER(T, YDY) 39 DO 3 JEQ=1,N 40 3 Y(6,JEQ) =DY(JEQ) 41 D0 109 JT-=,N.~L 109 Y(1,JT)-Y(2,JT) 43j Y ( 1YIAJTER, )- ( =( D0-DELY(JTERM) )*Y(2,JTERM) - I.D-30 "i -, CALL LEnr; (T,Y, DY) DO 103 J(C=1,N 4, 103 Y(7,J&LQ)-=~Y(JEQ) 4r/, C Y(bt.JEC) CZNTAINS F AT ( 1.D0+DELY(JLEQ))*Y C Y(7,jLE) CONTAINS F AT ( 1.D0-DELY(JEQ))*Y;,' LD 5 JEQJIN 5 PW (J E, J T hRN) = (Y(6 dJEQ )-Y( 7,JEQ ) ) 5. 1 / ( 2. DO*DELY(JTERM ) *Y (2,JTERM) +2. D- 30) 2 CONTINUE 53 JUMAT= 1 5, -4 20 ID 17 IC=I1,N 55 DO 17 IRU',lq 5,5 QMAT IR, IC) =-.SD0*H*PW ZR I! C) 57 i ( IRES0. IC) QMAT( IR, ZC)uQMAT(X IR IC)+1 DO 5~ ~17 C0NTINUE 59 CALL DLCIOMP(N,N, QMATJL) 60 IF(L(N). N.0O) G0 T0 22 61 WRITE(6r,100) 36

~> o- I u O I. I?::AT(I' ((I)) - ((A))/. IS SII.C;I:tTLi LA') t~> 03 1(F LA G = - 4;> ozl 4GO TO 23 6> b o22 CONTINUE > (0o C HEUSE DY VECTOR o7 DO 7 JEQ=1,N >, 65 7 DY(JEdQ)=1i*Y(5,JEQ) > 69' CALL SOLVE (N, N, QMAT, DY, L ) > 70 C DY NOW CONTAINS DELTA Y VECTOR > 7 1 JMAT=JNAT + 1 > 7 IF(JMAT-10) 21, 19, 19 >,5 c21 Abl=0. D0 ~> 7,~ VDO 10 JtLO1,N > 7z IF(DAiS(Y(2jJEQ)).LE.YMIN) GO TO 10 7> o ANA=LA S(DY( JEQ ) /Y(2,JE ) ) 77> 7 AM =DMAX 1 ( AbM, AMA ) > 7z 10 CONTINUE 79 IF(AM.EQ.0.DO0) GO TO 19, 90 IF(AM-EPS) 1, 1 1.12 -> w1 11 IF (AM-2. D- I *EPS) 18, 1 8,'1 9 > $2 12 H=5.1)-1.H GO TO 20 > 84 18 CONTINUE ~> ~: IF(Hi.GE.IHMAX) GO T0 19 6> 6 H=2. DO** 87 IF(H.GT.HiMAX) H=HMAX G0 TO 20 >'9 19 DO 14 JEQ=I,N 90 14 Y( 1,JEQ)=Y(2,JEQ)+DY(JEQ) > 91 T=T+H > 92 23 CONTINUE 93 fiETURN > 94 END.> 95 SUBROUTINE DECOMP(N, NDIM, A, IP) 96 REAL*8 A(NDIM,NDIM),T,DABS > 7 INTEGER IP(NDIM) 98 C 99 C MATRI.X TRIANGULARIZATION BY GAUSSIAN ELIMINATION. I 00 uG C INPUT.. 1> i1 C N = ORDER 0F MATRIX. 1> 02 C NDIM = DECLARED DIMENSION OF ARRAY A > 103 C A = MATRIX TO BE TRIANGULARIZED. > 1C4 C OUTPUT.. 105 C A(I,J), I.LE.J = UPPER TRIANGULAR FACTR, U'J > 106 C. A(I.J), I.GT.J = MULTIPLIERS = LOWER TRIANGUILAR > 07 C FAC TO,01R I-L > iC> C IP(K), K.LT.N = INDEX OF I<-THI PIVOT ROW. 109 C 1P(N) = (-1)**(NUMNBER OF INTERCHANGES) OR 0 S 11C C USE'SOLVE' T6 0BTAIN SOLUTION OF LINEAR SYSTEM!. > 1 1 1 G LETERDM(A) = IP(N)*A( 1 1 )*A(2,2)*...*A(N, ) > ii,1. IF I?(N)=0, A IS SINGULAR,. SOLVE WILL DIVIDE BY EF:,. > I 1 CU IrJTEhCHANGES FINISHED IN U 0ONLY PARTLY IN L 1> 1 1_ IP(N) = 1 > 11it) DO 6 K = 1,N - 1 i7 IF(K.EO.N) GO T0 5 > 1,i KPI = I<+1:> 111 b = K 1, DX8O 1 I = KP1,N' 1; - IF(DAIS(A(I,I)).GT.DABS(A(M,K))) r- = I > 1.1'CNIT INUE > 1:.: IP(K) ='! 1 I4 ~IF(Ii.NL.1K) IP(N)7= -IP(N) 37

>:25 T = A (M., ) > 126 A(M.,{) = A( KK) > 1L7 A(i, ]) = T ~>'"U''sX (T EIF(T.LQ.0.DO) GO TO 5 > 1"9 DO 2' I = KP1,N > 13 2 A(I,KI) = -A(I, )/T > 131 DO 4~ J = I<P1,N > 13'2 T = A(lOJ) > I'JJ A(bljJ) = A(i, J) > 1.'A ( ]{, J ) = T > 135 IF(T.LQ.O.D0) GO TO 4 > 13 6 DO 3 I = HlP1,N > 1;37 3 A IJ) A( I.J) + A(I,K)*T > 136 4 CONTINUE > 1,}9 5 IF(A(hK,K).EQ.O.DO) IP(N) 0 > 143 6 CONTINUE > 141 RETURN > 1'42 END > 143 C > 14 C > 1!.5 SULROUTINE SOLVE(N, NDIMj A, B, IP) > 146 REAL*8 A(NDIMNDIM),B(NDIM),T > 147 INTEGER IP(NDIM) > 1i. 4 > 1 4') C SOLUTION 0F LINEAR SYSTEM, A*X = B > i _ 1 C INPUT.. > 151 N = 0RDEJ. 0OF MATRiIX. >!15. C NDIM = DECLARED DIMENSION OF AR.AY A i153 C A = TRIANGULARIZ'ED MATRIX 0DTAINED FF0OM'DECOMP'. > 1'4 C L3 = RIGHtT HAND SIDE VECTOR. > 1$5 C IP = PIVOT VECTOR OBTAINED FROM'DECOMP'. > 156 G DO NOT USE IF DECAMP HAS SET IP(N)=O 0 > 1-.7 C OUTPUT.. > 1i:_. C B = SOLUTION VECTOR, X e > _'J60 Ir(N.EQ.1) GO TO 9 > 1 61 NMI = N- 1;> i 2 D0 7 K = lNMI > 163 KP1 = K+1 > 164 MV = IP(K) > 165 T = B(M) > ioc (M) = B(K) > i. 7 b(K) = T >':'.. DC 7 1 = KPl,N:>.';~?7 B(I) = SI( ) + A( I,i)*T >;7? DO $ Kr 7 = 1,NDIl > 17i l{1l = N -i13 > 1 7J BCK() = 1)(KC)/A(I,= K) ~ 1 74i T = -B(]() > 1'5 DO 8 I -=,:t.I1 > 176 8& 1(I) = i3(I) + A(r, I)*T > 177 9 D(1) = S(1)/A(1,I) > 1 76 RETURN > 179 END _.' J Z FILE 38

> C - -.....-SUthibIJT IN L DRi(X,Y, DY) --------- > 3 SUbROUTINE DER(X,Y, DY) >'4 C > 5 iMAPLICIT REAL*, (A-H i,0-Z) > o REAL*8 Y(8,5),DY(5). MWM,11(6), RC(3, 2, 8), RK(3, 8), 7 1 ZZ(6).,i(6),CP(6), TCP(4),TH(4) 8 DIMENSION L(5),PW(5,5),QQMAT(5,5),A(6, 2,6),YK(6) 9 CINIMON P, SUM0.,SUMH,HNlr 10 COMMON / 6ILDATA / M.A,RC > 1 C > 12 C EVALUATION OF MWM & RHO OF GAS MIXTURE 13 C 14 SY1=O. 0 > 1D 1 1 =1,4 16 Z(!)=D2EXP(Y(1,I)) > 717 SYMl=SYMsI+Z( I) > 1 CONTINUE > 19 Z (5) =(SUO-(3SYM+Z( I ) )/2.DO ~> ^O Z(6)=(SUMH-SYM+Z(2)-Z(3) )/2.DO > 21 MWM= 1. DO/(SYM+Z (5) +Z (6) ) 2> 2 T=DEXP(-Y ( 1, 5 ) ) 23 CD1=P/(82.05*T) 24 RH=CD 1 *MWM 25 C:= Zo C EVALUATION OF RATE CONSTANTS > 7 C 6>,; DO 3 I=1,8 > 29 DO 2 J=1;,2 > 30 2 RI(J, I ) =1C( I, I)*(T**RC (2,J,I))*DEXP(RC (3, JI)/T) > 3 1 IF(I.GE.3.0) GO TO I I > 32' HK( 3, I )=-RK( 1, 1 )/RtEK(2, I ) 3> 3 GO TO 3 3 34 11 K(3, I )=RX 1, I )*RK( 2, I ) > 35 3 CONTINUE' o C 37 C EVALUATION OF CHEMICAL RATE OF FORMATION OF SPECIES 36 C > 39 Y.1=K( 1,i )*Z( 1 )*Z(5 )-F.K(3, 1 )*Z(4)*Z(2) > 40 Y2=RK( 1,2)*Z(2)*Z(6)-RK(3,2)*Z(4)*Z( ) 41 Y3=RK( 13)*Z( 1 )*Z(3)-RK(3,3)*Z(6)*Z(4) 42 Y4=RK( 1,4)*Z(2)*Z(3)-RK(3,4)*Z(4)*Z(4) ~ 43 Y5=(RK(.S)*Z(i)*Z(*Z1 )-RK(3,5)*Z6)/RH0)*CDI > 44 Y6=(RK(1,6)*Z(2)*Z(i4)-RK(3,6)*Z(3)/RH0)*CD1 > 45 Y7=(RK(l 7)*Z(2)*Z(1)-RK(3. 7)*Z(4)/RH0)*CDI 46 Y8=(RK( 18)*Z(2)*Z(2)-RK(3, 8)*Z(5)/RH)*CD1 47 C > 46 C EVALUATION OF TIME RATE OF CHANGE OF SPECIES CQNCEZNrT"AT 49 C 5G YK( 1 ) =(-Y+Y2-Y3-2. DO*Y5-Y7)*RH 5 1 YK(2)= ( Y 1 -.Y2-Y4-Y6-Y7 -2. DO*Y8 )*RHO > S YK( 3 ) = ( -Y-Y4z+Y 6)* Ril 5 5 Y ( 4 ) = ( Y 1 +Y2+Y3+2. DO*Y4-Y6+Y7) *RH0 > 54 Y I( 5) = ( -Y1+Y8)*RIl0 ~ ~ YK((6)=(-Ya2+Y3+Y5>)*Ri0 56> D0 4 1=1,4 > 57 4 DY(I ) =YK ( I )/Z ( I) 58 C > 55 C EVALUATION OF TIME RATE OF CHANGE OF TEMbPERATURE,'( l,5) > C60 C > Iitl= O. O >., 3'a GI'CPI1= 0.~0 > 6_5 Q[C:0.0 > t~ IF (T.LT. l.D 3) GO TO 5 39

> o05 K=1 > 66 G0 TO 6 > 67 5 K=2 > 68 6 DO 8 I=1,6 > 69 CP(I)=A( 1K, I) > 70 H(I)=A(6,KjI)+A(1,K,I)*T > 71 DO 7 J=2,5 > 72 LI=J-1 > 73 TCP (J) =T**L 1 > 7.i TH(J)=TCP(J)*T/J > 7 CP( I )=CP( I )+A (J, I )*TCP (J) > 76 H( I )=H( I) +A(JK, I )*TH(J) > 77 7 CNiTINUiE CPl=CP+Z. ( I )*CP( ) > 79 HN=HM+Z (I )*li( I) > SO QM=QM+YK( I)*H(I) > ol 8 CONTINUE > 62 DY (5) =QM/(CPM*T) > o3 HM-i=HM* 1. 986 >.84 RETURN > &5 END,kEND 0F FILE 40

i L CbLOCK DATA 3 L; $SUrTibP:0C;;iM FOI DATA: M(I1),A(IJ1,K I, II).C(J2, K2, I2),4 C 5 C II=IDENTIFICATION NO. OF SPECIES FOR H,O H20,oH,0I2, u C & }I2 1;!i: R; PEC;TIVELY 7 k; I2=IDLNTIFICATION NO. OF CHIEMICAL REACTIONS > C JI=JITH COEFF. FOIR EVALUATION OF lf & CP 9 C J2=JSTHI CONTANT FOR ARRHENIUS EQtIIATION OF RATE AND 10 C EQUILIbRIUM GUiTANITS 1 1 G K1=INDEX F0R COEFF FT;I( EVALUATIOIN OF H & CP 1 C KI1=: T.GR. 10000 DEGREE K > 13 C 1<1=2: T.LE. 1000 DEGREE t4 > 1 4 G KS2=INDEX FOR RATE AND EQUILIbRIUM CONSTANTS 15 C; 1(2=1: FOR F0,W!ARD RATE CONSTANTS c.'. (2=2 FOR EQUILIBRIU CONSTiANT 17 C 1(2=3: F 10: bACKWARD RATE GONSTAINT 19 }IEAL*8 M,A, RCAA 20 DIMENSION M(6),A(6,2,6).RC(3,2r8),AA(6,2,2) 21 COMMON / BLDATA / M,A,I RC 2o EQUIVALENCE(A (, 1,5),AA ( 1, I)) > i i C C MOLECULAR WEIGHT OF SPECIES, C 27 DATA N/ 9 1 +0. 10079700D+01, +0.15999400D+02,+0. 18015340L+02, 30 2 +0. 17007370D+02.+0.31998800D+02+00.20159400D+01/ 3i J 3'~ C COEFF. FOR EVALUATION OF H & CP OF SPECIES > ~ DATA A/ j: G....-C CO.. 0EFF. FOR H (I1=1) 27 1 +0.25000000D+01, +000 OOOOOO0D+0O, +0. OOCOOOCOD+OO, 3o 2 +0.OOOOOOO00 +00+0. 00000000GD+0C,+0.25471627D+05, 39 3 +0. 25000000D+01, +O. OOOOOOOOD+00, +0. 00 CU(OJOD+00, > 64& k.0 4 +0. OOOOOOOOD+00+0. OOOOOOOOD+00 +0. 25471627D+05, i1 C --- COEFF. FOR 0 (11=2) > - 1 +0. 25420596D+01, -0. 275506 19D-04, -0.31!128 033D-C,,> 2 +0.4 1 674- 1 -O. 4306S5 1 5D- 1 5 +0.2 3C,S31+05 > 4 3 +0.29464287D+01, — 0. 16351665D-09 +0.2A 10316D-05, > 4 0. 16028432D 08, +0. 38906964D- 1 2 +0. 2' h1476t44D+05, 4 C.... —-.... COEFF. FOR H'20 (I1=3)?~7 1 stO. 2 7167633D+0 1 ~ +0. 2945 1374D-0'C, -O 80 C2243741D-0 6, > 4 +: -2 0. 1 O2G -G9 - -0.4 l 72145 1, - 0.2 C! ~D-t r+,? > ~;3 +0.40701 275D+0 1-0 -. I 1 084991S-C2,+('..!4: 1.'1! 0D-5 5,t 4 -0. 29637404D-08, +0. 8070'21C0 3l D- 1, -0. 3(' 79'7?c- D+, C.- --------- (COEi F. FOR OFH (-I I =4) -- 35 1 +0.29106427D+0 1 +0. 959316':''>'. - 3 3, -0C I L9.. w17C 2-, 53 2 +0.13756646D-10, +0. 14224542D- 15, +.0,535 1iL.4 ~4;3 +0.38375943+0 1 -0 10776858D-02, +0.96830378D-06 > 5 4 +0.18713972D-09,-0.22571094D-12, +0.36412823D+04/ 26 C 57 DATA AA/ 59 C.......... COEFF. FO; 02 (I1=5) -. --. 60 1 + O 36 21 9535D+,+ 1 + O. 7 3 6 1 8 6 I:D- 0 3 - O. 19 o 5 2' L - C, 61 2 +0. 36210 5 58D-! 10,-0. 2894z5627D-1 4 -0 14 v 9 -. I " 1 ":+C 4, 3 3 + 3 5f5985D+0-0 1 + 0 9 +87821 84DD-02 + 0. C 7C5454z:D-. 63 4 -0.676351 37D-0,+0 21555993D-11. J 10473226+.;4, 64 G.-.-...... CEFF. F0R H2 (1=6)....... 41

1 t+ 0. o 1 ) (C1 9 3 1 9L + () 1 + 0. 5 1 1 1 9 1 4 b - O., +O. 6 5 1 /4/1 2 ) L -, 7, t -,- -0. 3Z9C)'99 7 L- 1 O49+ 997D-10,+0. 3695/;5 j- 1/, -(). r773 r()42L,+0 j;, +O.3 b 7 45 1;D+0+00. 2676 5 200bO2,-0. t58 099162DC-05,. () ~t 4 O * 5b5I0 1 039 1 L)- O 8r - 0. 1 8 1 22739 D- 1 1,-0.93 904 74D+0 3/ 70 W) C COUNSTANT) F0R EVALUATI0N 0F RATE AND ErOIJILIBRIIUM > 71 o UGNSTANT. C DATA i'/ l 7. C ---— _........ I i+=l -+- (1 2=1). - -..... 1 LI1 +0.22000000D+15,+0000t)00000,+00,+0*83400000D+0, +0.3770603D+03, -0. 399330COD+00, +. 866 14309D+04, o C - - - - - - -+L2=Ufli+Hll (I2=2) - - 7'0 1 +0. 17000000L'+14,+0.O000000b+0,+0.47500000D+04, 2 +0. 1 72tb924DL +O 1, +0* 33668377LD-0 1, +0. 91 140673D+03, i 1 C........ —-- +H 1-;2( = H 1+0- H ( I 2 = -3) b I +0.84000000D+14,+0O.OOOOOOOOD+00,. +010100000D+05 Ob 2 +0. 2279984o D-01, +0. 2755986 7D+00, -0. 7908496 7D+04. 84 C ---------- 0+HI20=20}1- (I2=4) - 1 +0. 58000000D+ 14, +0. ~ 0000OO D+00, +.O 90 700000D+04, 2 +0. 1 3 1 43D-O 1, +O.24 1 99278DL+0. -O. 88 1 99665D+(4., U...-. -... -- + H+ +1=i2+M (I2=5) - - B —rsc I ~+0.10000000D+19,-O. 1o0000O00L+G1+0,+O.OOOCOOD+O, 2 +0.36 77904 1DtO +0. 436 1 0090;-02, +0. 52393204D+05., 90 --.-....... H+0I+["i=H20+M (I2=6) 91 1 +0.40000000+20,-0. 10000000D+01,+0.0000 OOOOOOOD+00'-,z c2 +0. 161 3 8405D+ 03-O. 27 1 2924 2 D + 0,+0.6302054Dt+05 93 C ---------- O0+H+M=0H+M (I2=7) - 94 1 +0 53000000D+ 16,+0.OOOOOOOOD+00 +0. 27800000D+G04 Lu 2 +0.21271 777D+O 1,-0. 29299637D-O 1 +0.5 1482088D+05, 96 C -------- O+a+M=02+M (I2=8)97 1 O +0.4000000OdD+18-0.IO 1OOO000D+01+0.00000OOOOOOOD+00O'9:)~ 2 +0.80164625D+03,-0.42867870D+00,+0 60143569D+05/ 99 C i O END END; I;.' F ILL 42

TABLE 1. RATE CONSTANTS OF NITRIC OXIDE MECHANISM (I) O+N2 = NO+N 14 0.09 -37750/T klf =0.627 x 10 T e 14 -167/T klb = 0.31 x 10 e lb (II) N +0 2 = NO + O k2f = 6.43 x 109 T e-3125/T k2b =0. 506 x 109 T115 e -9000/T (III) N + OH - NO + H k3 =4.2x 1013 k =1.6 1014 - -23, 850/T The above rate constants are the same as those used by Sarofim, et. a I

TABLE 2. CHEMICAL MODEL NUMERICAL SYSTEM OF H2-02 MIXTURE i Species 1 H 2 O 3 H20 4 OH 5 02 6 H2 j Reaction 1 H + 02 = OH + O 2 O + H2 = OH + H 3 H + H20=H2 +OH 4 O + H20 = 20H 5 H ++M =H2 +M 6 O+OH+M=H20 +M 7 O+H+M=OH+M 8 0 O+O+M =02 +M

TABLE 1. RATE CONSTANTS OF NITRIC OXIDE MECHANISM (I) O+N2 = NO+N 14 0.09 -37750/T klf =0.627 x 10 T e 14 -167/T klb = 0.31 x 10 e lb (II) N +0 2 = NO + O k2f = 6.43 x 109 T e-3125/T k2b =0. 506 x 109 T115 e -9000/T (III) N + OH - NO + H k3 =4.2x 1013 k =1.6 1014 - -23, 850/T The above rate constants are the same as those used by Sarofim, et. a I

TABLE 4. EQUILIBRIUM CONSTANTS OF H2-02 SYSTEM C B R.i =AT e Reaction No. A B C (J) 1 0. 3767 x 103 -0. 39933 0. 86614 x 104 3 0.4386 x102 -0.2756 0.'7908 5 x 104 4 0. 75872 x 102 -0.242 0. 882 x 104 5 0.27174 -0. 4361 x 10'2 -0. 52393 x 105 6 0. 61966 x 10-2 0.27129 -0.60302 x 105 7 0.8915 0.293 x 10-1 -0. 51482 x 105 8 0. 12474 x 10-2 0. 42868 -0. 60143 x 105

TABLE 5. COEFFICIENTS FOR EVALUATION OF MOLAR ENTHALPY OF SPECIES Species Temp. A A A A A3 AA Range1 2 3 4 A H 300- 100 K 0.25E1 0 0 0 0 0. 25471627E5 i=l 1000- 5000 K 0.25E1 0 0 0 0 0. 25471627E5 O 300- 1000 K 0.25420596E1 -0.27550619E-4 -0. 31028033E-8 0. 45510674E-11 -0. 43680515E-15 0. 29230803E5 i=2 1000- 50000K 0.29464287E1 -0. 16331665E-2 0.24210316E-5 -0.16028432E-8 0. 38906964E-12 0. 29147644E5 H20 300- 1000 K 0.27167633E1 0.29451374E-2 -0.80224374E-6 0. 10226682E-9 -0. 48472145E-14 -0. 29905826E5 i=3 1000- 50000K 0.40701275E1 -0. 11084499E-2 0. 41521180E-5 -0.29637404E-8 0. 80702103E-12 -0. 30279722E:5 OH 300- 1000~K 0.29106427E1 0.95931650E-3 -0.19441702E-6 0.13756646E-10 0. 14224542E-15 0. 39353315E4 i=4 1000- 50000K 0.38375943E1 -0.10778853E-2 0.96830373E-6 0. 18713972E-9 -0.22571094E-12 0.36412323E4 02 300- 10000K 0. 36219535E1 0. 73618264E-3 -0. 19652228E-6 0. 36201558E-10 -0. 28945627E-14 -0. 12019825E4 i = 5 1000- 5000 K 0.36255985E1 -0. 18782184E-2 0. 70554544E-5 -0.67635137E-8 0. 21555993E-11 -0. 10475226E4 H2 300- 10000K 0. 31001901E1 0. 51119464E-3 0. 52644210E-7 j -0. 34909973E-10 0. 36945345E-14 -0. 87738042E3 i=6 1000-5000 K 0. 30574451E1 0.26765200E-2 -0.58099162E-5 0. 55210391E-8 -0. 18122739E-11 -0. 98890474E3 2 3 4 5 T T T T Hi =T +A A A + A - A 1 1 2 253 5 46

x( V V 0 C Fuel r~, v,. T Tc Air - -_ Mixture =...7 Tr, x) Inlet Fully Developed Burner Section Flow Figure 1. Slow in the Exhaust System.

(4) 0. 81 0 00000 z 0. 4 0.22 0 ~0. 1 0.2 0. 3 0.4 own0 6o?0,80 0.4 ~~~~~~~~~~~~~~~~ ~~~~~~~~0.5 0. 6 0.7.. time, t, (second) Figure 2. versus t Profiles for Various SimpljJied Equilibrium Models 4-Air Mixture at 0 = l. O p 1 atm and Tc = 2227OX. (1) Newhall's Model (2) Marteney s Model (3) Zeldovichts Model (4) Zelvovichts Model with B >> 1 (5) 2-reaction Model (6) 3-reaction Model

1.0 1.0 XNO XNO 0. 8 NO NO X~~~~~o ip=~~~~~~~0.6 1. XNO | 1. 6 /2 and 3 reaction,= - 1.6 model 3 reaction model 0.4 — — 0.4 = 1.6 2 reaction model 0.2 -- 0.2 10103 10103 04 t, Time (sec) Figure 3. XN versus t Profile for CH -Air Mixture at -= 0. 6 and 1. 6 reaction and 3 reaction model) p - 1 atm.

105 104 T = 298 0K reactants 10 10 Treactants'10000NO (sec) 102 N-O 10 10-1 10o-2 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.08 2..0 2.2 2- 4 Equivalence Ratio, p Figure 4. Characteristic Time, TNO of NO Formation Versus Equivalence Ratio, q, Profiles for CH4-Air Mixture at 1 atm and TO = 298 and 1000~K.

~IM0001 put -I 86-'9186g:= sju ~r'tula Tr = d 41 aonlxjt Z9N-'O'- 3 IoJ "o'o'OlatI aouaATnb[ snsJ1A UOTI1LtU10, ON JO ~'oauaT.,L o!IsTjalo@Bly D'~ aTnS2Tf (ot e aouaoIA.nbA)'P"'Z Z o'Z 0Z 81 9'. 1 T T O'T. 8'0 90' O'o o 01 ~_o~' I I i I.I \\ I O1'-''01-;_?, r i "T O z z0

l~o r.256 C1aicuaated i by (.942 3 1. 6 (dT/dt) = 65 for x/vC 12- 2 X O dT 10 0. 2x10~~~~~~~~~~~~~ dt LN v 628 ~~~~~~~~~~~~~~~~~~~~~~~~~ C 314 ~2.(dT/dt) fo;unr 1 ~~~Nd iidT/ml ~~con stant o ii. xiv2 2x 10N t o~o. T / d t (d T2 =0.157 ~~~~~~~~~~dt - 5 (.0656 (.(344 x/ve019 0.4 c c ~O 0~~~~~~~~~~~~1,~2 0 /se )ue ~ a~ 5 to~~~~~~~~~~~~~~~~~~~~~~~~~-Td( l I Mstur at ~ e %NO Ver ~s d3/dt CrVes at VartoUs

1.0 dT dT dT dT = - a = constant dt dx 0.8 VT ~0. 80~ T 1-)l/2 V c C C 0. 6 tE- 0.4 0.2 iI I I I I I I I I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 xFCu = a C)t Figure 7. T versus Curve d vessI - = a = Constan

.__-L Radicals Induct ion Ma in Radicals 3800 -Period Reaction Re i 7 Region Regions 3600 - 3400 3200 -- \ H2 2800 2600 T K 2400 2 2200 03'-g 2000 - 0. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 160at = 0 and p = 1 atm 1400 --- 0 02 1200 - 1000 0 0. 0. 2 0. 3 0. 4 0.5 0. 6 0. 7 0. 8 0. 9 1. 1 tx 104 (sec) Figure 8. Concentration and Temperature Profiles of H2-O2 Mixture at = 1. 0 andp= l atm.

Ceramic Stack Pebble Bed -Probe Station Combustion Chamber Air Inlet CHa Inlet Figure 9. Combustion Chamber Modifications.

1/4" OD Tube Water Out Water In -by c.... >I......... Gas Aspiration Flow B —-.. - ~ O 1/ 16" OD Tube 5/ 32" OD Tub Figure 10. Water Cooled-Gas Probe.

2500 2000 0 E150O 1000 Equivalence Ratio 700 t 0 1 2 3 4 5 6 7 8 Axial Distance, X, feet Figure 12. CHg4Air Axial Temperature Profile; Stoichiometric, Lean, Rich.

3500 * Results of New Setup (Ceramic Exhaust Pipe) 3230 * 3000 2900 2650 ~250 A 200 150 Without Flame Holer 000 1 2 3 4 5 6 7 i Axial Distance, X, feet. Figure 11. CH4-Air Axial Temperature Profile, k = 1. 0.

2000 X = l ft. 1500 500 0.2'0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Equivalence Ratio, q Figure 13. CH4-Air, Temperature vs. Equivalence Ratio.