UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
A NONEQUILIL7.ILC1 ANALYSIS OF THE OTTO CYCLE ENGINE FROM THE V: E',.FOINT OF COMBUSTION AND EMISSIONS By SHESH KRISHNA SRIVATSA A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1974 ACKNOWLEDGEMENTS I would like to express my thanks to Dr. V. P. Roan, the Chairman of my supervisory committee, for suggesting this study, and for his guidance and constructive criticism during the course of this work. Thanks are also due to Professors R. B. Gaither, R. K. Irey, D. Kim, and U. H. Kurzweg for serving as members of the supervisory committee. Finally, the help of Mrs. Jeanne Ojeda in typing the manuscript is gratefully acknowledged. TABLE OF CONTENTS Page ACKNOWLEDGEMENTS ................................................. ii LIST OF TABLES..................................................... vii LIST OF FIGURES...................................................viii KEY TO SYMBOLS ..................................................... x ABSTRACT...........................................................xiii CHAPTERS: I INTRODUCTION........................................... 1 II SPARK IGNITION ENGINE COMBUSTION AND RESULTING EMISSIONS............................................... 3 2.1 Introduction... ..................................... 3 2.2 Crankcase Blowby............................. ...... 3 2.3 Evaporation Losses.......................... ....... 4 2.4 Formation of Pollutants............................ 5 2.4.1 Hydrocarbons........................... ... ... 5 2.4.2 Carbon monoxide ............................. 7 2.4.3 Oxides of nitrogen.......................... 7 2.5 Effects of Various Parameters on Engine Emissions.. 8 2.5.1 Hydrocarbons ................................. 8 2.5.2 Carbon monoxide............................. 9 2.5.3 Oxides of nitrogen.......................... 10 2.6 Emission Standards................................. 10 TABLE OF CONTENTS (CONTINUED) CHAPTERS: Page 2.7 Control Methods................................. 11 2.8 Description of the Combustion Process in an Otto Cycle Engine............................... 13 2.9 Laminar Flames: Introductory Remarks........... 14 2.10 Closing Remarks................................ 16 III LITERATURE REVIEW..................................... 17 References for Chapter III............................. 33 IV MATHEMATICAL ANALYSIS OF THE COMBUSTION PROCESS........ 37 4.1 Formulation of the Problem...................... 37 4.1.1 The conservation equations................ 37 4.1.2 Major assumptions for analysis........... 40 4.1.3 Boundary conditions...................... 43 4.2 Method of Solution .............................. 50 4.2.1 Initialization........................... 53 4.2.2 Details of computation................... 53 4.2.3 Step size and stability.................. 62 4.3 Chemical Kinetics............................... 67 4.4 Uncertainty in the Definition of Temperature.... 71 References for Chapter IV............................ 76 V MATHEMATICAL SIMULATION OF THE OTTO CYCLE............. 78 5.1 Introduction.................................... 78 5.2 Present Work..................................... 81 5.3 Major Assumptions for Analysis.................. 83 TABLE OF CONTENTS (CONTINUED) CHAPTERS: V 5.4 Theoretical Analysis............................ 1. Calculations of conditions at beginning of compression................. ............... 2. The compression process...................... 3. Combustion................. .................. 4. Expansion stroke analysis.................... 5.5 The Kinetics of Nitric Oxide Formation........... 5.6 Details of Calculations for Nitric Oxide Formation ........................................ 5.7 Injection of WaterAlcohol Mixtures.............. References for Chapter V.............................. VI RESULTS AND DISCUSSIONS................................ 6.1 The Combustion Process............................ 6.1.1 Concentration and temperature profiles.... 6.1.2 Unity Lewis number approximations......... 6.1.3 Twostage hydrocarbon combustion.......... 6.1.4 Closing remarks.................... ........ 6.2 Mathematical Model for the Otto Cycle............ 6.2.1 The conibustion process.................... 6.2.2 The expansion process..................... 6.2.3 Wateralcohol injection................... References for Chapter VI............................. VII CONCLUSIONS AND RECOMMENDATIONS....................... Page 86 86 88 90 102 107 114 115 119 122 122 122 131 132 136 137 137 141 150 157 158 TABLE OF CONTENTS (CONTINUED) APPENDICES: I II III BIOGRAPHICAL EQUILIBRIUM CALCULATIONS........................... References for Appendix I.......................... THE DETERMINATION OF TRANSPORT PROPERTIES.......... References for Appendix II......................... COMPUTER PROGRAMS.................................. References for Appendix III........................ SKETCH.............................................. Page 161 172 173 197 200 334 335 LIST OF TABLES Table: Page 4.1 Reactions considered in this analysis.................... 69 5.1 Reaction mechanism....................................... 103 6.1 Results for wateralcohol induction ...................... 151 LIST OF FIGURES Figure: Page 4.1a Staggered mesh scheme........... ....................... 51 4.1b Enlarged view of staggered mesh........................ 51 5.1 pV diagram for the Otto cycle ......................... 79 5.2 Chemical reaction times for hydrogenair and hydrocarbon air mixtures at one atmosphere and unity equivalence ratio...................................... 92 5.3 Schematic diagram illustrating the sequence of processes occurring during combustion................... 94 5.4 Mass fraction of fuelair mixture burnt as a func tion of crank angle..................................... 96 5.5 Experimental and theoretical nitric oxide formation rates for combustion of hydrogenair mixtures........... 112 5.6 Experimental and theoretical nitric oxide formation rates for combustion of propaneair mixtures............ 113 6.1 Temperature profile at different instants of time illustrating the development of the steadystate profile................................................ 123 6.2 Concentration profiles for the species CEH, CO, CO2, H 2 and 02 ............................................ 124 6.3 Concentration profiles for the species H2, OH, NO, NO2, and N 20 .......................................... 125 6.4 Concentration profiles for the species H and 0......... 126 6.5 Concentration profiles for the species CH3, CHO, and N. 127 6.6 Test of unity Lewis number approximation for methane and carbon monoxide in a methaneair flame.............. 133 6.7 Test of unity Lewis number approximations for methane and carbon monoxide in 0.1 atm. methaneoxygen flame... 134 viii LIST OF FIGURES (CONTINUED) Figure: Page 6.8 Steadystate temperature profiles for unity and nonunity Lewis numbers............................... 135 6.9 Pressure and temperature histories during combus tion process............................................ 139 6.10 Nitric oxide formation during combustion................ 140 6.11 Pressuretime and temperaturetime plots for the expansion stroke....................................... 142 6.12 Concentrationtime histories of the species CO, H, O, OH, and NO during the expansion stroke.............. 144 6.13 Concentrationtime history of species N during the expansion stroke ....................................... 145 6.14 Concentrationtime history of species 02 during the expansion stroke....................................... 146 6.15 Concentrationtime histories of the species H2 and NO2 during the expansion stroke........................ 147 6.16 Concentrationtime histories of the species CO2 and H20 during the expansion stroke........................ 148 KEY TO SYMBOLS a = Velocity of sound A = Preexponential (frequency) factor in Arrhenius rate constant C = Specific heat at constant pressure D.. = Binary diffusion coefficient for species i and j 1J DT,i = Thermal diffusion coefficient for species i E = Activation energy f = Mass fraction of residual burnt gases in the engine cylinder f. = External force per unit mass on species i 1 h = Heat transfer coefficient h. = Specific enthalpy of species i i h = Standard heat of formation per unit mass of species i at 1 temperature To (= 298 degrees K) k = Boltzmann's constant k = Forward rate constant kb = Backward rate constant K = Equilibrium constant H = Mach number; number of elementary reaction steps n = Number of moles per unit volume n. = Molar concentration of species i 1 N = Number of chemical species present N. = Molar flux of species i= n.V. p = Thermodynamic pressure p = Pressure tensor > q = Heat flux vector Q = Heat transfer r = Radius R = Universal Gas Constant u s = Entropy t = Time T = Temperature v = Mass average velocity of the gas mixture V. = Diffusion velocity of species i 1 w. = Rate of production of species i by chemical reaction 1 W. = Molecular weight of species i S= ole fraction of species i X. = Masole fraction of species i 1 Y. = Mass fraction of species i U = Internal energy y = Shear viscosity A = Thermal conductivity V.i = Stoichiometric coefficient of species i in jth elementary 1J reaction step 0 = Collision integral = Intermolecular potential function p = Density 0 = Collision diameter 6 = Crank angle Subscripts: b = Burnt gases BDC = Bottom dead center f = Flame element TDC = Top dead center u = Unburnt gases Superscripts: o = Designates standard state  = Written over a symbol designates an average quantity P = Product species R = Reactant species Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A NONEQUILIBRIUM ANALYSIS OF THE OTTO CYCLE ENGINE FROM THE VIEWPOINT OF COMBUSTION AND EMISSIONS By Shesh Krishna Srivatsa August, 1974 Chairman: Dr. V. P. Roan Major Department : Mechanical Engineering An analytical study of the thermodynamic processes occurring in an Otto cycle engine has been conducted with a view to obtaining a better understanding of the combustion process and a realistic prediction of the emission levels. Much of the work previously reported in this area has been experimental, the reason being the extreme complexity of any realistic solution to the governing equa tions. The availability of highspeed computing machines provides the necessary tool for such an approach. With such machines, theoretical calculations have assumed a new importance and in many cases are being used to obtain estimates which would require extensive experimental work. Theoretical analysis, by its very nature, involves the use of simplifying assumptions to reduce any real problem to a stage where it is solvable. Any unrealistic simplifying assumption will reduce the solution to one of academic interest only and the literature on any subject is replete with such examples. In this study, an analysis with a minimum of simplifying assumptions has been carried out and the results obtained are realistic and usable. xiii The study presented here has been divided into two parts. Firstly, a study of the combustion process has been carried out by numerically solving the set of partial differential equations govern ing onedimensional unsteady laminar flame propagation in a homogeneous fuelair mixture. Secondly, a mathematical model for the Otto cycle has been developed to predict the thermodynamic properties and the species concentrations throughout the cycle based on a nonequilibrium analysis and finite rate combustion. For the first part, an explicit finite difference method has been used to numerically solve the governing equations and thereby obtain a detailed description of the flame structure in terms of the distribution of the state properties and species concentrations through the flame zone. Most of the previous analyses have ignored the varia tions in the thermodynamic and transport properties with temperature, pressure and chemical composition and have also treated the chemical kinetics as being governed by a single overall reaction. These particular simplifying assumptions (particularly the last one which is rather unrealistic for any system of practical interest) have been removed in the computational scheme presented here. The scheme is general and considers variable transport and thermodynamic properties and a realistic reaction mechanism in terms of a set of elementary reaction steps. The validity of the procedure has been demonstrated by carrying out calculations for the methaneair flame and obtaining good agreement with the limited experimental data available. The mathematical model of the Otto cycle gives a detailed description of the cycle in terms of the thermodynamic properties and species concentrations at any stage in the cycle, including the xiv exhaust. The combustion process has been considered to proceed at a finite rate instead of the commonly made assumption that this process is instantaneous. Variations in the thermodynamic properties with temperature and chemical composition have been allowed instead of assigning them constant values. An accurate prediction of the concentrations of the various species at exhaust was the main objective in the development of the model. The model has been used to study the effect of the intake of an additive mixture of water and alcohols (methanol and ethanol) in different proportions and amounts on the emission levels. The model is quite flexible and permits the study of the influence of other factors such as exhaust gas recirculation, airfuel ratio, compression ratio, speed, ignition timing, etc., on the performance and emission levels. CHAPTER I INTRODUCTION The role of the internal combustion engine as a major contributor of the hydrocarbon gases, carbon monoxide, and oxides of nitrogen that contaminate the atmosphere has been recognized for some years. The formation of these pollutants is closely related to the combustion process and a realistic modeling of this process would help the study of the methods to reduce the emissions. A knowledge of the transient mechanism of high temperature oxidation of hydrocarbon fuels would lead to a better understanding of the combustion process and may ultimately help in the design of engines which will emit lower concentrations of pollutants. Nhile the amount of unburnt or partially burnt products in the exhaust is generally insignificant from an energy point of view, it is significant, and very much so, when considered from an environmentalist's point of view. Thus these studies leading to new designs might or might not improve the overall efficiency of the internal combustion engine (in fact a reduction in efficiency might be achieved), but they might help to reduce the emissions. The high temperature oxidation mechanism of some simple fuels like hydrogen and carbon monoxide have been widely studied and are fairly well understood, but the state of development for the oxidation of hydrocarbon fuels leaves much to be desired. Whatever little work has been reported is for lower hydrocarbons (like methane, etc.) and kinetic data and reaction mechanism postulations for some practical fuels (like isooctane) are virtually nonexistent. One of the methods available for the quantitative determination of chemical kinetic constants is a detailed analysis of the flame structure. A comparison of the flame structure predicted theoretically on the basis of a postulated reaction mechanism, with the experimentally measured data, provides a convenient way of checking the validity of the reaction mechanism. The influence of varying the rate constant parameters on the theoretical predictions permits the determination of influence coefficients for the different input parameters. A parametric study of this type is useful in arriving at some sort of a conclusion without much trial and error. To put this principle into practice, there remains the need for a general theoretical computational procedure, and it is the aim of this study to proceed in this direction by presenting a fairly general computational scheme for a realistic reaction mechanism. In addition to the study of the combustion process, a mathematical model for the entire Otto cycle has also been presented. The model can be used to study the effect of different parameters on the exhaust emissions. CHAPTER II SPARKIGNITION ENGINE COMBUSTION AND RESULTING EMISSIONS 2.1 Introduction In this chapter, the formation of the various pollutants, the factors affecting them, and the various control methods available are briefly discussed. Finally, a description of the combustion process together with some introductory remarks on laminar flames is provided. There are three important sources of emissions from an auto mobile: the engine exhaust, the crankcase blowby, and evaporation losses from the fuel system vents, the engine exhaust emitting by far the largest amount of pollutants. Thus most of the research on automobile pollution has been concerned with reducing the exhaust emissions. The role of the other two sources of emissions is discussed briefly for the sake of completeness. 2.2 Crankcase Blowby This refers to the gases escaping from the engine cylinders into the crankcase through the gaps between the sealing surfaces of the piston and the cylinder walls. This leakage which occurs mainly during the expansion and the compression strokes increases with increasing load (increased pressures and mass flow) and typically consists of about 85% unburnt fuelair mixture and 15% combustion products. In terms of hydrocarbon concentrations, this is about 600012000 parts per million. On uncontrolled engines the crankcase emissions account for about 20% of the total emissions from an automobile. Crankcase emission control systems are at present in a well developed state, so that this source of emissions has been virtually eliminated. All the control systems use the same basic approach: recycling of the gases from the crankcase back into the combustion chamber, and this is commonly referred to as 'positive crankcase ventilation' (PCV). 2.3 Evaporation Losses These losses occur from the carburetor and the fuel tank, and on uncontrolled engines account for about 1820% of the total hydro carbon emissions. Losses from the carburetor occur mainly after the engine is shut off when the temperature of the fuel system rises due to the heat stored in the engine. The losses depend upon the operating conditions immediately before engine shutdown, the fuel volatility, carburetor bowl capacity, and the ambient temperature. The losses from the carburetor are generally small when the engine is running, and the fuel tank accounts for the majority of the evaporation losses. The losses from the fuel tank are dependent on the fuel volatility, its temperature, what fraction of the tank is filled with fuel, and the ambient temperature. Fuel evaporation loss control devices store the vapors from the fuel tank and the carburetor. The vapors are subsequently burnt in the engine combustion chamber. The various controls in use are based on different modifications of this basic principle. Evaporation loss controls have been in effect since 1972. As mentioned before, these two losses (crankcase and evaporation) contribute about 20% each to the total emissions of an uncontrolled vehicle. The controls developed for minimizing the pollutants emitted thereby are in a fairly good state of development so that these two sources of emission have been virtually eliminated. All modern auto mobiles (1972 and later year models) are equipped with these controls so that presently the engine exhaust accounts for almost all the pollutants. 2.4 Formation of Pollutants We shall now look at the formation of the various pollutants which occur in the exhaustour primary concern here. 2.4.1 Hydrocarbons The presence of hydrocarbons in the engine exhaust shows that combustion failed to occur in certain parts of the engine cylinder. Calculations (chemical equilibrium or finite rate) reveal that the concentrations of the hydrocarbons in high temperature combustion products are negligibly small. The source of the unburnt hydro carbons is not this fairly homogeneous part of the combustion chamber. The major source of the unburnt hydrocarbons is the quench layer adjacent to the relatively cool combustion chamber walls, wherein the oxidation of the fuel is only partially completed. Further, due to improper mixing, the local mixture composition may not be within the flammability limits and this would partly account for the unburnt hydrocarbons. In an Otto cycle engine the quenching is by far the more important of the two. The quench layer may be considered to consist of two zones: one known as the total quench zone which is immediately adjacent to the wall and where virtually no chemical reactions occur, and the other known as the partial quench zone where the temperature is high enough for cool flame and slow oxidation reactions to take place, but not at a fast enough rate to cause complete combustion. After flame propagation, some of the hydrocarbons in the quench layer mix with the burnt gases and experience partial oxidation. The availability of oxygen and the degree of turbulence strongly influence the extent of this oxidation. Also all the hydrocarbons present in the quench layer are not expelled during exhaust and some remain as residuals in the cylinder for the next cycle. Reactions continue in the gases expelled into the exhaust pipe causing further oxidation of the unburnt hydrocarbons. For these reasons, the amount of hydrocarbons actually present in the exhaust gases is considerably less than the hydrocarbons present in the initial quench layer. Besides the quench layer the crevices between the piston rings and other dead spaces and pockets in the combustion chamber which prevent flame propagation also contribute to the exhaust hydrocarbons. The hydrocarbons in the exhaust are a complex mixture of oxygenated and partially oxidized compounds together with products produced by the cracking of the fuel; thus the exhaust hydrocarbons are different in composition from those produced by crankcase blowby and evaporation losses which contain mainly the fuel itself. 2.4.2 Carbon monoxide This is a product of incomplete combustion and results when sufficient oxygen is not available for combustion, either on an overall or a local basis. Significant amounts of carbon monoxide may be present in the combustion chamber immediately after the combustion is complete even when sufficient oxygen is available. This can be shown by a chemical equilibrium analysis. However, as the gases cool during the expansion stroke, much of this carbon monoxide gets converted to carbon dioxide, so that except for fuelrich mixtures, the amount of carbon monoxide in the exhaust may be of the order of about 1% by volume or about 10,000 parts per million. The conversion of carbon monoxide to carbon dioxide during the expansion stroke cannot be predicted from chemical equilibrium considerations since the chemical reactions involved are not rapid enough to reach equilibrium. As such, the actual concentration of carbon monoxide at exhaust is higher than that if chemical equilibrium were to prevail. 2.4.3 Oxides of nitrogen A brief discussion is provided here; details may be found in Chapter V. The most important of the oxides of nitrogen of concern here is nitric oxide (NO). The formation of nitric oxide, a high enthalpy species, occurs due to the existence of high temperatures during the combustion process. The process of formation is slow relative to the time scale of other processes occurring in a typical internal combustion engine, and so has to be considered on a rate dependent basis and not on a chemical equilibrium basis. During the expansion process, the concentration of nitric oxide is essentially frozen at an early stage and hence the concentration at exhaust is higher than that if chemical equilibrium were to prevail. 2.5 Effects of Various Parameters on Engine Emissions 2.5.1 Hydrocarbons The hydrocarbon emissions depend primarily on the airfuel ratio and engine load conditions. The emissions decrease with an increase in the airfuel ratio to a certain extent until the mixture is so lean that regular flame propagation is hindered. This is commonly referred to as the 'misfire' condition. With further increase of the airfuel ratio, the hydrocarbon emissions increase. By increasing the turbulence and having heated inlet manifolds, regular flame propagation can be extended to very lean mixtures. In terms of the operating condi tions the influence of the airfuel ratio may be summed up as: hydrocarbon emissions are a maximum during engine warming, idle and deceleration. However, since the amount of air intake under these conditions is small, the absolute amount (gm./hr.) of emissions is not necessarily a maximum. The emissions also increase with increasing load and compression ratio and early spark timing. With greater spark retard, combustion is completed later during the expansion stroke and is spread over a longer time. The exhaust temperatures are also higher, thus promoting oxidation in the exhaust system. The result is lower emissions, but the performance may be adversely affected. The emissions decrease with speed to a certain extent. At very high speeds, the ignition voltage available at the spark plugs decreases due to the short time available for coil buildup. This causes only partial combustion, increasing the concentrations of unburnt hydro carbons. The recently introduced capacitivedischarge ignition systems maintain a high enough voltage irrespective of the operating speed. Increasing the average operating temperature and freeing the combustion chamber from deposits serve to decrease the hydrocarbon emissions. The quenching effect and hence the emission of hydrocarbons are decreased by a decrease in the surface area to volume ratio within the combustion chamber, so that a larger part of the mixture is in the flammable zone. This is achieved by increasing the compression ratio, increasing the stroke to bore ratio and modifying the geometry to obtain compact combustion chambers avoiding dead spaces and pockets. Increasing the combustion chamber surface temperature also serves to decrease the extent of the quenching effects. 2.5.2 Carbon monoxide Like hydrocarbons, the carbon monoxide emissions also depend on the airfuel ratio and engine load conditions. The carbon monoxide emissions are a maximum during idling and deceleration when the air fuel ratio is towards the fuelrich side. As with hydrocarbons, the absolute amount (gm./hr.) of the emissions under these conditions need not necessarily be at its maximum value. The carbon monoxide emissions also increase with load and during acceleration periods. The concentration is generally at a minimum during moderate power cruising. Other factors which affect the hydrocarbon emissions also affect the carbon monoxide emissions in a similar manner. 2.5.3 Oxides of nitrogen A detailed discussion is provided in Chapter V; here only brief mention is made of the effect of the important parameters. The nitric oxide emissions depend primarily on three factors: (i) The emissions are a maximum at an airfuel ratio corresponding to about 10% excess air (as compared to stoichiometric), and decrease for richer or leaner mixtures. (ii) Higher combustion temperatures favor larger concentrations of nitric oxide. (iii) The emissions are greater, the greater the time available for the formation of nitric oxide, this being a strongly timedependent process. 2.6 Emission Standards The standards set of 1976 vehicles are (based on constant volume sample cold start test average with a constant volume sample hot start test, both with the Federal 22minute driving cycle): 0.41 gm./vehiclemile for hydrocarbons, 3.40 gm./vehiclemile for carbon monoxide, and 3.00 gm./vehiclemile for oxides of nitrogen. Oxides of nitrogen have to be cut to 0.4 gm./vehiclemile for 1977 vehicles. For the sake of comparison, the emission levels without any exhaust controls for a typical automobile are given below: Hydrocarbons: 850 ppm (approximately 11 gm./vehiclemile) Carbon Monoxide: 3.4% by volume (approximately 80 gm./vehiclemile) Oxides of Nitrogen: 1000 ppm (approximately 4 gm./vehiclemile) With control techniques in current use the emission levels are approximately (1973 standards based on constant volume sample cold start test): 3.4 gm./vehiclemile for hydrocarbons, 39.0 gm./vehiclemile for carbon monoxide, and 3.0 gm./vehiclemile for oxides of nitrogen. [These figures would be somewhat different (lower for hydrocarbons and car bon monoxide) with the testing procedure employed for the 1976 standards.] 2.7 Control Methods A brief mention of the various methods used to control the exhaust emissions is made here. (1) Leaner airfuel ratios tend to decrease both the hydrocarbons and carbon monoxide emissions; oxides of nitrogen are also reduced after a certain airfuel ratio. However, at these airfuel ratios if regular flame propagation cannot occur, the former two emissions may actually increase rather than decrease. Presently, research is being conducted to study the extension of the lean limit and also how to cope with the generally adverse effects on driveability and some other operational characteristics of the engine at extremely lean airfuel ratios. (2) Secondary air pumped either into the engine cylinders directly during the expansion stroke or later in the exhaust manifold brings about a reduction in the hydrocarbon and carbon monoxide emissions by promoting more nearly complete oxidation. The nitric oxide emissions are not significantly affected, they may increase slightly. (3) Design modifications: These include improved carburetor performance, a fast acting and more accurate choke, an inductive or electronic igni tion system with modified timing, reduced compression ratios, air preheating, etc. Some of these call for a compromise between the operational characteristics and the emission levels. Previously, the best design was considered to be the one that resulted in the best operational characteristics, and this is generally different from the condition for lowest emission levels. (4) Postcylinder exhaust treatment of the emissions prior to their exit from the exhaust system can be done in several ways besides secondary air injection mentioned above. The most important of these is the use of catalytic converters to promote the reactions in the exhaust system in a direction which favors the oxidation of the unburnt or partially burnt products, and also to decompose nitric oxide into molecular nitrogen and oxygen, the two processes generally being carried out in different stages. Extensive research has been done in this area but still some problems remain before the system can be put in operation on an extensive commercial basis. The high initial cost involved, the interference of the lead in the gasoline with the catalyst action, and the associated reduction of effective lifetime are some of the problems to be considered. Other modes of exhaust treatment include flame afterburners, but these do not result in the reduction of nitric oxide emissions, and besides require considerable auxiliary equipment. (5) Exhaust gas recirculation serves to reduce the nitric oxide emissions, but the carbon monoxide and hydrocarbon emissions are adversely affected as are some of the operational characteristics of the engine. There are several other control techniques which have been suggested and the reader is referred to the appropriate literature on the subject. 2.8 Description of the Combustion Process in an Otto Cycle Engine In a spark ignition engine, the high temperature of the electric spark causes the ignition of a small portion of the fuelair mixture in the immediate vicinity of the spark plug gap. The result is a more or less spherical volume of hot gas consisting of burnt products which are rapidly equilibrated due to the reactions being very fast at the high temperature involved. This causes a nearly spherical flame to originate from the point of ignition. Spherical flames are truly one dimensional systems by virtue of their geometry, whereas flat flames are quasionedimensional systems, there being small but negligible changes in the various state properties in the direction perpendicular to that of propagation. The fuelair mixture in the combustion chamber is exposed to the flame, layer by layer, until combustion is complete as far as practicable. In practice, combustion is never totally complete due to the quenching effects of the walls and other reasons. The entire combustion process is extremely rapid, e.g., for combustion spread over 40 degrees of crank revolution at 2000 revolutions per minute, the process would be completed in approximately 3.33 milliseconds. In spite of the rapidity of the process, at any given instant of time, the fraction of the total fuelair mixture actually being burnt is not significant if normal combustion is occurring. This would not be the case when detonation and knocking occur. Therefore, under normal circumstances, the combustion process occurs with uniform release of energy and may be approximated in the initial stages by the model of a onedimensional laminar spherical flame propagating into a premixed fuelair mixture. The assumption of a laminar flame is made to render the problem solvable (in general, turbulence would be present in varying degrees), there being no convenient numerical techniques to handle a turbulent flame for a general case. As the flame propagates, it would deviate from its spherical shape, the deviation depending on the geometry of the combustion chamber walls with respect to the point of ignition. This deviation can be accounted for by considering the general three dimensional flame equations in their timedependent form, and a solution so obtained would give the flame shape and position as a function of time. Even with modern electronic computers this problem is almost insurmountable, and would hardly give more insight to the actual combustion process than the solution of a onedimensional problem. 2.9 Laminar Flames: Introductory Remarks Laminar premixed flames represent a unique case of wave propagationthe propagation of a selfsustaining exothermic reaction wave into a mixture of combustible gases. Heat conduction, diffusion, and a set of rapid exothermic chemical reactions all are significant in the propagation of such a flame. Thus this is one of the earliest combustion problems that has necessitated the simultaneous consideration of both the fluid dynamic and chemical kinetics aspects for its solution. Laminar flame velocity (also referred to variously as the burning velocity or the normal combustion velocity), which is of primary importance in the study of premixed flames, is defined as the velocity at which unburnt gases move through a combustion wave in the direction normal to the wave surface. This definition refers to the velocity of unburnt gas relative to the flame front as determined at some point far enough removed from the flame, such that the influence of the flame itself at the point is negligible. The propagation characteristics of steady, adiabatic, one dimensional laminar premixed flames are governed primarily by two properties of the systemthe activation energy of the rate controlling step and the Lewis number, where the Lewis number represents the ratio of the energy transported by conduction to that transported by diffusion. The geometry is of secondary importance, and therefore, within reason able limits, the laminar flame velocity may be defined independent of a particular experimental apparatus. It will be seen later that the flame velocity is an eigenvalue of the conserv ,. equations. A laminar flame is made up of a series of fairly welldefined regions. The thickness of the flame zone is rather arbitrary, similar to the thickness of a boundary layer. In the first region, the unburnt gas undergoes appreciable changes in composition and temperature due to the transport processes of diffusion and thermal conduction, but the reaction rates are negligible. This is referred to as the preheat zone. After this there is a relatively thin primary reaction zone wherein all the important rapid chemical reactions occur. This is the luminous zone in hydrocarbon flames. In addition to this, there may be one or more secondary reaction zones whose spatial dimensions are an order of magnitude larger than those of the primary reaction zone, and wherein slow reactions (e.g., carbon monoxide afterburning in hydrocarbon flames) and radical recombination reactions occur. At low enough pressures the radical recombination reactions and the afterburning may occur in distinctly separate zones; on the other hand, with increasing pressure, the secondary reaction zone may partly merge into the primary zone. The existence of these zones and carbon monoxide afterburning has been confirmed by experiments. 2.10 Closing Remarks An introduction to the study of the combustion process in and the emissions from an Otto cycle engine has been presented. These ideas will be developed further in the succeeding chapters. Chapter III presents a review of the extensive literature on the methods of solution of the laminar flame propagation equations. Chapter IV presents a mathematical analysis of the combustion processthe formulation and solution of the laminar flame propagation equations. In Chapter V a mathematical model has been developed to obtain a detailed thermodynamic analysis of the Otto cycle. Results and discussions are presented in Chapter VI and finally Chapter VII is used for some closing remarks. CHAPTER III LITERATURE REVIEW The problem of laminar flame propagation in premixed gases has been studied fairly extensively. The problem was first studied by Mallard and LeChatelier (1883). Their theory has been referred to as the thermal theory in the sense that they considered thermal effects to be of primary importance and the chemical kinetics to be secondary. Later on, Taffanel (1913) and Daniell (1930) demonstrated for a highly simplified model that the flame propagation velocity is proportional to the square root of the ratio of the thermal conductivity to the specific heat at constant pressure and to the square root of the rate of the reaction. An excellent review of the early literature (up to about 1950) on the steady state flame propagation in premixed gases has been presented by Evans [1]. Williams [2, p. 95] has presented a review with references to more recent literature (up to about 1960). The assumption of Lewis number equal to unity (normal diffusion) has been widely made in flame analysis. This assumption implies that the energy fluxes due to diffusion and conduction balance each other, and so the enthalpy is constant throughout the flame zone; hence the enthalpy equation need not be integrated. This simplifies the problem considerably and even permits closed form solutions for a few simple cases. The assumption of unity Lewis number is, however, not valid for many gases encountered in combustion problems, two notable examples being hydrogenair and propaneair mixtures (Le = 3.30, and Le = 0.49, respectively). Assuming that the activation energy of the rate controlling step and the heat of reaction are large and that all the reaction occurs at the hot boundary, Zeldovich and FrankKamenetski [3] were able to obtain a closed form solution for the burning rate eigenvalue. Semenov [Sokolik, 4] has modified their solution for nonnormal diffusion by multiplying the flame velocity as calculated by the ZeldovichFrank Kamenetski formulation by the factor 1/vLe for first order reactions and by 1/Le for second order reactions. This analysis is also limited by the requirement of large or infinite activation energies. Von Karman's zeroth and first order approximations [5] also yield closed form solu tions, the former being a special case of the ZeldovichFrankKamenetski analysis wherein only the first term of a series expansion for the eigenvalue is retained. Some of the other analytical solutions are those of Boys and Corner [6, 7] and Adams [8]. Wilde [9] empirically modified Adams' solution thereby obtaining better agreement with 'exact' numerical solutions. Hirschfelder [10] has obtained two approximate solutions in one of which the eigenvalue is obtained as a root of a quadratic equation and in the other as a root of a third order equation. The accuracy is comparable to that obtained with Wilde's solution. Williams [2] has given a graphical comparison of the accuracy of these approximate methods. Spalding [11] has proposed a centroidd rule' which provides a convenient way of estimating the flame speed to an accuracy of about 1% for the case in which the reaction rate is an explicit function of temperature. This is true for unity Lewis number of the reactants and radicals for the following class of reactions: (a) single step reactions, (b) reactions satisfying the steady state approximation, and (c) explosive chain reactions. The method gives an expression for the flame speed in terms of the centroid of the reaction rate function without resorting to the solution of the governing dif ferential equations. Graphs are provided for estimating the flame speed for nonunity Lewis number. These show that when the Lewis number of the radicals is greater than unity the flame speed in creases, but only slightly. Small values of Lewis number (of the radicals), on the other hand, greatly reduce the flame speed. For the reactants, a Lewis number greater than unity decreases the flame speed. Rosen [12, 13] has suggested a variational method employing the RayleighRitz method for obtaining successively better estimates of the functional and eigenvalue. He has also given a plausible derivation of the centroid rule proposed by Spalding [11]. The centruid rule method has been extended for the case of nonnormal diffLusion with bimolecular reactions (A+BC) in a stoichiometric twocomponent system by Adler [14]. The flame speed has been shown to be inversely proportional to the Lewis number of the reactants for large activation energies. Spalding [15] has compared several approximate solutions for the case of unity Lewis number and temperature explicit reactions and arranged them in order of merit as: Wilde (best), Adams, Von Karman Penner, ZeldovichFrankKamenetski, and Corner. Wilde's method has been extended to study the effect of the diffusion coefficient of reactants and radicals on the flame speed and fairly good agreement with 'exact' solutions is reported. Adler [16] has considered temperature explicit reaction rate functions and obtained limits on the burning rate eigenvalue independent of the type of reaction considered for a fixed centroid of the reaction rate function. For the limiting case of high activation energy the upper and lower limits are close, so that the flame speed may be obtained with good accuracy. Favin and Westenberg [17] have considered onedimensional steady premixed spherical laminar flames and solved the equations for the case of unity Lewis number and a unimolecular reaction process. They have integrated the equations from the cold to the hot boundary and considered the flameholder temperature, rather than the mass flow rate, as the eigenvalue. Adler and Spalding [18] have obtained analytical solutions for the temperature gradient, radical concentrations, and the flame speed for the case of nonbranching chain reactions with negligible chain breaking. De Sendagorta [19] has obtained a numerically exact solution for one value of the reduced activation temperature E/RTb and used this to show that good accuracy could be obtained by using an exponential approximation for the species flux fraction in analytical solutions. His solution covers both first and second order single step reactions and is considerably more accurate than the other analytical solutions cited above. Nonnormal diffusion has been considered by Wilde [9], de Sendagorta [19] and Von KarmanPenner [20] utilizing KarmanPohlhausen profile methods. Millan and da Riva [21] have obtained a solution for the distribution of radicals in laminar flames based on the steady state approximation. A general criterion for the applicability of this approximation is given and it is shown that this assumption is a first order approximation to the actual solution, and this solution can be improved by a method of successive approximations. Spalding and Adler [22, 23] have considered onedimensional laminar flames with an enthalpy gradient (defined as positive for heat loss at the cold boundary) which is useful when heat losses by conduction or radiation are significant, and obtained numerical solutions. For unity Lewis number, they have shown that the flame speed decreases or increases as the enthalpy gradient decreases or increases from zero, the effect being more, the higher the activation energy. The flame speed decreases to zero for a certain critical negative value of the enthalpy gradient. Heat losses and flame propagation limits due to heat losses have also been considered in the analyses presented by Mayer [24], Spalding [25], Adler [26], and Idler and Kennerley [27]. From these analyses it can be concluded that for a given amount of heat loss, a given combustible mixture has two possible flame speeds, the lower one of which is normally unstable. With increasing heat loss, the two speeds approach each other and finally coincide for a critical value of the heat loss. For still greater heat loss, the flame speed is imaginary. The critical value of heat loss is identified with the flammability limit. Some other studies of nonadiabatic flames are given in References 28, 29, and 30. With the advent of modern electronic computing machines, the emphasis has gradually shifted to obtain more meaningful solutions with less drastic assumptions, by using numerical methods. These studies involve less stringent assumptions, are more accurate, and include better estimates of thermodynamic and transport properties, and also consider the effects of chain reactions instead of the classical unimolecular process R P (R = reactant, P = product), considered in most of the previous analytical solutions. Firstly, brief mention is made of the early numerical solutions. Von Karman's second approximation [5] involves an iterative procedure with numerical integration. The approximate solution of Boys and Corner [6] can be improved by numerical iterations. Some of the other iterative procedures requiring numerical integration are those due to Klein [31] and Johnson and Nachbar [32]. The latter one is the most accurate of the simple approximate procedures and it permits the determination of the upper and lower bounds for the eigenvalue. Using their iteration scheme successive upper and lower bounds can be obtained until the eigenvalue converges to any desired accuracy. Generally the first approximation is accurate enough for practical purposes so that it is seldom necessary to use the iterative scheme. Friedman and Burke [33] have obtained numerical solutions for the case of a hypothetical first order irreversible reaction with exponential dependence of rate on temperature over a wide range of two parameters: dimensionless activation temperature (E/RTb) and the Lewis number. The theoretical analysis presented by Giddings and Hirschfelder [34] was the first one to consider branched chain reactions. The assumption of a steady state condition was made to approximate the radical concentrations and the governing equations were solved by the integral equation approach devised by Klein [31]. Menkes [35] has solved the problem for a cylindrical flame by transforming the differential equation into an integral equation which is solved numerically by a method of successive approximations. The cylindrical flame was investigated to get an insight into the stability analysis. It was found that the flame radius decreases and the specific mass flow increases rapidly with increasing inlet temperature and at constant flame thickness. However, if the flame radius is decreased, keeping the inlet temperature constant, the flame thickness increases rapidly resulting in a decrease of the specific mass flow rate. Similar results have been obtained by Menkes [36] by an approximate analytical solution. Jain and Ebenezer [37] have obtained analytical solutions for a cylindrical flame, wherein the reaction rate is represented by a step function, and have studied the accuracy of the thinflame" approximation. A thin flame is one whose thickness is small compared to its radius of curvature, so that the mass rate of burning per unit area can be taken as equal to that in a plane flame. Spalding [38] has obtained numerical solutions for the case of a steady laminar spherical flame by stepwise integration from the cold boundary. The model used consists of a combustible gas emerging from a point source into an atmosphere of adiabatic combustion products. The results have been compared with those obtained from the thinflame theory and the two are shown to be asymptotically equal for large mass flow rates and high reaction rates. For the same model Spalding and Jain [39] have obtained analytical solutions for the cases of (a) step function reaction rate curves and (b) Adamstype reaction rate curves based on the assumption of a constant temperature gradient in the reaction zone. It is shown that the radius of the flame depends more on the position of the centroid of the reaction rate curve rather than on its specific shape. The same results have been confirmed by means of a resistance network analog computer by Spalding, Jain and Samain [40]. Jain and Kumar [41] have employed the method of matched asymptotic expansions and obtained analytical solutions for the limiting case of large activation temperatures. They have also obtained numerical solutions for a wide range of values of both the reduced activation temperature and Lewis number, for both first and second order reactions. Based on the numerical solutions, a reason ably accurate rule for the flame speed eigenvalue in terms of Lewis number and the centroid of the reaction rate function is suggested; it is linear in Lewis number for first order reactions and quadratic in Lewis number for second order reactions. There is good agreement of these results with those of Spalding [11, 15] and de Sendagorta [19]. However, the present solution agrees with Adler's [14] solution only for Lewis numbers around unity, and there is considerable disagreement when the Lewis number differs significantly from unity for which a possible explanation has been given. Besides the various methods mentioned above, other recent methods have sought not just to determine the laminar flame speed but also the detailed flame structure regarding the distribution of temperature and the various chemical species over the thickness of the flame zone. A detailed analysis of the flame structure gives a better insight into the mechanism of flame propagation and also permits the study of chemical reactions from the viewpoint of mechanism and kinetics. The difficulties involved in such a detailed analysis of the flame zone may be summed up into three categories: (a) a reliable postulation of the reaction mechanism for the oxidation of the fuel; (b) a reliable estimate of the transport properties, especially at temperatures encountered in flames and for free radicals; and (c) the absence of suitable numerical techniques to handle the complicated system of equations in a convenient way. Basically two mathematical formulations have been widely used in the detailed analysis of the flame zone. The first one which is a trial and error method is due to Hirschfelder [42]. This method consists of setting up the differential equations for the conservation of mass and energy and the diffusion equations in a coordinate system which is at rest with respect to the flame front. The boundary conditions are specified at the burnt and unburnt ends. In practice, integrating these equations to obtain the concentrations of species and state properties in the flame zone requires a knowledge of the flame speed. Since this is not known a priori, the problem reduces to an eigenvalue problem; a solution to the equations exists for only one unique value of the flame speed, known as the eigenvalue of the conservation equations or the eigenvalue of the twopoint boundary value problem. Assuming various flame speeds, the equations are numerically integrated from the burnt to the unburnt end. The solution which agrees best with the boundary conditions is the one assumed to be correct. Convergence criteria for this method are not available. For more complicated cases, where the chemical kinetics is governed by a set of reaction steps, the assumption of flame speed alone is not sufficient to start the integra tion of the conservation equations, and one is confronted with a multi eigenvalue problem in nonseparable equations [1, 43, 44]. The number of eigenvalues is equal to the number of linearly independent species flux fractions in the flame equations. This multieigenvalue problem has, however, been reduced to the determination of a single eigenvalue, without loss of generality, in a method introduced by Campbell, Heinen and Schalit [44]. They have obtained numerical solutions for a one dimensional flame with idealized kinetics and transport properties. Three examples using this trial and error method have been given by Hirschfelder [42]: (a) the unimolecular decomposition of hydrazine, (b) the bimolecular decomposition of nitric oxide, and (c) the two step chain mechanism describing the decomposition of ozone. Agreement with experimental results is fair* for the first and third cases; no It should be noted that the solution method itself, even if extremely powerful, would not give results in good agreement with the experimental data if the accuracy of the reaction kinetics and transport property data is poor. Thus the judgment of accuracy should rather be based on the solution of the same problem with the same input data by different methods and comparison with experiments. experimental data were available for comparison in the second case. The second method is due to Spalding [45] and consists of setting up the differential equations for conservation of mass and energy and the diffusion equations in their timedependent or unsteady form. Assuming arbitrary initial profiles for the state properties and the species concentrations, the equations are solved by forward marching finite difference methods until the various profiles of the species concentrations and the state properties converge to a steady state. A big advantage of this method is that almost any arbitrary initial profile can be used for starting the solution; however, by using profiles which closely resemble their steady state shape, the computation can be cut down considerably. For the temperature and the major species, Sshaped profiles asymptotically approaching the boundary values are commonly used. The radicals are usually assumed to be at zero concentration initially. The assumed profiles will develop into the steady case if enough energy is stored in them. If this initial energy is insufficient then the gradients of temper ature and species concentrations will decrease monotonically and the flame is extinguished. Thus this method can also be used to obtain an estimation of the minimum ignition energy, if required. No iterations are required with this approach. The computational advantage or disadvantage of either method has not been conclusively established as yet. This method has been demonstrated by Spalding for the case of the hydrazine decomposition flame. The conservation equations reduced to two partial differential equations and these were solved simultaneously by Schmidt's graphical method for the solution of the unsteady heat conduction equation. Calculations were carried out for two different flame temperatures and fairly good agreement between theoretical and experimental values were obtained. Spalding's method has been used by various researchers for solving different flames. The later techniques used are basically refinements of the finite difference solutions and the inclusion of chainbranching reactions and other improvements. A review of some of the recent work using this method is presented below. Zeldovich and Barenblatt [46] have reviewed previous work on the steady state propagation of flat flames. They have also considered the effects of the Lewis number, heat loss and chain reactions on the velocity of propagation of flames using Spalding's timedependent approach. Using a perturbation technique the stability of the steady state solution has been examined. They have shown that for nonunity Lewis number, the total enthalpy reaches a maximum or a minimum in the flame zone, depending on whether the Lewis number is greater or smaller than unity. The same result was also obtained earlier by Lewis and Von Elbe [47]. Heat loss was found to decrease the velocity of propagation and there was found to exist a certain critical value of the heat loss parameter 'b' [q = b(tt )], above which a steady state solution did not exist. Adams and Cook [48] have used the timedependent method to study the hydrazine decomposition flame which has been fairly extensively studied. As before, the problem was reduced to the simultaneous solution of two partial differential equations: one involving the temperature and the other the mass fraction of hydrazine. A central difference formula was used for approximating the space derivatives and a forward difference formula for the time derivatives. The equations were then numerically solved by explicit forward marching steps in time until the profiles of temperature and mass fraction reached a steady state. Considering chain breaking reactions of both second and third order, it was shown that the dependence of flame speed on temper ature is insensitive to changes in pressure and the mechanism of the chain breaking reactions, while the dependence on pressure is sensitive to the mechanism assumed. DixonLewis [49] has implemented Spalding's timedependent method for the hydrogenoxygen flame. To reduce the computational work, some simplifications regarding the transport and thermodynamic properties were made. The chemical kinetics was treated in terms of a set of elementary reaction steps and thus this analysis avoids some of the misleading results of previous analyses which considered the chemical kinetics to be governed by a single reaction or a simple unrealistic mechanism. Important conclusions regarding the influence of some parameters on the hydrogenoxygen flame have been made in this analysis. For the sake of computational stability, the standard 2 1 stability condition of At/(A)2 < 1 for the heat conduction equation had to be modified to At/(AI) 2 1/30 for this problem. Detailed transport property calculations have been made [50] and considered in the analysis of the hydrogenoxygennitrogen flame by DixonLewis [51] using the RungeKuttallerson procedure (RungeKutta method with variable step size). Spalding et al. [52] have utilized the mathematical similarity of the equations governing steady twodimensional boundary layer flows with the onedimensional unsteady flame equations. The computational procedure developed by Patankar and Spalding [53] for the study of the hydrodynamical and heat transfer aspects of steady twodimensional boundary layer flows has been adapted to the study of the onedimensional unsteady hydrazine decomposition flame. Both these phenomena are governed by sets of simultaneous parabolic partial differential equations. The assumptions of unity Lewis number and constant and uniform specific heats were made in the solution method, although neither was necessary for the computational procedure. Initial concentration profiles of the form yA = Au (10w 15w + 6w ) were used for the major species, and the radicals were assumed to be initially at zero concentration. Here, w = (' )/(u b), where k = stream function, and the subscripts b and u refer to the burnt and unburnt ends, respectively. A variable step size for the space dimension was used and the equations were solved by an implicit method after linearization at each step. In the absence of formal stability criteria, the largest time step size was found by trial and error. Wilde [54] has treated this problem as a boundary value problem and formulated it in three different ways: (a) a set of first order timeindependent ordinary differential equations for species conserva tion, diffusion and thermal flux; (b) a set of second order time independent ordinary differential equations; and (c) a set of time dependent parabolic partial differential equations. The first for mulation was solved by quasilinearization, also known as Newton RaphsonKantorovich linearization, and then integrating by using algorithms suitable for stiff equations. The second formulation was also linearized and then the equations were written in finite difference form, the result being a set of linear algebraic equation, block tridiagonal in form. For the last formulation, explicit finite difference methods were found to require a very small time step size due to the presence of the term involving the rate of production of species by chemical reactions. Considerable time saving was obtained by using forward difference formulas in time, implicit representation of the space derivatives and an averaging of the chemical rate term at the present and future time step. The second formulation was faster per iteration than the first one, but required closer initial estimates and a greater number of iterations for convergence. The last method was longer in running time than the other two but was surer in convergence. Three flames were studied using this method: ozone decomposition, hydrogenbromine flame, and the hydrogenoxygen flame. The calculated flame velocities were about twice the cor responding experimental values. The discrepancy was explained by blaming the unreliable kinetic and transport data. Bledjian [55] has provided a general formulation which allows for the consideration of spherical and cylindrical flames in addition to plane flames. An important difference between plane and nonplanar flames should be mentioned here. For nonplanar flames the spatial variable (radius) appears explicitly in the gradient terms while for plane flames the spatial variables appear only in the derivatives. Thus for the simple cases, the space variable can be eliminated and the governing equations can be reduced to a single equation with temperature as the independent variable and the mass flux as the dependent variable. An added advantage of this is that the range of integration is now finite in the independent variable. The method of lines was used to solve the governing equations in Bledjian's formulation [55]. This consists of converting the set of partial differential equations to a set of ordinary differential equations by discretizing all but one of the independent variables at a time, the resulting equations being then integrated using RungeKutta's fourth order method. Difficulties were found when the set of equations was stiff. The step size for the next integration cycle was estimated based on the local truncation error, as calculated by using Richardson's extrapolation formulas. Calculations were carried out for two flames: the hydrazine decomposition flame and the ozone decomposition flame. The results for the first case agreed well with Spalding's calculations [45] and for the second case good agreement was obtained with experimental results. Estimations of the minimum ignition energy for the two cases have also been obtained. For the sake of computational stability, the time step size was restricted by the condition (X/pC )At/(A,)2 1/30 instead of the standard stability condition (X/pC )At(ALp)2 1/2 for the heat conduction equation. condition REFERENCES FOR CHAPTER III 1. Evans, M. W., "Current theoretical concepts of steadystate flame propagation," Chem. Rev., 51, 363 (1952). 2. Williams, F. A., "Combustion Theory," AddisonWesley (1965). 3. Zeldovich, Y. B. and FrankKamenetski, D. A., "Theory of uniform flame propagation," J. Phys. Chem. (U.S.S.R.), 12, 100 (1938). 4. Sokolik, A. S., "Selfignition, flame and detonation in gases" (Translated from Russian), Jerusalem (1963). 5. Von Karman, T., "The present status of the theory of laminar flame propagation," Sixth (Int'l) Symp. on Combustion, Reinhold Publishing Corp., New York (1956). 6. Boys, S. F. and Corner, J., "The structure of the reaction zone in a flame," Proc. Roy. Soc. London, A197, 90 (1949). 7. Corner, J., "The effect of diffusion of the main reactants on flame speeds in gases," Proc. Roy. Soc. London, A198, 388 (1949). 8. Adams, E. N., Rept. CF957, Univ. of Wisconsin Naval Research Lab. (1948). 9. Wilde, K. A., "Improved approximate solutions of flame equations governed by simple chemical reactions," J. Chem. Phys., 22, 1788 (1954). 10. Hirschfelder, J. O., "Theory of flames produced by unimolecular reactions," Phys. Fluids, 2, 565 (1959). 11. Spalding, D. B., "I. Predicting the laminar flame speed in gases with temperature explicit reaction rates," Combustion and Flame, 1, 287 (1957). 12. Rosen, G., "An action principle for the laminar flame," Seventh (Int'l) Symp. on Combustion, Butterworths Scientific Publications, London (1959). 13. Rosen, G., "Generalization of the laminar flame action principle for Arrheniustype rate functions," J. Chem. Phys., 32, 311 (1960). 14. Adler, J., "The prediction of laminar flame speeds in stoichiometric mixtures with nonnormal diffusion," Combustion and Flame, 9, 273 (1965). 15. Spalding, D. B., "II. Onedimensional laminar flame theory for temperature explicit reaction rates," Combustion and Flame, 1, 296 (1957). 16. Adler, J., "The limits of the eigenvalue of the laminar flame equation in terms of the reaction ratetemperature centroid," Combustion and Flame, 3, 389 (1959). 17. Favin, S. and Westenberg, A. A., "The theory of a spherical premixed flame," Combustion and Flame, 4, 161 (1960). 18. Spalding, D. B. and Adler, J., "Flame propagation by non branching chain reactions," Combustion and Flame, 5, 123 (1961). 19. De Sendagorta, J. M., "Propagation velocity and structure of single step reaction flames," Combustion and Flame, 5, 305 (1961). 20. Von Karman, T. and Penner, S. S., "Selected Combustion Problems," AGARD, Vol. I, Butterworths, London (1954). 21. Millan, G. and Da Riva, I., "Distribution of radicals in laminar flames," Eighth (Int'l) Symp. on Combustion, Williams and Wilkins, Baltimore (1962). 22. Spalding, D. B. and Adler, J., "Onedimensional laminar flame propagation with an enthalpy gradient," Combustion and Flame, 4, 94 (1960). 23. Adler, J. and Spalding, D. B., "Onedimensional laminar flame propagation with an enthalpy gradient," Proc. Roy. Soc. London, A261, 53 (1961). 24. Mayer, E., "A theory of flame propagation limits due to heat loss," Combustion and Flame, 1, 438 (1957). 25. Spalding, D. B., "A theory of inflammability limits and flame quenching," Proc. Roy. Soc. London, A240, 83 (1957). 26. Adler, J., "Onedimensional laminar flame propagation with distributed heat losses: Thinflame theory," Combustion and Flame, 7, 39 (1963). 27. Adler, J. and Kennerley, J. A., "Steady laminar propagation with conductive heat loss," Combustion and Flame, 10, 191 (1966). 28. Chen, T. and Toong, T., "Structure and propagation of laminar flames near a heat sink," Combustion and Flame, 4, 313 (1960). 29. Yang, C. H., "Burning velocity and structure of flames near extinction limits," Combustion and Flame, 5, 163 (1961). 30. Adler, J., "Downstream boundary conditions of nonadiabatic onedimensional laminar flames," Combustion and Flame, 11, 85 (1967). 31. Klein, G., "A contribution to flame theory," Phil. Trans. Roy. Soc. London, A249, 389 (1957). 32. Johnson, W. E. and Nachbar, W., "Deflagration limits in the steady linear burning of a monopropellant with application to ammonium perchlorate," Eighth (Int'l) Symp. on Combustion, Williams and Wilkins, Baltimore (1962). 33. Friedman, R. and Burke, E., "A theoretical model of a gaseous combustion wave governed by a first order reaction," J. Chem. Phys., 21, 710 (1953). 34. Giddings, J. C. and Hirschfelder, J. 0., "Flame properties and the kinetics of chainbranching reactions," Sixth (Int'l) Symp. on Combustion, Reinhold Publishing Corp., New York (1956). 35. Menkes, J., "A contribution to the theory of laminar flames with radial symmetry," Combustion and Flame, 4, 1 (1960). 36. Menkes, J., "Analytical solution of a flame with cylindrical symmetry," Eighth (Int'l) Symp. on Combustion, Williams and Wilkins, Baltimore (1962). 37. Jain, V. K. and Ebenezer, J. S., "Analytical solutions of a steady cylindrical flame," Combustion and Flame, 10, 391 (1966). 38. Spalding, D. B., "The theory of steady laminar spherical flame propagation: equations and numerical solution," Combustion and Flame, 4, 51 (1960). 39. Spalding, D. B. and Jain, V. K., "The theory of steady laminar spherical flame propagation: analytical solutions," Combustion and Flame, 5, 11 (1961). 40. Spalding, D. B., Jain, V. K., and Samain, M. D., "The theory of steady laminar spherical flame propagation: analog solutions," Combustion and Flame, 5, 19 (1961). 41. Jain, V. K. and Kumar, R. N., "Thegry of laminar flame propagation with nonnormal diffusion," Combustion and Flame, 13, 285 (1969). 42. Hirschfelder, J. O., Curtiss, C. F., and Campbell, D. E., "The theory of laminar flame propagation IV," J. Phys. Chem., 57, 403 (1953). 43. Campbell, E. S., "Theoretical study of the hydrogenbromine flame," Sixth (Int'l) Symp. on Combustion, Reinhold Publishing Corp., New York (1956). 44. Campbell, E. S., Heinen, F. J., and Schalit, L. M., "A theoretical study of some properties of laminar steady state flames as a function of properties of their chemical components," Ninth (Int'l) Symp. on Combustion, Academic Press, New York (1963). 45. Spalding, D. B., "The theory of flame phenomena with a chain reaction," Phil. Trans. Roy. Soc. London, A249, 1 (1956). 46. Zeldovich, Y. B. and Barenblatt, G. I., "Theory of flame propagation," Combustion and Flame, 3, 61 (1959). 47. Lewis, B. and Von Elbe, G., "Combustion, flames and explosions in gases," 2nd Ed., Academic Press, New York (1961). 48. Adams, G. K. and Cook, G. B., "The effect of pressure on the mechanism and speed of the hydrazine decomposition flame," Combustion and Flame, 4, 9 (1960). 49. DixonLewis, G., "Flame structure and flame reaction kinetics I. Solutions of conservation equations and applications to rich hydrogenoxygen flames," Proc. Roy. Soc. London, A298, 495 (1967). 50. DixonLewis, G., "Flame structure and flame reaction kinetics II. Transport phenomena in multicomponent systems," Proc. Roy. Soc. London, A307, 111 (1968). 51. DixonLewis, G., "Flame structure and flame reaction kinetics V. Investigation of reaction mechanism in a rich hydrogen+ nitrogen+oxygen flame by solution of conservation equations," Proc. Roy. Soc. London, A317, 235 (1970). 52. Spalding, D. B., Stephenson, P. L., and Taylor, R. G., A calculation procedure for the prediction of laminar flame speeds," Combustion and Flame, 17, 55 (1971). 53. Patankar, S. V. and Spalding, D. B., "Heat and mass transfer in boundary layers," 2nd Ed., Intertext Books, London (1970). 54. Wilde, K. A., "Boundaryvalue solutions of the onedimensional laminar flame propagation equations," Combustion and Flame, 18, 43 (1972). 55. Bledjian, L., "Computation of timedependent laminar flame structure," Combustion and Flame, 20, 5 (1973). CHAPTER IV MATHEMATICAL ANALYSIS OF THE COMBUSTION PROCESS 4.1 Formulation of the Problem In this chapter the case of laminar premixed flame propagation in a stoichiometric mixture of methane and air, initially at a pressure of 10 atmospheres and at a temperature of 628 degrees K, which approx imately is the condition in an Otto cycle engine at the time of ignition, has been considered. The timedependent approach of Spalding [1] has been used here. The reason methane has been considered is because of the availability of the reaction kinetics data. The method is general, and if the reaction mechanism and kinetic data for the oxidation of any other fuel are known, the same method can be used to obtain the detailed flame structure. The variations of the transport and thermodynamic properties with temperature, pressure and chemical composition have been accounted for as far as the available data permit. 4.1.1 The conservation equations The flame zone is generally a few hundred mean free paths thick, and so the concepts of continuum analysis may be used here. The general conservation equations for multicomponent reacting ideal gas mixtures are given in various references [2,3]. These equations are; (1) Overall continuity: t + V. () = 0 (1) 38 Momentum: SN + (vV)v = (Vp)/p + g Y.f. i=l Energy: N aU  > V p + p*vVu = Vqj p: (V) + p Z Y .iV. i=l Continuity of species: 5Y. w. i i  + vVY. [V(pYV)]/ , 1t i p ( )] i=1,2,..... N In addition to the above conservation equations, other equations describing the fluxes of mass and energy, the chemical kinetics of the system and some constitutive relations are necessary. These are (again in their general form): (1) Diffusion in multicomponent mixtures is governed by: N N X.Xij ij + VX. .(V V.) + (Y. X.)E + i j i j 1 j=l Di. j jD i 1 1 p j A =1 N X ff (ff D% (VT] P Di Y Y J=l j 3 i=1,2,......... .N N Y.V = 0 i=l The derivation of equation (5) assumes that the medium under consideration can be treated as a dilute gas. However, there is and experimental evidence that this result also applies quite well in dense gases and liquids [3, p. 718]. (2) Heat flux vector: N N N rX.D q = VT h p h.Y.V. + R T (Vi V.) (7) I I i u j=l WiDij The same comment as above applies to this equation also. (3) Chemical kinetics for a set of M elementary reactions between N species: SR P M R N tXpn Jk N X p k k=l ( 2 j=l kN kbk N l u ' k=1 u (8) (4) Thermally perfect gas: P = pR T/W (9) u N S= X X.W. (10) i=l1 (5) Caloric equatiai of state: T h. = h + C .dT i=1,2,........N (11) S 1i p, ' T Implicit in these equations is the assumption that every point in the flow field has a thermodynamically definable state, i.e., it is possible to attach significance to the concepts of local temperature, pressure, density and chemical composition (see Sec. 4.4 for un certainties in the definition of temperature). 4.1.2 Major assumptions for analysis The general equations are now simplified for the case of one dimensional laminar flame propagation in spherical coordinates using the following assumptions. (1) The kinetic energy associated with the motion of the gases is small compared to the heat of the reaction and may be neglected. Also, since the velocities and velocity gradients are small, viscous effects may be neglected. Thus the momentum equation can be replaced by the assumption that the pressure is constant through the flame. In fact, solving the momentum equation for the pressure change across the flame yields: V M S Ap = P, P = p M ; M = (12) u u From experimental observations it is seen that in general M < 0.02, V /V < 15, and y < 1.667, which yields Ap 0.01 p . u b u u Hence for ordinary flames the process may be considered essentially isobaric. This assumption would not hold for detonation waves. (2) Radiative heat transfer is negligible. Although relatively high temperatures ( 25000K) are encountered in flames, radiation is not important since the gases have low emissivity and absorptivity and so the radiative heat transfer during the reaction time is small compared to the total energy flux. With flames containing solid carbon this is generally not the case. (3) The Soret and Dufour effects are negligible. (4) Any diffusion caused by pressure gradients is negligible. Since assumption (1) already considers the process to be isobaric, this particular assumption is not necessary. This term may be significant in strong detonation waves wherein tremendous pressure gradients may be established. (5) External body forces are neglected. For very hot flames wherein ionization occurs to a significant degree, this term may have to be considered if an external electric field is also present; the external force on an ion is equal to the product of the ionic charge and the local electric field strength. This gives rise to what is termed as 'forced diffusion'. (6) The chemical kinetics is assumed to follow Arrhenius equation. Thus: kk A= kTb exp(E /RT) (13) The backward reaction rate constant kbk is then obtained from a knowledge of kfk and the equilibrium constant KCk: kfk K k (14) Ck kbk There is considerable controversy [4,5] regarding the validity of equation (14); it will nevertheless be used in the absence of a more accurate relation. (7) The gas mixture is thermally perfect, i.e., it obeys the ideal gas equation. Nitrogen is a major constituent (> 70%) of the mixture and it closely approximates thermally perfect behavior for the conditions considered in this study. Other constituents are present in relatively smaller concentrations, so as not to affect the overall thermally perfect behavior. The equations finally simplify to: (1) Overall continuity: D + (r2pv ) = 0 (15) t 2 Br r r subscriptt r denotes the radial direction) (2) Energy equation: 9T DT 1 1 3 2 DT 1 N h i t r Dr pC 2 3r r2 r C i ir ri p r p i=l N 1 C h.w, (16) pC iJ i i p i=l (3) Species continuity: BY. BY. w. S+Vr rv 1 2 (r ) i=1,2,..... N (17) 9t r r p p 2 r i ir r (4) Diffusion equation: ax. N X.X ( V. ) i=1,2,......N (18) Dr jjr ir J=l I j ; I N and Y.V. = 0 (19) j jr j=1 3 Jr (18) represents a set of (Nl) linearly independent equations, since N X. Sar i=l1 The unknowns in these equations are: Y., Vir, T, v X. and P are determined once Y. and T are known respectively, using the following two equations: Y./w. X = 1 1 i=1,2,.......N (20) I (Y./W ) j=1 and p = p' (the ideal gas equation) (21) R T u 4.1.3 Boundary conditions The next step is to set up the boundary conditions at the hot and cold ends for the set of partial differential equations given above. a. lne cold boundary difficulty. This has been discussed considerably in the literature. Basically the difficulty arises from the incompatibility of the energy conservation equation with the boundary condition it is subjected to at the cold end. The model for the flame which assumes both 'steady state' and 'adiabatic' conditions does not represent truly the physical processes occurring and these conditions contradict each other for a realistic reaction rate law, wherein the reaction rate at the cold boundary is small but finite [6]. To satisfy the conditions of both adiabatic and steady state requires any one of several assumptions. The assumptions are mathematically equivalent but have different physical meanings. Also the validity of these assumptions on physical grounds is open to question, although these assumptions have been used widely to solve the flame equations. The assumptions are: (1) This relaxes the adiabatic condition at the cold boundary by allowing the existence of a small but finite heat sink at the boundary in the form of a flame holder. It was introduced originally by Hirschfelder et al. [3]. This was done because a solution does not exist when the very small heat loss at the cold boundary is neglected; the neglect of heat transfer (due to lateral conduction and radiation) at other points in the flame zone does not, however, affect the existence of a solution. Mathematically, this is equivalent to relaxing the condition (T/Dr) = 0 at the cold boundary, and instead allowing for (3T/r) = E, where E is a small positive quantity. The value of E determines the 'strength' of the flameholder in terms of the heat abstracted by it and has to be greater than a certain critical value to permit the existence of a solution. This assumption is physically unrealistic for a freely propagating flame which does not have a flame holder or heat sink propagating with it. (2) The second assumption compromises another physical aspect by allowing for the existence of a temperature (greater than the initial temperature of the reactants) below which the heat release due to chemical reactions is identically zero. This temperature has widely been referred to as the 'ignition temperature' [7]. The choice of this temperature is again not arbitrary but depends on the particular flame being studied and it can be determined from the governing equations. It permits the adiabatic assumption (3T/3r) = 0 at T = T to be satisfied; however, now this occurs at r = o, as u required for a strict mathematical model [6]. The use of the ignition temperature concept can be explained by relaxing the adiabatic condition. It is assumed that the heat loss from the cold boundary and the heat released by chemical reactions (not identically zero) at the cold boundary reach an equilibrium at the temperature T Further, assuming the rate of the chemical reaction to be extremely slow (the initial mixture is in a retastable state), the change in concentration with time at the cold boundary can be neglected. Thus it is as though no reaction has taken place, the heat released at the cold boundary having been conducted away. (3) The third assumption involves the use of truncated kineticswherein the reaction rate goes to zero as the temperature approaches its cold boundary value. Tlis has been shown to be equivalent to the use of an ignition temperature which is very close to, but not equal to the cold boundary temperature. Again, this assumption is not physically possible. However, for most systems of interest, the reaction rate at the cold boundary is so small as to be immeasurable, and so it can be considered to be essentially zero. The reaction rate law is then written as kf = ATn exp(E/RTbT) (22) T T where T = u (23) b u as compared to the actual law, k = ATn exp(E/RT) (24) The factor T varies from zero at the cold boundary to unity at the hot boundary. The change in magnitude of the reaction rate due to this modification becomes smaller as the temperature approaches its hot boundary value and this is where most of the chemical reaction and its associated heat release occur. In the calculations carried out here for the methaneair flame with an initial temperature of 628 degrees K, it was found that the chemical reaction rate at this temperature was lower than the smallest number that can be handled by the IBM 370/165 computer [i.e., exp(176)]. In fact, the reaction rates were ignored at temperatures below 650 degrees K for all the elementary reaction steps, and above 650 degrees K, the reactions for which E/RTbT was greater than 120.0 were ignored. The condition E/RTb T 120.0 was used instead of E/RTbT 2 176.0 to avoid the occurrence of a number smaller than exp(176) later during the computations. Again, when integrating from the hot to the cold boundary, the integration was stopped arbitrarily for (Fk Fkl)/Fk j5 0.005; (F = dependent variable; temperature and mass fractions) the subscript 'k' refers to the grid point in the rdirection. This had to be done since, when carrying out numerical calculations, the condition (8F/ar) = 0 can strictly be achieved only at r = c. b. Hot boundary. When only the initial (unreacted) and final (reacted) states of the gas are considered, it is seen that the gas has passed through the flame zone without any net gain or loss of heat. All that has occurred is a change in composition due to chemical reaction and all the heat of the reaction has gone to raise the temperature of the gas. This adiabatic condition is not fulfilled for flames which lose significant amounts of heat by other processes, e.g., radiation, heat losses to flameholders, etc. However, within the framework of the assumptions used in this analysis, the adiabatic assumption is satisfied. Again, the adiabatic condition is not satisfied for any element of gas within the flame zone, since it loses or gains heat by conduction and diffusion of species. Therefore, for a constant pressure flame, since the heat transfer equals the change in enthalpy, the enthalpies (per unit mass of mixture) of only the initial and the final states are equal. For the special case of unity Lewis number, the constant enthalpy condition is satisfied throughout the flame zone. Since the heat fluxes due to diffusion and thermal conduction are equal in such a case, any element of gas within the flame zone does not suffer any net gain or loss of heat and so goes through the flame zone at constant enthalpy. The condition of equal enthalpy at the hot and cold ends is utilized to obtain the boundary condition at the hot end. Since the gas mixture is in chemical equilibrium at the hot end, the concentra tions of the various species are obtained by simply calculating the equilibrium concentrations corresponding to an isenthalpic combustion process, starting with a mixture of given concentrations and state properties at the cold end, i.e., h(Tb,p) = h(Tu,p) (25) which is known from the initial conditions at the cold end. Thus N h(Tu,p) = C (hXiu )TT (26) i=l u where the h. are obtained from the JANAF Tables. 1 Considering a Taylor series expansion about Tbk, the kth estimate for the value of the burnt gas temperature, and neglecting higher order terms, h(Tp) = h(Tbk) + (Tb Tbk) (h (27) where the partial derivative is evaluated at (Tbk,p) and N h(Tbk,P) = I (hXi )T=T (28) i=l b b bk h(Tb,p) h(Tbk,p) T = p(Tbk) (29) P k Tbk+l = Tbk + ATbk (30) A first estimation of the temperature at the hot end is obtained from: AhX Ah fuel T = T + (31) b u C (T) (31) P u where Ah = heat of combustion (kcal/gm.mole) C (T ) = specific heat at constant pressure of the mixture at Pu temperature T (kcal/gm.mole OK) N = Cpi (Tu)Xi (32) i=l u The Ci are again obtained from the JANAF Tables. The equilibrium concentrations of the species are calculated at the estimated value of the burnt gas temperature by the method given in Appendix I. Then equation (30) is used to obtain the next estimate of this temperature and the process is repeated till convergence is obtained to within a predetermined accuracy: A11 T (33) ATbk Tbk  (33) = 0.001 was used in the calculations. This determines the boundary conditions at the hot end: T = T (34) X. = X. i=1,2,.......N i ib 4.2 Method of Solution An explicit finite difference method was used for numerically solving the partial differential equations (1519), and obtain the profiles of the state properties and the species concentrations through the flame zone. A scheme involving a symmetrically staggered mesh in the rt plane was used. This method has been used by Nash [8] for the solution of the threedimensional boundary layer equations. The staggered mesh scheme uses a set of secondary points in addition to the primary collocation points, the former lying at the centers of the rectangles formed by the latter, as shown in Fig. 4.1. In this method, values of the dependent variables are at first calculated at the secondary points [at time (t+At/2)] using the values associated with the primary points at time t Then using values associated with both the primary points at time t and the secondary points at time (t+At/2), values of the dependent variables at the primary points at time (t+At) can be obtained. This is more accurate than proceeding directly from time level t to level (t+At). The difference equations for the two steps are readily obtained by double Taylor series expansions about the appropriate points, and are [8]: 1orI( c 2 (F + F ) + {Ar + T AtM + B + e c 2 A B [3A B (35) ~7~N 9 At/2 FAr Ar/2 t  Fig. 4.1(a). Staggered mesh scheme, rt plane t T(i At  0 o primary points o secondary points E(i,k+l) 1 1 C: (i + k + ) 2D: (i 2 i 1 D: (i k + ) 2' 2 Fig. 4.1(b). Enlarged view of staggered mesh T' At _L, A(i,k) B(i+l,k) .Ar II. I  I ^ I C I py ~   .~    1I    FE = FA +I (1p)(FD + FC) + X(lp)Ar L  SD D l cJ + (l+p)At D+ J + e2 (36) Here F denotes any one of the dependent variables. The quantities (8F/)t) are obtained from the governing partial differential equations (16) and (17), the rderivatives in these equations and in equations (35) and (36) being expressed in finite difference form. A fourth order dif ference formula is used for this purpose: lar1 2 F 1 (F F ) ( F ri, 12Ar i2,k i+2,k 3Ar il,k i+l,k 1 4 5F + I (Ar)4 (37) or In equations (35) and (36), el and e2 denote the discretization or truncation errors, and contain higher order terms of the double Taylor series expansions. e. = ( X )(Ar)2 2 (At)2 t + At(Ar)2 + (38) 1 2 d2F 1 2 D2F e = (l X) (Ar)  (l0)(At) 1 2 8 F + (4X4Xh + p)At(Ar)2 + ...... (39) 8 tr2 X and p are constants which govern the accuracy and the stability of the computation. The pair of values X = 1/8 and p = 0 serve to minimize e2 and also assure stability [8]. With these values the twostep scheme reduces to a onestep scheme: FC = (FA + FB) + Ar B + At (40) FE = (FD + FC) + Ar D C + At D (41) (FDD C DC F(r) = Fu r > ru 4.2.1 Initialization Sshaped profiles for the dependent variables Y. and T were assumed to start the integration. Second order curves were chosen and these were of the form: F(r) =F r r u u (r r)2 F(r) = F + (F Fu) u r r rb +6 26 b2 (42) (rrb) F(r) = F (Fb Fu) rb+ 6 r > rb 26 F(r) = F rb r An arbitrary value of 6 = 2 (r rb) = 0.005 cm. was chosen. It is immaterial what value of 6 is chosen since the final profiles are independent of the starting ones. 4.2.2 Details of computation The initial profiles having been fixed, the integration process can now be started. Expressing (DX./Or) in finite difference form in 54 accordance with equation (37), the diffusion equations (18) and (19) reduce to a set of N simultaneous algebraic equations with the N diffusion velocities V. (j=l,2,....N) as the unknowns. These equations were solved using GaussJordan elimination with maximum pivot strategy. Thus the diffusion velocities were obtained at each grid point for a given time level t. Before proceeding further, two other methods for the determina tion of diffusion velocities will be mentioned. These methods are approximate, but involve lesser computation than the method mentioned above, and may therefore be used where conditions warrant their use. (1) In flames involving air as oxidant, one component (namely, nitrogen) is present in great excess so that it is possible to simplify the problem by regarding each of the other species, one at a time, as forming a binary mixture with the main component. In general, the effective binary diffusion coefficient D. for the diffusion of species i in a im mixture m is given by [9]: N S(1/nD)(X.N. X.N.) 1 _= 1 1 . 1 1 (43) nD. N mN. X. ) N. 1 j=l For the diffusion of trace component i=1,2,3,...N (ifj), in nearly pure species j, D. z D. im li and so the molar fluxes and hence the diffusion velocities are obtained by solving the set of algebraic equations: X. N N. = nD. X. N (44) 1 im ar 1 (44 j=1 Ni = n Vi (45) (2) This approximation may be used when species i is still a trace, but there is no single component which may be considered to be in excess compared to all the others. For such a case the diffusion velocity of the trace species is obtained from: dX 1 i ri dr Vi = 1dr (46) X. Xj/D. Sj i J j which shows that the quantity X./Dij is an effective diffusion jfi coefficient for the species i and the mixture. The above equation is not valid for the diffusion of a species present in significant amounts. Such species may be treated in two ways. Often a species which is present in excess is not significantly affected by the chemical reaction (e.g., nitrogen) so that its concentration is more or less constant through the flame zone and hence its diffusion velocity may be neglected. On the other hand, for species which are chemically reactive, and not present as a trace, the diffusion velocity may be computed by difference from equation (19), which is N SYiVi = 0 i=l Getting back to the main problem, it is seen that the mass average velocity v is as yet unknown, and a suitable profile for it was assumed subject to vru = 0, i.e., the unburnt gases are at rest. This need be done only for integration over the first time step; for succeeding time steps, the values of vr from the preceding time step were used for a first approximation. With this assumed profile for Vr, (3T/at) and (3Yi/3t) (i=1,2,....N) were calculated at each grid point for a given time level t, using equations (16) and (17). The derivatives in r occurring in these equations were approximated, again by the fourth order difference formula (37); values of Vjr (j=1,2,....N) as calculated above were used, and the JANAF Tables were used to obtain the values of the state properties h and C The w. (i=1,2,....N) were p 1 calculated according to equation (8) and modified by truncating the reaction kinetics as discussed in connection with the cold boundary difficulty. Using the ideal gas equation, it may be shown that: SP, N ti (47) P W 1 (47) Si=i i and thus the overall continuity equation may be rewritten in the form: Sr 1 Y. 1 2 v 2 lp + T tJ r dr (48) pr i=l i r u (for given time t) satisfying r = 0. The values of (DT/at) and (gY./8t) calculated above with an assumed profile for v were now used in equation (48) to obtain better values for vr, which were then used to obtain improved values of (8T/at) and (OY./8t). The iteration can be terminated when the vr values converge to within a predetermined accuracy. Actually the values of (DT/8t) between two successive iterations were compared for checking convergence. This was done since (T/8t) is used to determine the temperature at the next time level (equations 40, 41) and so convergence of (8T/t) is of greater importance than that of v which is not used anywhere else in the calculations. The con vergence of (8Y./ t) (i=1,2,....N) may also be tested, but these were not found to vary significantly and the convergence of (aT/9t) generally assures the convergence of (gY.//t) for all i. A convergence criterion of the type At {[T] (49) ( k _k1 was used with E = 0.0005. This assures the errors in the calculated temperatures due to uncertainty in the term (3T/8t) to be less than E/2 each time either equation (40) or (41) is used. The subscripts in the above convergence condition refer to the iteration number and this condition is to be satisfied at each grid point in the flame zone before proceeding further. Equation (48) involves quadrature; this was done using the trapezoidal rule with the same space interval Ar as used for the other computations. Simpson's rule along with Newton's 3/8 rule may also be used (IBM subroutine QSF). The difference between the two was found to be negligible since the space interval Ar is small. Hence the trapezoidal rule which needs less computer running time was used. The term involving (8Y.i/t) in equation (48) has only a small influence on the final profile of v and may be neglected to a first approximation. This would amount to making the assumption that the average molecular weight of the gas mixture is constant through the flame zone. This is approximately true for methaneair flames in which a single constituent, nitrogen, accounts for more than 70% of the concentration at any point in the flame zone. For the sake of completeness, this term was included in the calculations carried out here, although its effect is small. The iterative scheme used here for the solution of equation (48) is similar to Picard's method of successive approximations for the solution of a Volterra integral equation [10]. This is seen from the analysis given below. Equations (16) and (17) may be written as: T DT = r + G(r) (50) ot r or DY. DY. v + H. (r) i=1,2,.....N (51) at r ar r where ( ) 1 1 @ r 2 X BTj pC 2 ar (r2 G(r) pr p r 1 C P ah. YiV i ar 1 r ar N i I h.w. pCp 1 1 p i=1 H.(r) = i (r2 pY.Vi ) i=1,2.....N 1 p p 2 Dr 1 r r (53) At a given instant of time the right hand sides of (52) and (53) are functions of r only. Equation (48) now reduces to: r 1 f pr2 r vpr r u N Y. i=l 1 1 T N H.(r) +\  W. i=l 1 + (r) dr T (54) Again, for a given instant of time, N aY V N 1 _ _W. Dr i=l 1 N H.(r) w. i=1 1 1 aT T ar G(r) T 2 pr = A (r) SA (r) = A3(r) (55) where A A2 and A3 are functions of r only. So (54) finally reduces to: r vr(r) = A3r) A3(r) [vr(r)Al(r) + A2(r)]dr r u (52) (56) 60 This is nothing but a Volterra integral equation of the second kind, and the iterative technique described above is Picard's method of successive approximations. The first approximation to v (r) in this method is given by: r rv(r) A3(r) A3(r)A2(r)dr r u 1 0K { 1 i 1 (r YiVi} Pr Ii=l pr r  N N 3h N + 1 1 1 V i 1 2 p r p i=l p i=l (57) Having finalized the values of the time derivatives (DT/Dt) and (OYi/t) at a given time instant t, equation (40) was used to proceed to the secondary points at the time level (t+At/2) and by a similar process the values of the dependent variables at time level (t+At) were obtained using equation (41). The mole fractions were then calculated from the mass fractions using equation (20) and the mean gas density was obtained from a knowledge of the temperature and the mean molecular weight using the ideal gas equation. Calculations in the rdirection at the cold boundary were stopped at any given time level when it was found that the relative changes in the temperature and the species mass fractions were less than a predetermined value, i.e., F Fk k F (58) k E = 0.005 was used in the calculations. This process was repeated until the dependent variables converged to their steady state values. The transport properties in the conservation equations were determined (Appendix II) at each grid point every 75 time steps; the error introduced by not calculating the new values of the transport properties at every time step would cause small errors in the transient state solutions, but would not affect the steady state solution. However, since the estimation of the transport coefficients is quite often approx imate, the error introduced by not calculating these coefficients at every time step, might in fact be lower than the error involved in the calculations based on the approximate theories. The last step involves the computation of the laminar flame velocity based on the unburnt gas density. This was done by considering a control volume which is large compared to the thickness of the flame zone, and equating the net rate of production of species i due to chemical reactions to the rate of outflow of species i. The result is: rb f w.dr r u = (59) uu Y. Y. lu ib The values of the flame velocity, S as calculated from the above equation converged as the state properties and species concentra tions approached their steady state values only when species i was considered to be one of the major species (i.e., CO2, H20, CH4, 02); however, for each one of the species a different value of the flame velocity was obtained. The reason for this are the unreliable kinetic data. Further, this type of behavior has been observed before [11]. The flame velocity calculations based on the other species did not converge to a steady value, presumably because the term (Yiu Yib) involved in the denominator of equation (59) is small for these species and small uncertainties in this term could cause large errors in the calculation of S These values might, however, converge if u computations are carried on further. In the calculations reported in Ref. 12, it was found that the flame velocity based on the unstable species needed more time steps to converge to a steady value. 4.2.3 Step size and stability The Sshaped profiles used initially were at least qualitatively correct for the major species (CO2, CH 02, and H20). However, this was not so for the intermediate species which have a maximum concentra tion somewhere in the flame zone and decrease towards both the hot and the cold boundaries. For this reason initial time steps were restricted to a very small value until the profiles of the intermediate species approached their correct form, after which larger step sizes could be used. A starting step size of At = 1.OE08 sec. was used and this was doubled every 50 steps until it reached a maximum value of 32.0E08 sec. The step size had to be restricted to this maximum value to satisfy the 2 stability condition (X/pC )At/(Ar) < 1/2. Since X, p, and C vary from point to point, the combination of values of and Cp which from point to point, the combination of values of A, p, and C which would give the least upper bound for At was used to show that At < 1.28(Ar)2 which gives At < 32.0E08 sec. with Ar = 5.0E04 cm. max max Thus this explicit scheme is conditionally stable. The use of these values of At, however, caused instability in the integration of the term involving the chemical rate of production of species. So this term was integrated separately with an appropriate value of step size. The integration can be done by RungeKutta's fourth order method for a set of first order differential equations. The equations i. dn. dt fi(nn2'.....N) i=1,2,.....N (60) form a set of what is termed as 'stiff'* equations, and such equations need very small step sizes and hence prohibitively long computer running times with RungeKutta's method. Gear [13] has developed an algorithm for dealing with stiff equa tions, and this was used in the integration of the reaction rate equations (details of the method in Appendix III). Instead of starting the integra tion by a first order method, as proposed by Gear, a higher method was used. This necessitated the calculations of up to six derivatives of ni, since the maximum order for which the algorithm is devised is six. Generally, it was found in this problem that a starting fourth order method resulted in the maximum value of the step size. The order for subsequent time steps was decided by the program itself. The following equations were used to obtain the derivatives. So called because this behavior was first noticed in connection with servomechanisms having a tight coupling between the driving and the driven member. dn. M di (i i)[Pj P2j] dt 1 i 1i 1i 2j j=1 d n. M 2 (v.. v R)[P P P .] dt2 j=l 1j Ij3j P2jP4j dt j=1 d ( P R +2 p 2 P t j (1 i ii 3j lj 5j P2j4j 2j6j dt j=1 d n. M R 3 SVP R p3 7= (.. .)[P .P + 3P .P .P + P .P 4 13 13 ij )[Plj 3j lj3j5j + 7 dt j1  P p3 2j 4j  3P 2P .P 2j 43 6j (60a) (60b) (60c) (60d)  P2jPj ] 5 d n. 4 2 2 __ ( P )[p P + 6P P P + 3P P + 4P 2P dt5 j=1 i j i lj3j l+ 3j 5j + 53I j + lj 3jP7j dt j=1 4 2 2 + .P p P.4 6P.P 3 2 P 4P2.P.P. lj j P2j 4j 6P2j 46j 3P2j 6j 2jP4jP8j  P2jP j I 2j 10j (60e) d n. M i P 5 3 2 1t 16 j 1 3j 1j 3j 5 j j3j 5j dt j=l + o10P .P2 + 10P .P .P. + 5P .P. + P .P 1 P 3j 7j + 1P 5j7j jj + PjP 1j 5 3 2 P.P 10P.P .P 15P .P .P 2j 4j 2j 4j 6j 52jP4j 6j O1P p2 P 5P P(60f) 2j 4j 8j 2jP6j 8j 5P2j lOj 2j 12j N Pj = kfj i= N P2 = kb 1 2j i=l N 3j r=1 N P = P5j E =l N 7j =l N P9j = zj Q4Z k=i N Plj= 1 R kj li R Vj Q2Z R j Q3Z N (60i) P4j = 1 Q N (60k) P6j = Z j Q2 N (60m) P8 = OJ R=1 P V j Q3 N (60o) P10j= k 1 j Q4 Rj k JQ5e 1 dn QI n dt N (60q) Pl2j = I 3 =i (60s) Q2Z 1 2 "Jl P vj Q5Z 2 1 d2n n dt2 dn 3 3 dn dt J 2dt dn] 2 d j Idn 4 dt . 4 S1 d 4 "n dt4 where R V.. [n ij V. . [ni] j (60g) (60h) (60j) (60( ) (60n) (60p) (60r) (60t) d2 dt2 2 Q3~, 3 nk 6   'U 3 1 d3n Sn dt3 (60u) d2n d t2 2 dt 2 2 dt 2 3 2 n d3n 3 d 3t dt3 4 dn 2 dt n (60v) dn} 2 dt ) 12 3+ n k 24 dn 5 60 dnl 3 d2n 30 dn d n 2 20 d 2 d3n S0 + +20 5 5 dt 4 dt J 2 3 d dt 2 3 dt 3 d 2 n d nn d dt n dt n dt 10 dt d3n] 5 dn d4n 1 dn 2 1d J [dt P 2 dt + (60w) n2 3 2 dt 4 n dt5 n dt dt n dt dt The calculation of these derivatives requires a knowledge of the temperature to determine the values of the rate constants involved in these equations. The temperature was assigned its value at the beginning of a time step of the main integration and regarded as constant over this main time step. The computations might be repeated, however, using a suitable average value once an estimate has been obtained for the temperature at the end of the time step. Instead of using a simple arithmetic average T = 1 (T + T2), an average of the 1 nature exp(E/TmR) = [exp(E/RT1) + exp(E/RT2)], which in fact averages the rate constant as opposed to the temperature might be used. (E = activation energy of the rate controlling step.) This iteration across a time step would not significantly alter the results, and so this was not done in this study. The term h.w. in the energy equation was treated similarly. Since the temperature was assigned its value at the beginning of a main time step and regarded as constant over the main time step, the h. (i=1,2,......N) are also constant over this main step and so, N N h.w. dt = h w. dk (61) i=l x i=l So the result of the previous integration was used here instead of resorting to a separate integration of the term X h.w.. Again, a mean value for the h. can be used for improved accuracy. 4.3 Chemical Kinetics The conservation equations contain terms involving the rate of production of species due to chemical reactions. These terms must be evaluated before the equations can be solved. The evaluation of these terms requires a knowledge of the reaction mechanism and the rate constant parameters for the different elementary reaction steps. Rate constants for the elementary reactions occurring in the COHN system have been tabulated by several authors. However, there is considerable disagreement over the values of the rate constants from author to author. Again, there is disagreement over the reaction mechanism proposed for the oxidation of methane. In this analysis a total of 16 species were considered. It was found that it is possible to write a total of 84 (42 forward and 42 reverse) bimolecular elementary reaction steps for these 16 species. In addition to these, several third order reactions involving a third body are also possible. The kinetic data for all these reactions are not available; moreover it is not necessary to consider all the reactions for which data are available. A detailed analysis of all the bimolecular reactions for which data are available was made to determine the contribution of each reaction towards the total rate of production of each of the species. This was done over the temperature range 60024000K. On the basis of this, reactions which contributed less CO2, CO, N2, H20, H, OH, O, NO, NO2, N20, N, 02, H2, CH3, CHO, CH4. than 0.1% to the total rate of production of any of the species over the entire temperature range were eliminated from further consideration. The contribution of each reaction from an energy standpoint was not tested; it being assumed that the reactions which contribute less than 0.1% to the total rate of production of any of the species would not cause any significant contribution to the energy equation also, unless their heats of reaction are at least an order of magnitude more than those of the reactions which were not eliminated from consideration. A similar analysis was also done for the third order reactions. Finally, a mechanism consisting of 25 elementary reaction steps (Table 4.1) was used for the analysis. The first 14 reactions are the ones which Seery and Bowman [14] have suggested in their proposed scheme for the oxidation of methane. The remaining reactions are concerned with the N0 system. The kinetic data used in this analysis have by no means been firmly established, and suggested error limits for some of the rate constants are as high as 25%. The error is greater the higher the temperature, since most of the experimental data are gathered in the temperature range 3001000 degrees K and involve extrapolation to higher temperatures. Also the oxidation mechanism of methane has not been firmly established and much work remains to be done on the effects of pressure, temperature, mixture composition, and the presence of inerts on the oxidation mechanism. In the absence of better data and a more definite reaction mechanism, the data referred to above were used. The method developed however is general, and as 69 Table 4.1 Reactions considered in this analysis k = ATn exp(E/RT) Forward rate constant Reaction E/R CH4 + M = CH3 + H + M 4 3 OH + CHO = CO + H20 0 + CH = HO + CHO CHO + M = CO + H + M 0 + CH4 = OH + CH3 4 3 CH4 + H = H2 + CH3 OH + CH4 = H20 + CH3 4 2 3 OH + CO = CO + H H + 0 = OH + 0 2 0 + H2 = H + OH 0 + H20 = OH + OH H + H20 = H2 + OH M + H20 = H + OH + M 0 + CH = HCO + H N + 0 = NO + N NO + 0 = 0 + N N + 0 = NO + NO NO + 0 = 0 + NO M + NO = 0 + NO + M 2 0 + N20 = NO + NO NO + NO = NO + 02 NO + N2 + 02 = N20 + NO2 OH + M = 0 + H + M H + NO = OH + N N2 + OH = H + N20 1.5x109 1.0x1014 2.0x10 13 2.0xl013 1.7xl013 6.3x1013 2.8x1013 5.6x1011 14 2.24x10 1.74x1013 8.4x104 14 1.0x1014 22 2.0x10 1.0x1014 1.36xl014 1.55x109 4.55x1024 1.0x1012 l.lxlO16 12 3.24x102 6.0xl014 2.30x1016 1.0x1021 5.0x1011 4.0x1014 Units: E/R in degrees K N rate constants in (cc/g.mole) /sec. where N = order of the reaction Ref. 0 0 0 0.50 0 0 0 0 0 0 0 0 1.0 0 0 1.00 2.50 0 0 0 1.50 0 1.50 0.50 0 50300 0 0 14400 4380 6350 2500 543 8450 4750 9120 10200 61500 0 37700 19450 64300 22700 32500 32140 4950 17885 50500 2500 8000 Reaction Table 4.1 (Continued) Reactions which were not considered because of negligible contribution Forward rate constant Reaction A n E/R Ref. N + NO = NO + N 0 N2 + 02 = 0 + N20 NO + H = H20 + NO OH + NO = H + NO2 0 + H = OH + OH 0 +N20= NO2 + N N + M = N + N +M NO2 + NO = NO + NO + 02 NO + N + 0 = N20 + NO2 0 + 0 + M = 0 + M CO + NO = NO + CO2 NO + NO = NO + N NO + M = 0 + N + M NO + M = N +N 02 + M N2 + 0 = NO2 + N 2 2 2 4.7x101 4.xl19 4.0x1019 2.0x109 1.0xl014 14 l.Oxl14 3.7x104 2.97x1021 13 1.0x1013 16 2.3x016 1. 0Oxl14 5.0xl011 1.00x10 17 2.27x107 6.00x1014 2.70x1014 0 0 0 0 0 0.5 1.5 0.5 0 0 0 0 0.50 1.50 1.00 41650 46600 9000 600 35000 20900 112500 15500 17885 0 13900 49300 74900 52600 60600 Reactions for which kinetic data were not available CO2 + N2 = CO + N20 CO2 + OH = 02 + CHO CO2 + N CO + NO CO + OH = 0 + CHO CO + CHi = CH3 + CHO N2 + NO = N20 + N H20 + N = NO + H2 CO2 + H = O + CHO CO2 + 0 = CO + 02 CO2 + 2 = CO + H 0 CO + H2 = H + CHO N2 + H20 = N20 + H2 H20 + 0 = 02 + H2 Most of these reactions are not of much significance from a practical standpoint. Table 4.1 (Continued) more is known about the oxidation mechanism and the rate constants, the new information can be readily incorporated to obtain more realistic solutions. In the computer program developed any number of species and elementary reaction steps may be considered, the only limitation being the computer time and storage available. 4.4 Uncertainty in the Definition of Temperature It is now generally realized that there does exist a distinct possibility that the various chemical species in a flame zone may not be in equilibrium with respect to their various degrees of freedom. Thus the concept of a 'temperature' to define the thermodynamic state of the gases in the flame zone is rendered difficult. This problem is briefly discussed here; details may be found in the appropriate literature. The statistical definition of temperature is based on the MaxwellBoltzmann equilibrium energy distribution function. It is defined as a parameter characteristic of a system at equilibrium and it specifies the distribution of energy between the particles comprising the system. Normally the definition is based on the distribution of the translational kinetic energy. Since polyatomic molecules have in general other modes of energy storage, it is possible to analogously define various temperatures, e.g., Tib' Tro T elec. For a system at complete equilibrium all these rot' elec definitions give a unique value of temperature, the various degrees of freedom being at equilibrium with each other. The kinetic theory definition of temperature for a system composed of different chemical species is expressed in terms of the mean particle energy, averaged over all the species [3, p. 455]: 3 kT= 1 1 n. ( mV2 2 jl 2 i i j=1 Temperature for the other degrees of freedom may be defined analogously. For a multicomponent system not at equilibrium, it is possible that not only are the various degrees of freedom of a particular species not in equilibrium with each other, but that there is nonequilibrium between the different chemical species, i.e., it may not be possible to define a unique translational (or rotational or vibrational) temperature for the entire system. In a laminar premixed flame, the temperature of the gas is raised extremely rapidly (more so in detonations and shock waves). When this happens the different degrees of freedom would in general come to equilibrium at different rates. Translational equilibrium is attained the fastestin about five mean free paths, as demonstrated by shockwave experiments in gases with no internal degrees of freedom (Ne, Ar, Kr). The temperature of primary concern here is the trans lational one, since this is used in the equation of state for transport property calculations and in the evaluation of reaction rate constants. Since the assumption of a continuum has already been made, i.e., changes in macroscopic properties occur over distances, large compared to a mean free path, the concept of a local translational temperature should be valid. The rotational degrees of freedom in diatomic and polyatomic molecules also come to equilibrium at a rapid rate though not as fast as the translational mode. It is the vibrational degrees) of freedom which poses problems. Equilibration or relaxation times for the vibrational degrees of freedom are in general much larger than for the translational or rotational degrees. The reason for this is that energy differences between quantized vibrational energy levels are much greater than for the translational or rotational levels. As a result, when the temperature is being raised rapidly, the number of collisions needed to excite the next higher vibrational energy level may be significant, of the order of 5000, compared to 50 for translational and rotational equilibrium. For triatomic molecules (N20, CO2) it may be more, about 250050,000. Thus the translational and rotational temperatures of a gas being heated rapidly follow their equilibrium values closely while the vibrational temperature may still be its coldgas value. In such cases a strict unique definition of temperature is rendered impossible. It may be neces sary to define two different temperaturesone translational and the other vibrational. Thus the definition of temperature is linked to the requirement of equilibrium. If there is nonequilibrium between the different chemical species or the different degrees of freedom of a particular species, the MaxwellBoltzmann equilibrium distribution function ceases to hold, and so, strictly speaking, a temperature cannot be defined for the system. To use the concept of temperature for a nonequilibrium system, the assumption generally made is that even though the system as a whole is not at equilibrium, local equilibrium exists at any particular point in the system. This boils down to assuming that in any elementary volume (whose dimensions are 'small' compared to those of the system) the particles of the different species are in equilibrium with respect to each other and their individual degrees of freedom, so that the MaxwellBoltzmann energy distribution may be used to define a local temperature. This is how the use of temperature in flame analysis is justified. In flames where the transfer of energy between the vibrational and translational energy modes is extremely slow, the need for separate definitions of translational and vibrational temperatures may destroy the concept of a local temperature also. Such cases have been shown to exist by spectroscopic means [20]. The chemical reaction rates and the transport properties are evaluated at values of temperature based on a static system where local equilibrium has been achieved. For the case just mentioned these evaluations may be somewhat different due to nonequilibrium effects; it is not known whether the difference is significant in all cases. Some work on the high temperature oxidation of hydrogen and the dissociation of oxygen has been reported. The results show that the slowness of vibrational relaxa tion does change the processes to a certain extent; the chemical reactions occur at 'temperatures' different from their local equilibrium values. For chemical reactions which are extremely rapid (halflives of 10 microseconds or less) the classical assumption of equipartition 75 of energy is invalid. Observations for some rapid reactions have revealed that some of the reaction products are formed almost entirely in higher vibrational levels (e.g., 14th or 15th) and the transition to the equilibrium ground state occurs rather slowly [21]. REFERENCES FOR CHAPTER IV 1. Spalding, D. B., "The theory of flame phenomena with a chain reaction," Phil. Trans. Roy. Soc. London, A249, 1 (1956). 2. Williams, F. A., "Combustion Theory," AddisonWesley (1965). 3. Hirschfelder, J. 0., Curtiss, C. F., and Bird, R. B., "Molecular theory of gases and liquids," John Wiley and Sons, New York (1964). 4. Pritchard, H. 0., "The kinetics of dissociation of a diatomic gas," J. Phys. Chem., 65, 504 (1961). 5. Rice, 0. K., "On the relation between an equilibrium constant and the nonequilibrium rate constants of direct and reverse reactions," J. Phys. Chem., 65, 1972 (1961). 6. Berlad, A. L. and Yang, C. H., "On the existence of steady state flames," Combustion and Flame, 3, 447 (1959). 7. Von Karman, T. and Penner, S. S., "Selected combustion problems," AGARD, Vol. I, Butterworths, London (1954). 8. Nash, J. F., "An explicit scheme for the calculation of three dimensional turbulent boundary layers," Trans. of ASIE, Journal of Basic Eng., 94, 131 (1972). 9. Bird, R. B., Stewart, W. E., and Lightfoot, E. N., "Transport Phenomena," John Wiley and Sons, New York (1960). 10. Tricomi, F. G., "Integral Equations," Interscience Publishers, New York (1970). 11. Bledjian, L., "Conputation of timedependent laminar flame structure," Combustion and Flame, 20, 5 (1973). 12. Spalding, D. B., Stephenson, P. L., and Taylor, R. G., "A calculation procedure for the prediction of laminar flame speeds," Combustion and Flame, 17, 55 (1971). 13. Gear, C. W., "The automatic integration of ordinary differential equations," Comm. of ACH, 14, 176 (1971). 14. Seery, D. J. and Bowman, C. T., "An experimental and analytical study of methane oxidation behind shock waves," Combustion and Flame, 14, 37 (1970). 15. Baulch, D. L., Drysdale, D. D., and Lloyd, A. C., "Critical evaluation of rate data for homogeneous gas phase reactions of interest in high temperature systems," Dept. of Physical Chemistry, The University, Leeds, England 4lay 1968July 1970). 16. Bahn, G. S., "Reaction rate compilation for the HON system," Gordon and Breach Science Publishers, New York (1968). 17. Bowman, C. T., "An experimental and analytical investigation of the high temperature oxidation mechanisms of hydrocarbon fuels," Combustion Sci. Tech., 2, 161 (1970). 18. Newhall, H. K., "Kinetics of engine generated nitrogen oxide and carbon monoxide," Twelfth Symp. (Int'l) on Combustion, The Combustion Institute, Pittsburgh, Pa. (1969). 19. Ellis, G. and Bahn, G., "Engineering selection of reaction rate constants for gaseous chemical species at high temperatures," The Combustion Institute, Western States Section, Paper 6227 (1962). 20. Broida, H. P., "Temperature: Its measurement and control in science and industry," Vol. II, Reinhold Publishing Corp., New York (1955). 21. Strehlow, R. A., "Fundamentals of Combustion," International Textbook Co., Scranton, Pa. (1968). CHAPTER V MATHEMATICAL SIMULATION OF THE OTTO CYCLE 5.1 Introduction The Otto cycle may be envisioned as a sequence of thermodynamic processes, during which the thermodynamic properties and the chemical composition of the working fluid undergo changes. The cycle starts with the intake of a fuelair mixture at state 1 (Fig. 5.1). The mixture is compressed until state 2, where it is ignited with an electric spark. Process 23 represents the combustion and 34 is the expansion process at the end of which the gases are expelled from the engine cylinder and a fresh charge drawn in. The early theoretical analyses made several assumptions to simplify the problem and make it solvable without a great deal of computational work. For example, it is common to assume that processes 12 and 34 are isentropic and that the combustion process 23 and exhaust process 41 occur at constant volume. Further simplifications are obtained through the assumption that air is the working fluid throughout the cycle and its properties (specific heats) are those of atmospheric air (known as the air standard cycle). Recent analyses have removed several of these simplifying assumptions and modeled the processes occurring in the Otto cycle engine on a more realistic basis. This has largely been possible due to the availability of highspeed electronic computing machines. A review of the early 1 1 1 ?I ~`~31rrrrcnra~cI Fig. 5.1. pV diagram for the Otto cycle  r~ I' f  ~1Ur*LII ^ II__ 1__1_ work is given in Ref. 1. Among some of the recent work considering the complex sequence of processes to some detail are those of Edson [2], Spadaccini [31, and others. Edson's analysis was made on the assumption of chemical equilibrium at every step in the cycle; the main aim, however, was to study the influence of compression ratio and dissociation on the thermal efficiency. This assumption of chemical equilibrium has been made in several analyses due to the lack of chemical kinetics data. Spadaccini's analysis [3] has been one of the most detailed analyses made so far. In this analysis the compression process is assumed to be isentropic and the specific heat is assumed to be a linear function of temperature for this process. Combustion is assumed to occur at constant volume and the composition and thermodynamic properties of the mixture after combustion are found on the basis of a chemical equilibrium analysis. For the expansion process, a nonequilibrium finite rate kinetics analysis with allowances for heat transfer is made considering a total of 13 chemical species to be present. It will be seen later that an equilibrium analysis of the expansion process does not convey the true state of affairs with respect to the concentrations of the species at exhaust. Since the present interest is to study the exhaust emissions, a nonequilibrium analysis is a must to obtain a realistic estimate of the emission levels. Other work in this area has been referred to at the appropriate places in this chapter. 5.2 Present Work Here again the compression process has been assumed to be isentropic. The combustion process is initiated sometime before top dead center is reached and completed shortly after top dead center. Instead of assuming the combustion process to be instantaneous as done in many of the previous analyses, here the process has been assumed to proceed at a finite rate. Using an assumed burning rate law which gives the mass fraction of the total fuelair mixture burnt as a function of the crank angle, a stepbystep analysis of the combustion process has been made to obtain the state properties and the chemical composition of the working fluid at the point where combustion is complete. After this, for the expansion process analysis, the same reaction mechanism as used by Spadaccini has been used here; however, different mathematical techniques have been used to solve the governing equations. The concentrations of the various species at exhaust are thus obtained and so this model can be used to study the influence of several parameters such as airfuel ratio, exhaust gas recirculation, compression ratio, speed, ignition timing, etc., on the exhaust emissions. Due to the approximations involved in the theoretical analysis, and the fact that the processes occurring in the exhaust pipe have not been considered in this analysis, agreement with the experimentally measured values may not be exact. However, the trend in the results will definitely show up in this theoretical analysis and so experimental runs need be made to obtain results only for a few combinations of the different parameters which seem promising from the emissions viewpoint, after a preliminary detailed theoretical analysis. The only empirical input here is the burning rate law. The validity of the burning rate law can be checked by comparing the calculated thermodynamic properties (i.e., the pressure and the temperature) with the experimentally measured values. The law may then be modified, if necessary, to obtain better agreement with experimental results. The assumptions made here agree qualitatively with the experimental results of Refs. 4 and 5 and other published literature; quantitative agreement of the mass fraction burnt as a function of the crank angle could not be checked due to differences in the operating conditions. A rigorous theoretical analysis involving the solution of the general threedimensional timedependent conservation equations for a flame propagating into a premixed fuelair mixture would yield the mass fraction burnt as a function of the crank angle. However, with present day computational techniques, this type of an analysis of the problem is almost impossible. Moreover, the analysis is beset by difficulties due to the limited knowledge of the complex nature of the processes occurring in the engine: the reaction mechanism and the rate constant values, degree of turbulence in the combustion chamber, problems arising from the fuelair mixture not being homo geneous and the fuel not being completely vaporized, heat transfer at the combustion chamber walls and the associated quenching effects on the flame, together with a host of other difficulties. The effect of the injection of alcoholwater mixtures on the emissions has been studied using the model. In general, many of the control techniques which have been suggested involve some treatment of the exhaust gases. The injection of wateralcohol mixtures strives to reduce the formation of the pollutants in the engine itself rather than offer some treatment of the exhaust gases. The hydrocarbon emissions result mainly from the heterogeneous processes occurring in the quench zone and are influenced by several other factors. Due to the lack of knowledge of the complex processes occurring, it was not possible to predict the hydrocarbon emissions using this model. However, it should be remembered that control techniques which result in the reduction of carbon monoxide also generally result in the reduction of hydrocarbons. So the latter need not be studied separately, but can be determined by subsequent experiments after the preliminary theoretical analysis. 5.3 Major Assumptions for Analysis Before describing the mathematical model of the Otto cycle, a discussion of the assumptions made in the analysis is in order. In this section the assumptions relating to the entire cycle are mentioned. Assumptions for the specific processes are mentioned in the appropriate places in this presentation. (a) The fuel is assumed to be completely vaporized and mixed homo geneously with air. If the fuelair mixture is not homogeneous, it will be locally fuelrich at some points and locally fuellean at others, though it may be overall stoichiometric. Compared to stoichiometric combustion, the fuelrich parts would yield larger amounts of unburnt hydrocarbons and carbon monoxide while the fuellean parts would yield a larger amount of the oxides of nitrogen. Also the energy release would be less since a lesser amount of the fuel is actually burnt. The result would be (as compared to the homogeneously mixed case): (i) relatively low mean exhaust temperature, and (ii) higher concentrations of carbon monoxide, unburnt hydrocarbons, and oxides of nitrogen. Similar conclusions would hold when the mixture was overall fuellean or fuelrich, and not homogeneous. Thus the theoretical calculations assuming perfect mixing represent the upper limit of performance and the lower limit of the emissions. (b) The gas mixture is assumed to be thermally perfect. Nitrogen is a major constituent (> 70%) of the mixture and it closely approximates thermally perfect behavior for the conditions encountered in a typical Otto cycle engine. Other constituents are present in relatively smaller concentrations, so as not to affect the overall thermally perfect behavior. The only constituent likely to be considerably thermally imperfect would be a fuel like isooctane. However, since its concentration in the mixture is small (mole fraction < 0.02), the results would be affected only slightly. (c) Reactions in the unburnt gas region are neglected, i.e., there are no precombustion reactions. The maximum temperature in the unburnt gas region is generally less than 1000 degrees K for normal combustion, so this is a fairly good assumption. Some reactions (i.e., low temperature pyrolysis, slow oxidation, cool flames) do occur in this region in practice, but the rates of these reactions are significant only at temperatures above about 1400 degrees K. (d) The quench layer is assumed to occupy a negligible fraction of the total volume. The quench layer is a major source of the hydro carbons and products of incomplete combustion present in the engine exhaust. This assumption implies that only the homogeneous chemical and physical processes occurring in the bulk of the gases are considered and the heterogeneous processes occurring near the combustion chamber walls are neglected. (e) There is no heat flow or mixing between adjacent elements of gas during combustion, i.e., transport phenomena of viscosity, thermal conductivity (also radiation and convection) are neglected. Transport phenomena are of importance only in the flame element where appreciable gradients of temperature and species concentrations exist, and these have been considered in the analysis of the flame zone presented earlier in this study. In the unburnt gas region, gradients in these properties are generally absent, since this region is fairly homogeneous. In the burnt gas region, gradients in these properties are present, but are not appreciable enough to cause significant fluxes of heat and mass. Besides, in this study a bulk analysis has been made considering the entire burnt gas region to be at a uniform temperature and concentration. (f) No pressure gradients exist in the combustion chamber. Although the pressure changes with time, it is spatially uniform within the combustion chamber; any pressure gradient introduced would be 
Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID ECBB9YIUD_NM3LSR INGEST_TIME 20120308T21:45:23Z PACKAGE AA00003949_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 