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 
COMPUTATION OF HIGHSPEED REACTING FLOWS By JAMES KEITH CLUTTER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1997 To my wife, Paula ACKNOWLEDGEMENT I would like to express my gratitude to those individuals who have helped me accomplish this goal. Dr. Bruce Simpson, Russ Klugg and everyone at Eglin AFB made it possible for me to start this effort while on active duty. Also, the continued financial support from the Armament Directorate of Wright Laboratory, Eglin AFB, has been essential as well as the computer resources they have put at my disposal. The collaboration and interaction with Gregg Abate and Rudy Johnson have both helped steer my research and provided key technical input. The guidance I have received from Dr. Dave Belk has been essential and I can not thank him enough for his encouragement and for serving on my committee. I would like to thank Dr. Hartmuth Lehr for providing his experimental data and for his inputs on the topic of shockinduced combustion. Also, Dr. Greg Wilson has generously provided both reference material and insights that are greatly appreciated. Dr. Mark Janus has been a good friend and provided both encouragement and education into the finer arts of motorsports and statistics. Several individuals at the University of Florida have helped me see this effort through. First and foremost is Dr. Wei Shyy who has been both a mentor and friend. Completion of this effort would not have been possible without his guidance, encouragement, and tolerance. I thank Dr. Ed Milton for helping me begin this endeavor and Dr. Renwei Mei, Dr. Corin Segal, and Dr. James Klausner for serving on my committee as well as Dr. David Mikolaitis for his inputs and suggestions in the area of kinetics. I would also like to thank my colleague, Dr. Venkata Krishnamurty, who has provided valuable insight into the modeling of turbulence. My family has provided endless support whether it be emotional or dropping in to help with the kids during times of peak studying and work. I can not thank my mom and my motherinlaw enough as well as my fatherinlaw, Richard Cooper, who recently passed away. I must thank my father, the eternal student and consummate mechanic, for instilling in me a desire to learn how things work and the moxie to take them apart. My children, Courtney and Cooper, are without a doubt the finest children in the world and have been immensely understanding and patient. They have also proven to be delightful distractions and have helped keep me focused on what is truly important. Finally, there is my wife, Paula, whose strength and encouragement are immeasurable. She has never ceased to support me whether it meant her understanding in me quitting my job to return to school or balancing her career and graduate school along with mine. I will never be able to express all she means to me. TABLE OF CONTENTS ACKNOWLEDGEMENT ............................................ iii ABSTRACT ....................................................... vii CHAPTERS 1 INTRODUCTION ............................... ............ 1 1.1 Reactive M unition ................................. ......... 1 1.2 Gun Muzzle Blast .......................................... 6 1.3 Phenomenological Scales .................................... 11 1.4 Summary ....................................... ........ 14 2 SHOCKINDUCED COMBUSTION ............................. 17 2.1 Reaction Mechanisms for Hydrogen and Oxygen .................. 21 2.2 Previous Computational Studies of Lehr's Experiments ............. 22 2.3 Summ ary ................... .................... ........ 33 3 GOVERNING EQUATIONS AND NUMERICAL TECHNIQUES ...... 37 3.1 Governing Equations ....................... ............ .... 37 3.2 Equation of State and Temperature ............................ 41 3.3 Thermal Conductivity and Viscosity ................... ....... 43 3.4 Flux Vector Splitting Scheme ........................ ....... 44 3.5 Chemical Reaction Source Term ............................... 47 3.6 Special Treatment For Chemical Length Scales ................... 49 3.6.1 Transformed Governing Equations ........................ 50 3.6.2 Logarithmic Interpolation ............................... 51 3.6.3 Error Estimation Based on OneDimensional Model Equations .. 53 3.6.3.1 Global Order of Accuracy Analysis ................. 57 3.6.3.2 Current Application of Scaling Factor Scheme ......... 61 3.6.4 Nonphysical Wave Speeds .............................. 63 3.7 Numerical Solution Methods ................................ 64 3.7.1 FullyImplicit Formulation ................ .............. 64 3.7.2 ExplicitImplicit Splitting Scheme ........................ 67 3.7.3 ODE Solution Issues ................................... 69 3.8 Summary .............................................. 71 4 TURBULENCE MODEL ....................................... 72 4.1 Averaging Techniques and Averaged Variables .................... 75 4.2 Averaged Equations ........................................ 77 4.3 Turbulence Closure ................................ ..... .. 79 4.4 ke Model ................................................ 81 4.5 Wall Function ............................................. 83 4.6 Turbulence Effects On The Reaction Source Term ................. 85 4.7 Summary ................................................ 86 5 SHOCKINDUCED COMBUSTION RESULTS AND DISCUSSIONS ... 88 5.1 Comparison of Solution Methods .................... ........ 94 5.2 Source Term Treatment ...................................... 96 5.3 Effect of Reaction Mechanism on the Computed Flow Field ......... 104 5.3.1 Pressure Effects on Reaction Rates ........................ 123 5.3.2 Summary of Reaction Mechanism Effects .................. 128 5.4 Computations Using Turbulence Model ......................... 129 5.5 Summary ...................................... ......... 136 6 GUN BLAST AND FLASH LITERATURE REVIEW ................ 138 6.1 Gun Muzzle Blast ........................................... 138 6.2 Gun Muzzle Flash ......................................... 142 6.3 Summary ................................................. 143 7 GUN MUZZLE BLAST RESULTS AND DISCUSSION .............. 145 7.1 Large Caliber Howitzer Simulations ............................ 145 7.1.1 Flow Field Model and Boundary Conditions ................ 146 7.1.2 Experimental Overview ................................ 150 7.1.3 Computational Results ................................. 151 7.2 Small Caliber Sound Suppressor Simulations .................... .160 7.3 Summary ................................................ 165 8 CONCLUSIONS AND FUTURE WORK ......................... 168 8.1 ShockInduced Combustion .................................. 168 8.2 Gun M uzzle Blast ................................ ...... 171 8.3 Summary ................................................ 172 APPENDIX JACOBIANS AND EIGENVECTORS ...................... 173 REFERENCES ................................................... 177 BIOGRAPHICAL SKETCH ........ ............................. 185 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATION OF HIGHSPEED REACTING FLOWS By JAMES KEITH CLUTTER May 1997 Chairperson: Dr. Wei Shyy Major Department: Aerospace Engineering, Mechanics and Engineering Science A computational study has been conducted for highspeed reacting flows relevant to munition problems, including shockinduced combustion and gun muzzle blast. The theoretical model considers inviscid and viscous flows, multispecies, finite rate chemical reaction schemes, and turbulence. Both the physical and numerical aspects are investigated to determine their impact on simulation accuracy. A range of hydrogen and oxygen reaction mechanisms are evaluated for the shockinduced combustion flow scenario. Characteristics of the mechanisms such as the induction time, heat release rate, and second explosion limit are found to impact the accuracy of the computation. On the numerical side, reaction source term treatments, including logarithmic weighting and scaling modifications, are investigated to determine their effectiveness in addressing numerical errors caused by disparate length scales between chemical reactions and fluid dynamics. It is demonstrated that these techniques can enhance solution accuracy. Computations of shockinduced combustion have also been performed using a kE model to account for the turbulent transport of species and heat. An algebraic model of the temperature fluctuations has been used to estimate the impact of the turbulent effect on the chemical reaction source terms. The turbulence effects when represented with the current models are found to be minimal in the shockinduced combustion flow investigated in the present work. For the gun system simulations, computations for both a large caliber howitzer and small caliber firearms are carried out. A reduced kinetic scheme and an algebraic turbulence model are employed. The present approach, which accounts for the chemical reaction aspects of the gun muzzle blast problem, is found to improve the prediction of peak overpressures and can capture the effects produced by small caliber firearm sound suppressors. The present study has established the numerical and physical requirements for simulating highspeed reacting flows. Also, key parameters useful in quantifying the roles of these aspects have been assessed. CHAPTER 1 INTRODUCTION In the present work, a computational capability is developed to analyze high speed reacting flows arising in munition problems. The particular problems focused on are a munition concept and the occurrence of gun muzzle blast. These problems are difficult to analyze given the fact many physical and chemical processes are involved. First a description of each problem is presented to highlight these various processes and frame the motivation of the current research. During these descriptions the phenomenological issues such as the fluid dynamics, chemical reactions, and turbulence are highlighted. This is followed by a description of how the time and length scales of these various aspects impact computational models. Finally, the specific scope of the present thesis is defined and the approach adopted to solve these problems is presented. 1.1 Reactive Munition The first munition problem of interest is a weapon concept dependent on the mixture of fuel and oxidizer and shock initiation and is similar to the concepts discussed by Gore et al. (1993). This class of weapons differs from conventional munitions in the manner by which they derive their destructive force. This destructive force has to be generated by some energy source and in typical munitions, such as general purpose bombs, the source is a condensed phase explosive. This explosive contains the chemical compounds needed to generate the force once initiated by a fuse. Alternate munition concepts are based on the idea of generating the destructive force from alternative means. This may be done to enhance the energy source beyond that which is carried to the target but another focus is to use materials on target as energy sources. This is highly attractive for munitions to be launched from both ships and aircraft given the fact that the less material that must be carried to the target, the lighter and smaller the munitions system must be. This directly impacts logistic concerns such as storage as well as factors such as the fuel required to deliver the munitions. This current study is part of an effort to study the various physical and chemical phenomena which directly impact the effectiveness of new munition concepts. Another facet of the effort is to quantify the modeling tools needed to simulate this class of weapons and to aid in the development of engineering codes useful in design and evaluation of these various concepts. The range of concepts and targets which fall in this category is quite extensive and a complete review is beyond the scope of this document. Therefore, one target and one munition concept are highlighted here. However, it should be noted that the underlying physical and chemical processes to be presented in the current context are common to all targets and munition concepts. As mentioned, one of the premises of various new munitions is to take advantage of materials on target which can be compromised and used to defeat the target. One component common to all air vehicles whether missiles or aircraft is onboard fuel. Given the particular vehicle, the fuel can be stored in various states. For brevity, here we consider an air vehicle carrying fuel in a fluid state. For example, most aircraft contain numerous fuel cells in the wings and fuselages. Depending on the stage in the mission, the fuel cells contain some level of liquid fuel as well as fuel vapor. Therefore, the munition concept used to help define the important physical and chemical phenomena to be studied centers around using the onboard fuel in a combustion process as the source of energy to generate the destructive force. Even though there can be an ample supply of fuel on target and an ample supply of oxidizer, the ambient air, if the required conditions are not achieved, the munitiontarget system can remain in a fuel or oxidizer rich state and no effective combustion is initiated. One plausible system concept to attack the fuel cells is a weapon that generates a "cloud" of objects designed to penetrate and destroy the fuel cells. The focus of the current project is to study various key phenomena to aid in the development of an engineering level model with which to quantify the destruction output once this "cloud" intercepts the fuel cell. Continuing with the concept of attacking the fuel cell, the target can be thought of as a control volume bounded by a solid wall with specific characteristics such as material and thickness. In the context of the larger issue of weapon effectiveness, criteria which need to be quantified and are dependent on output from any model of the munitiontarget interaction are factors such as whether the walls are damaged such that this damage is rapidly transferred throughout the entire air vehicle. Two tragic examples of such a process are the recent Valujet and TWA Flight 800 accidents. In the first case, the predominant theory is that the damage was initiated by the failure of one or several oxygen generators (McKenna, 1996). In the TWA case, as of this writing, it is still unknown whether the damage to the aircraft had been initiated by an explosive device or a system malfunction. Regardless of the initiation process, damage in both cases had been tremendous due to the fact that conditions were such that the destruction process was sustainable. In any combustion process, which both these cases entailed, sustainment is dependent on the correct supply of fuel and oxidizer as well as the required system settings. The first goal of this effort is to better understand the relationship between such system settings and the physical and chemical phenomena which conspire together to generate global effects such as the destructive force output. The second goal of the present work is to quantify the attributes of the computational model required to simulate these phenomena. Our fuel cell scenario is depicted in Figures 1.1 and 1.2 and the steps leading to destruction are as follows. The fuel cell contains some amount of liquid fuel (F1) as well as fuel vapor (Fv). At stage 1, the projectile intercepts the fuel cell with a shock front (Sf) due to the projectile flight through the ambient region which is composed of air, an oxidizer (0). After the projectile enters the fuel cell, the structure of the established shock (Sf) changes somewhat due to the properties of the fuel. There may also be a precursor compression (Sp) wave initiated by the impact of the projectile with the fuel cell wall. The strengths of the established shock in stage 2 and precursor compression wave depend on the characteristics of the wall. If the wall is relatively thin, then the momentum transfer from the projectile to the wall is negligible before the projectile penetrates. On the other hand, for a stronger wall, the momentum transfer will be larger and the velocity of the projectile will decrease substantially during penetration. This would result in a stronger precursor compression wave and dramatically weaken the shock already established by the projectile. Fuel Cell Projectile Cloud Figure 1.1. Fuel cell scenario showing projectile cloud. Single projectile highlighted in Figure 1.2. After penetration, as the projectile and shock proceed through the cell, the heat from the shock can vaporizes any liquid fuel. This fuel vapor mixes with the oxidizer (0) which has penetrated into the fuel cell and the mixing is enhanced by the wake of the projectile. The projectile may be constructed out of an ablative material formulated to enhance the reaction process. In that case, the projectile changes phase from a solid (Ms) to eventually a gas (Mg). This compound reacts with the fuel vapor (Fv) and oxidizer (0) to produce a final product (P). The final step leading to destruction is for the reaction to reach a state such that it will progress through the fuel cell and combine with the same process initiated by the other projectiles. The overall effectiveness of this or any weapon concept is a function of the target and weapon system parameters. The target parameters here include primarily the fuel chemical composition, phase, and state which directly impact key aspects such as how much total energy is needed to compromise the fuel volume. These aspects are controlled by a combination of detailed phenomena such as chemical time scales. The second of the target I 0i I aSf I I FI (1) Projectile intercepts fuel cell. Shock (Sf) due to flight through ambient air. Fuel cell contains some level of liquid fuel (FI) and fuel vapor (Fv) Fv O P SP F1 i L  ** jI I FCI (2) Precursor (Sp) generated during penetration Heat release behind shock (Sf) produces fuel vapor (Fv) Oxidizer (0) penetrates cell IF, (3) Solid projectile (Ms) ablates (4) Reaction front (Rf) proceeds producing gas (Mg) through cell Fuel vapor (Fv), catalysis (Mg) and oxidizer (0) mix and react producing product (P) Figure 1.2. Sequence of events surrounding an individual projectile. parameters is the container characteristics such as wall material and thickness. In this example, these parameters enter into the analysis process somewhere between the modeling of the munitiontarget interaction assessment and a larger lethality model. Essentially, some sort of effective destruction power is generated and the effectiveness of this power in damaging the fuel cell wall and surrounding structures has to be quantified and this damage output is fed into a higher level lethality model. The second class of parameters is those associated with the munition. The particular munition concept used in the current discussion is a scenario where the fuel cell is infiltrated by a certain number of projectiles. Therefore, the munition parameters are first the number of projectiles and their dispersion pattern. Next are the parameters of the individual projectiles such as their size, shape, material composition and velocity. The importance of these parameters on the system effectiveness is manifested in how they drive the important detailed physical and chemical phenomena. From an application context, the quantities of interest are the global effects such as the destructive force generated during the munitiontarget interaction. In the fuel cell scenario one key global effect is whether combustion is achieved and sustained. The correlation between system parameters and these global effects is controlled by the detailed physical and chemical phenomena. It is clear from the description of the fuel cell problem that among the key physical phenomena is the mixing of the fuel and oxidizer. Given the flow scenario, turbulence plays a key role in the mixing as well as the transport of heat. Among the processes of interest are phase change and the reactions between the fuel and oxidizer. The modeling of these various phenomena is the focus of the present work. 1.2 Gun Muzzle Blast The second munition problem used in the current study is the simulation of gun muzzle blast. This problem contains many of the same physical and chemical processes discussed in the description of the reactive munition problem and requires much of the same modeling tools. However, here computational efficiency is paramount given the desire to develop an engineering tool useful in the design of muzzle devices. The primary concern in this problem is the reduction of the pressure levels generated when a conventional gun is fired. In the case of large caliber systems, the reduction is sought to alleviate loads the support vehicle must sustain. In small caliber guns, the concern is the associated hearing loss from repeated use and the ability to go undetected for certain applications. This problem has been studied by several investigators. The primary tool of investigation has been experiments with a limited computational effort. A literature review is provided in a later chapter. Here the problem is described to highlight the need for further research. The ballistic process determining blast is depicted in Figure 1.3 and can be divided into two domains. The first is the interior ballistic cycle where the solid propellant in the rounds is ignited and burned. This results in a buildup of chamber pressure which accelerates the projectile down the barrel. Following the projectile is a volume of gun gas which is created from the spent propellant. This gas is characterized by both high pressures and temperatures as well as some level of particle loading. The point in the barrel where the particle loading is not changing and the solid propellant has been combusted to essentially an equilibrium state can be considered as the boundary between the interior ballistic cycle and the launch cycle (point A in Figure 1.3). This is the point where the current code imposes specific boundary conditions. The launch cycle encompasses various physical occurrences such as the projectile exiting the barrel and the progression of the precursor and main blast wave from the barrel into the ambient environment. The precursor blast wave denotes the shock created by the flyout of the projectile while the main blast wave denotes the high pressure gun gas which exits the muzzle with the projectile. The level of blast experienced is also directly related to parameters such as the particular propellant used in the rounds. Another physical phenomenon which often occurs when a conventional gun is fired is muzzle flash. The occurrence of muzzle flash is a result of further combustion of the gun gas with the ambient air. This chemical reaction process not only generates flash but also Solid Propellant Projectile Large Grains S I..... V (1) (3) mall Particles Gun Gas / Air Interface r Precursor Blast Wave Main Blast Wave ./ Precursor Blast Wave (2) (4) 1 Figure 1.3. Interior ballistic and launch cycle of a typical conventional gun. alters the overpressure signature. Even though the focus here is on the overpressure, it is worthwhile describing the mechanics of the occurrence of flash to obtain a full appreciation of the problem. The gun gas which is discharged during the launch cycle is always high in temperature, always contains some level of solid particle loading, and is typically rich in hydrocarbons which can react with the oxygen in the ambient air. Flash associated with conventional guns can be classified in three main regions; primary, intermediate, and secondary flash. These regions not only occur at different stages in time after the gun firing, but also in different spatial regions of the field downstream of the gun muzzle. These regions are represented in Figure 1.4. The primary flash is located at the exit of the muzzle and encompasses a small spatial region and is low in intensity. However, the intermediate flash is a more extensive region of greater intensity and occurs just downstream of the shock disk which forms some distance from the muzzle exit. The location of the intermediate flash region is a direct result of the fact that it is the heating of the gun gas by the shock disk which initiates the combustion producing the flash. Another influence of the shock disk is the 1 AZ I AA reduction in the velocity of the gas which increases the residence time of the gas and promotes the initiation of the chemical reactions. In addition to the heating from the shock disk and the reduction in the gun gas velocity as it crosses the shock, a certain level of oxygen must be mixed with the gun gas before the flash occurs. This is why there is typically a delay between the firing of the gun and the occurrence of the intermediate flash. This delay is directly dependent on the time it takes to entrain oxygen from the ambient air and to mix it through the volume of gun gas toward the line of fire. The entrainment process is highly dependent on the turbulent processes which operate along the interface between the gun gas and ambient air. Therefore, any alteration to the gun system such as the addition of a muzzle device will impact the mixing process and in turn the intermediate flash. The level of flash generated from the intermediate region is typically repeatable from firing to firing for a given gun system, i.e. gun caliber, ammunition, and muzzle device. The caliber of the gun and the ammunition fired determine the volume of the gas discharged as well as the gas properties when it exits the muzzle. These gas properties set the characteristics of the shock disk as well as the conditions upstream of the disk such as the fuel richness. Along with the influence of any muzzle device mentioned earlier, these parameters determine the level of flash in the intermediate region. Secondary flash is the final and typically the largest and most intense region of the three. Like the intermediate flash, the secondary flash is initiated by the mixing of the oxygen rich air and gun gas and also requires a source of energy to initiate the reaction. And as before, the residence time of the mixture must be such to allow combustion to initiate and progress. In the case of the intermediate flash, the shock disk provides the energy source as well as reduced the velocity of the gas to ensure a high enough residence time. In the case of secondary flash, it has been found (Klingenberg and Heimerl, 1992) that if the temperature in the intermediate flash region does not reach a certain level, dependent on the ammunition, then secondary flash does not occur. This supports the theory that the heat released in the Primary Barrel Gas Intermediate Flash Primary Flash Interface of Gun and Air Main Blast Wave Figure 1.4. Depiction of flash regions and the development of muzzle blast and flash. ,I < intermediate flash is the energy source for initiating the secondary flash which occurs some time after the intermediate flash and further downstream. The residence time criterion is met by the ever increasing volume of gun gas. As the cloud of gas grows larger, the velocity throughout the cloud is reduced. Therefore, as depicted in Figure 1.4, the secondary flash covers a larger area of the field. Secondary flash from shot to shot is not as repeatable as intermediate flash. This is due to the difference in the initiation process. In the case of intermediate flash, the initiation is driven by the shock disk which is less influenced by small perturbations in both the ambient field and the particular round fired. The initiation of secondary flash is much more dependent on the transport of both species and heat in a much larger region, further from the gun muzzle. This makes it more accessible to ambient influences such as winds. To simulate all the intricacies of flash requires a detailed model of the chemical composition of the gun gas. Some modeling efforts have been carried out and these are reviewed later. Since the focus of the current study is to predict the overpressure in regions somewhat removed from the gun system, the key phenomenon associated with the reactions which cause flash that must be reproduced is the energy released from the combustion. This is due to the fact that the combustion occurs behind the main blast wave. This deposition of energy behind this wave helps sustain the magnitude of the main blast wave as it progresses out into the field. This is key since the peak pressure at any location in the field is a function of the main blast wave magnitude. Therefore, the current model employs a reduced reaction mechanism in an attempt to better reproduce the blast levels without greatly increasing the computational requirements. However, the characteristics of the flash processes are to prove advantageous in evaluating the model. 1.3 Phenomenological Scales The description of the problems above makes it clear that there are a multitude of physical and chemical processes key in highspeed, chemically reacting flows. Associated with this myriad of mechanisms are multiple time and length scales. The time and length scales refer to the characteristic time and length over which effects due to certain physical and chemical mechanisms occur. Consequently, it is difficult to investigate this class of problems both theoretically and experimentally. Therefore, computational models are routinely used to study reacting flows and this is the methodology of the current work. Any models used to simulate the various detailed phenomena involved must correctly reproduce all phenomena and their corresponding scales, or at the least the role the phenomena play in the larger problem. As for the numerical aspects of such models, the requirement to resolve a wide range of scales determines simulation parameters such as the computational time and computer memory needed to carry out a useful simulation. This is due to the fact that when solving a problem using computational models, the governing equations are usually solved at discrete points in time and space. Therefore, the temporal and spatial discretization must be done with both the time and length scales of interest in mind. A key parameter of chemically reacting flows is the flow speed or the fluid dynamic time scale relative to the chemical processes. A nondimensional parameter useful in classifying problems is the Damk6hler number defined as the ratio between the fluid dynamic and chemical time scales, ) = Tf/Tc. For lowspeed flows, the chemical processes typically occur orders of magnitude faster than the fluid dynamic processes. This implies ) > 1 and adds validity to modeling assumptions such as local chemical equilibrium or thin flame zones. From a modeling standpoint, such assumptions allow the use of a simplified chemical reaction model. For instance, as a result of the equilibrium assumption, the finite rate chemical processes can be ignored and the presence of intermediate species may be neglected. Even if the flow is not in chemical equilibrium but the reactions can be considered fast, either a partial equilibrium (Correa and Shyy, 1987) or a progress variable (Libby and Williams, 1994) concept can be adopted. For the case of highspeed flows as in the case of shockinduced combustion, such assumptions are generally not valid and a model of the finite rate chemical reactions accounting for the detailed distribution of intermediate species is needed. This has direct implications as to the complexity of the computational model. Many combustion problems have to address the issue of turbulence. This is due to the fact that the mixing characteristics of turbulence are routinely present with high Reynolds number flows. Also, in many cases turbulence is used in reactive systems to enhance the mixing and subsequently the combustion process. As in all turbulent flows, the aspects due to turbulence must first be compared to the mean flow and with the addition of reactions, the relative scales of the turbulence to the chemical processes are important. Turbulence is known to possess a multitude of scales, particularly the large scale eddies which contain what is labeled the turbulent kinetic energy (k) and the small scale eddies which work to dissipate this energy. The dynamics of turbulence consist of various processes by which the large eddies break down into small eddies and thereby transfer energy between the large and small eddies. Depending on the relationship between the large and small scales of turbulence and those of the chemical reaction some conclusions can be drawn as to the impact of turbulence on the combustion process. A parameter useful in the analysis process is the ratio of the largeeddy time scale to the chemical time scale, the largeeddy Damk6hler number (I) = tj/Tc). The largeeddy time scale, t, = I/ k, is defined using the characteristic length scale (1) and velocity (k) of the eddy. Using this parameter, the regime of interest runs from slow chemistry, small dt, to fast chemistry, large D/. For the case of slow chemistry, the turbulent processes are working on a scale smaller than the chemical process. This produces enhanced mixing of reactants or even intermediate species resulting in the reaction zone being distributed in the flow field. In this case, the impact of turbulence on mass and thermal diffusion must be considered. For the case of fast chemistry the reactions take place on a scale much smaller than the large eddies. This results in distinct thin flame fronts which lie within the large eddies and are transported and strained by the turbulence. Such cases can be modeled using the laminar flamelet concept which is reviewed by Libby and Williams (1994). For the highspeed flows to be studied here, it is probable that the turbulent and chemical time scales are comparable. Therefore, one goal of the current study is to incorporate the effects of turbulence on the chemical processes in the computation. 1.4 Summary Here the problems to be studied have been described to identify the important physical and chemical processes which must be modeled. Given the highspeed nature of the flows, many assumptions used to reduce the complexity of models employed in lowspeed combustion are not applicable. Some of the length and time scales have been identified and how these impact the numerical aspects of the model is discussed later. The first primary focus of the current effort is to quantify the complexity of the models needed to accurately simulate the munition problems of interest. The second focus is to identify the impact both the phenomenological and numerical aspects of these models have on the accuracy of these simulations. The remainder of this document describes the specific models and computations to be carried out. Chapter 2 describes in detail a shockinduced combustion scenario to be used as a representation of the munition problem consisting of the fuel cell scenario. This is done in part due to the availability of experimental data associated with the particular scenario. Also, this scenario has been studied by previous investigators using a wide variety of physical, chemical, and numerical models. A literature review of these earlier works is given to help draw clear comparisons between these efforts and the current study. Using the code developed here, many of these various models can be tested within a single computational framework to ascertain the correlation between these various models and simulation accuracy. This not only allows judgments to be made as to the validity of the various models but also aids in quantifying the fidelity needed in a computational tool for the munition concepts of interest. The governing equations of the current class of munition problems are presented in Chapter 3 along with a detailed description of the computational code developed here. It is to become evident that there exist numerous modeling options for the various facets of the code from both a phenomenological and numerical point of view. Where appropriate these options are discussed. The assumptions which drive the model selection are detailed as well as compelling reasons for the current model selection based on both the physical and chemical aspects of the problems to be solved and the computational requirements the solution process must meet. Also discussed in Chapter 3 are the various numerical techniques used in the present computational code. Chapter 4 describes the turbulence model used in the current effort as well as the various terms which are impacted by the turbulence aspects of the flow scenarios. The model selection for the representation of turbulence is driven in a large part by the computational requirements of the code developed here. Various model modifications for addressing attributes of the highspeed flows to be solved are also discussed. Results for the shockinduced combustion scenario are presented in Chapter 5. Among the data reviewed is analysis as to how the chemical model and numerical aspects of the code impact accuracy. Computations using both inviscid and turbulence fluid dynamic models are presented to both test the current code but also the ascertain to what extent the turbulence aspects of the flow impact the combustion scenario. The current results are reviewed in the context of previous findings both to support the conclusion of the present effort and to fill in those aspects of solving this class of flows not directly addressed in earlier studies. A literature review of pertinent work in the area of gun muzzle blast is provided in Chapter 6. Many of the conclusions and findings of these earlier efforts have helped steer the development of the current code used to simulate this munition problem. It is to be made evident that the models to date have ignored some of the physical and chemical processes and it is the focus of the present effort to quantify the impact these processes have on the generation of muzzle blast. This is followed in Chapter 7 by a description of how the model discussed in Chapter 3 is used to simulate the muzzle blast scenario. Chapter 7 also presents results using a computational model analogous to those used in earlier efforts and compares these simulations to results from the current model. Finally, Chapter 8 summarizes the findings of the present study into both the shockinduced combustion and gun muzzle blast problems. The accuracy of the code developed here is discussed as well as the computational requirements. Areas into which more focus is needed are identified in the context of the munition problems of interest. Also, in a more general sense, those issues in need of more attention to increase the capability of simulating highspeed reacting flows for full scale engineering problems are identified and some suggestions for the future direction are made. CHAPTER 2 SHOCKINDUCED COMBUSTION The fuel tank scenario described earlier contains a multitude of physical and chemical processes. The predominate processes are the shock dynamics, mixing, phase change, and chemical reactions between the fuel and oxidizer. For the current look into this class of problems, a premixed combustion scenario is studied. The shock dynamics are explicitly present and the mixing aspects are manifested in the importance of the distribution of intermediate species and heat. The chemical processes are represented by the chemical reactions between the fuel and oxidizer. Aspects also important to the phase change process are present. These include the competition between the fluid dynamics and chemical time scales and the dependence of the process on local field properties. The combustion scenario to be used in the present research for studying the reactive munition problem is the shockinduced combustion scenario studied experimentally by Lehr (1972). In his experiments, Lehr shot projectiles into stoichiometric hydrogen and air mixtures as well as hydrogen and oxygen mixtures. Projectiles of various geometries have been used but here only data from shots using hemispherically nosed projectiles are studied. Lehr's primary findings are that the flow field features generated by the shockinduced combustion can either be one of steady combustion or a situation in which unsteady combustion occurs. The case of unsteady combustion can be further subdivided into low and high frequency regimes. The regime achieved is found to be directly related to the ratio of the projectile velocity to the detonation velocity of the mixture. The data available from Lehr's study are qualitative in nature and are in the form of shadowgraphs. A representative image from the three regimes has been provided by Lehr (personal communication, 1995) and is shown in Figure 2.1. The shock in all three images is visible as is a second front due to the chemical reactions and labeled the energy release, or reaction front. All three of these cases are for a stoichiometric mixture of hydrogen and air at a temperature of 293 K and a pressure of 320 Torr. As is clear by the images in the figure, for those cases where the velocity of the projectile is lower than that of the detonation velocity of the mixture (VD), an instability is seen in the reaction front. As the velocity approaches the detonation velocity, the frequency of the instability increases until the combustion becomes stable at velocities greater than VD. It should be noted that at the other end of the velocity spectrum, for low Mach numbers, no Vp = 1685 m/s Vp/VD = .82 M= 4.18 freq = 150 kHz Vp = 1931 m/s Vp/VD = .94 M = 4.79 freq = 720 kHz Vp = 2605 m/s VPVD = 1.27 M = 6.46 freq = steady Figure 2.1. Example shadowgraphs from Lehr's Experiments. combustion is initiated. Several investigators have proposed theories as to the important phenomena driving the instability and have supported their theories with computations. Among these is the work of McVey and Toong (1971), Alpert and Toong (1971), Matsuo et al. (1993), and Matsuo and Fujiwara (1993). There is also some question among the investigators as to when the actual transition from unsteady to steady occurs. Sussman (1994) presents arguments that the transition does not occur exactly at VpVD = 1. Many investigators have simulated the cases presented in Figure 2.1 and a review is provided in the next section. However, here some general comments are needed on the experimental data themselves as well as any computational effort to simulate these shockinduced combustion cases. The shadowgraph from Lehr for M=6.46 is shown in Figure 2.2 with points used to denote the shock and reaction fronts. It has been verified (Lehr, personal communication, 1995) that hemispherically nosed projectiles with cylindrical bodies have been used in the experiments. However, in the image the nose of the projectile does not appear to be a true hemisphere. One possible source for this alteration is that the projectile is simply yawed. Such a yawing would also be evident in a corresponding change to the projectile's baseline seen in the images and this is not the case. There is a slight alteration to the baseline in the images but not nearly on the order of the alteration of the nose. After consultation with Lehr it has been verified that the alteration is due to the optical system (Lehr and Maurer, 1969) and he cautions against using the data for precise measurements. However, the data can be used qualitatively to benchmark the computational models to see if the flow field features are being captured. And given the fact the current study is focused on the interplay of the fluid dynamics and chemical processes, Lehr's experiments offer a solid set of data in which the scales of each of these aspects are comparable. As is evident by now, the successful simulation of a flow scenario such as Lehr's experiments (Lehr, 1972) or the reactive munition problem requires a multitude of models for the various phenomenological aspects. Each model not only embodies the particular Figure 2.2. Shadowgraph for M=6.46 with the shock and reaction front denoted. phenomenology framework used to represent either the physics or chemistry but also their numerical traits. For instance, the reaction process between hydrogen and oxygen in reality occurs through countless intermediate processes. However, to represent this phenomenon with a numerical model, a distinct number of reaction steps must be selected and reaction rates for these steps must be specified. The selection of a particular model, i.e. the number of steps and rates, inherently specifies the chemical scales produced by the model. These scales may or may not match those which exist in reality. If a discrepancy exists between the computation and experimental data, the model can be altered until a good match is made. This can increase the confidence in the ability of the model to reproduce key phenomenological characteristics such as the chemical time scales in the case of the reaction model. However, such alteration should be done with caution. For example, many computations have been carried out by previous investigators for the shockinduced combustion scenario captured by Lehr (1972). In these earlier studies, several reaction mechanisms have been used as well as a variety of grids. In most cases, an individual study uses a single reaction mechanism and the grid is adapted or other numerical aspects are altered until good agreement with experimental results is achieved. The accuracy of the computation is dependent not on a single model but the combination of all models. It is to become clear that for this class of problems, the phenomenological and numerical aspects must be considered in conjunction. 2.1 Reaction Mechanisms for Hydrogen and Oxygen One of the predominate differences among the earlier computational studies into the shockinduced combustion scenario is the representation of the chemical reactions. Therefore, here a description of various reaction mechanisms used for hydrogen and oxygen is provided. As detailed earlier, the chemical processes in the problems of interest have time scales such that a finite rate chemistry model is needed. The chemical reaction processes involved are incorporated into the governing equations by way of the source terms to be defined later. There have been many reaction mechanisms proposed for the reaction between hydrogen and air. Here, just those used in the previous computational studies to be reviewed and those used in the current study are presented. They range in complexity from a 2step global model proposed by Rogers and Chinitz (1983) to a 32step model proposed by Jachimowski which is based on his original 33step model (Jachimowski, 1988) but incorporates a slight modification he and others suggest (Oldenborg et al., 1988). These mechanisms are presented in Tables 2.1 through 2.7. The constants are used to calculate the rate of the individual reaction using the Arrhenius expression P Ean kn = AnTb exp[ E (2.1) which is discussed in more detail later. In the review, a 19step mechanism is mentioned. This mechanism is composed of the first 19 steps of the Jachimowski. Note, the two 8step mechanisms presented use the same steps but different rate constants. Table 2.1. Twostep hydrogenoxygen reaction mechanism from Rogers and Chintz (1983). n Reaction An bn Ean / Ru 1 H2 + 02 = 20H (8.917) + 31.433/(28.950) x 1047 10 2449.77 2 H2 + 20H = 2H20 (2.0 + 1.333/(P0.833() x 1064 13 21400.87 ( = equivalence ration, set to 1 in current computations Units are moles, seconds, centimeters, and Kelvin Table 2.2. Sevenstep hydrogenoxygen reaction mechanism used by Shang et al. (1995) based on eighteenstep mechanisms from Drummond et. al. (1986). n Reaction An bn Ean / Ru 1 H2 + 02 = 20H 1.7 x 1013 0 24233. 2 H + 02 = OH + O 1.42 x 1014 0 8254. 3 H2 + OH = H + H20 3.16 x 107 1.8 1525. 4 O + H2 = OH + H 2.07 x 1014 0 6920. 5 20H = H20 + O 5.5 x 1013 0 3523. 6 H + OH + M = H20 + M 2.21 x 1022 2 0. 7 2H + M = H2 + M 6.53 x 1017 1 0. Units are moles, seconds, centimeters, and Kelvin 2.2 Previous Computational Studies of Lehr's Experiments Both steady and unsteady combustion cases captured in shadowgraphs by Lehr (1972) have been studied by several investigators using numerical techniques. The actual Table 2.3. Eightstep hydrogenoxygen reaction mechanism from Moretti (1965). n Reaction An bn Ean / Ru 1 H + 02 = OH + O 3.0 x 1014 0 8810. 2 O + H2 = OH + H 3.0 x 1014 0 4030. 3 H2 + OH = H + H20 3.0 x 1014 0 3020. 4 20H = O + H20 3.0 x 1014 0 3020. 5 H2 + M = 2H + M 1.85 x 1017 1 54000. 6 H20 + M = OH + H + M 9.66 x 1018 1 62200. 7 OH + M = O + H + M 8.0 x 1016 1 52200. 8 02 + M = 20 5.8 x 1016 1 60600. Units are moles, seconds, centimeters, and Kelvin Table 2.4. Eightstep hydrogenoxygen reaction mechanism from Moretti (1965) with rate constants proposed by Evans and Schexnayder (1980). n Reaction An bn Ean / Ru 1 H + 02 = OH + O 1.2 x 1017 .91 8315. 2 O + H2 = OH + H 7.5 x 1013 0 5586. 3 H2 + OH = H + H20 2.0 x 1013 0 2600. 4 20H = O + H20 5.3 x 1012 0 503. 5 H2+M=2H+M 5.5 x 1018 1 51987. 6 H20 + M = OH + H + M 5.5 x 1018 1.5 59386. 7 OH + M =O + H +M 8.5 x 1018 1 50830. 8 02 + M = 20 7.2 x 1018 1 59340. Units are moles, seconds, centimeters, and Kelvin parameter range studied by Lehr exceeds that presented by Lehr (1972) and these additional data have been provided by Lehr (personal communication, 1995). Here a brief review of previous numerical studies is provided. In addition to the impact the reaction mechanism has on the computed results, the numerical aspects including solution method and grid resolution strongly affect the computed flow fields. Therefore, a description of these facets of each study is included in the review along with key results. Yungster et al. (1989) have computed the M=6.46, hydrogenair case captured by Lehr (1972) as well as the shots into hydrogen and oxygen presented by Lehr (1972). The 8step reaction mechanism for hydrogen and air presented in Table 2.4 by Moretti (1965) is used but the rate constants are replaced by values given by Evans and Schexnayder (1980). Table 2.5. Eighteenstep hydrogenoxygen reaction mechanism from Drummond et al. (1986) used by Ahuja and Tiwari (1993). n Reaction An bn Ean / Ru 1 02 + H2 = 20H 1.70 x 1013 0 24245.93 2 02 + H = OH + O 1.42 x 1014 0 8258.22 3 H2 + OH = H20 + H 3.16 x 107 1.8 1525.76 4 H2 + O = OH + H 2.07 x 1014 0 6923.81 5 20H = H20 + O 5.50 x 1013 0 3524.85 6 OH + H + M = H20 + M 2.21 x 1022 2 0. 7 2H + M = H2 + M 6.53 x 1017 1 0. 8 H + 02 + M = HO2 + M 3.20 x 1018 1 0. 9 HO2 + OH = H20 + 02 5.00 x 1013 0 503.55 10 HO2 + H = H2 + 02 2.53 x 1013 0 352.48 11 H02 + H = 2 OH 1.99 x 1014 0 906.39 12 H02 + O = OH + 02 5.00 x 1013 0 503.55 13 2HO2 = H202 + 02 1.99 x 1012 0 0. 14 HO2 + H2 = H202 +H 3.01 x 1011 0 9416.39 15 H202 + OH = H20 +HO2 1.02 x 1013 0 956.75 16 H202 + H = H20 + OH 5.00 x 1014 0 5035.50 17 H202 + O = H2 + OH 1.99 x 1013 0 2970.95 18 H202 +M = 20H 1.21 x 1017 0 22911.53 Units are moles, seconds, centimeters, and Kelvin The governing equations are solved using a fully coupled equation set and the total variation diminishing (TVD) algorithm by Yee and Shinn (1989), also known as the "point implicit TVD MacCormack." Grids of resolution 42x44 and 32x32 are used with the only difference reported between the solutions being the "thickness of the captured shocks." Good agreement between the computed shock and reaction front location and the experimental value is reported. However, only a limited region encompassing the nose is modeled. Lee and Deiwert (1989) have computed the M=6.46 case using the F3D code by Ying (1986) which is based on the implicit flux vector splitting method. The governing equations are solved in a looselycoupled approach which assumes that in the Jacobians of the fluxes, the effective specific heat ratio and pressure are not affected by the changes in species mass fractions. The reaction mechanism used is the same as that of Yungster et al. (1989). Some of the calculations are carried out with the rate coefficients reduced one and two orders of magnitude in an effort to quantify the impact of reaction rates on the flow field. They do find Table 2.6. Nineteenstep hydrogenoxygen reaction mechanism used by Oldenborg et al. (1990). n Reaction An bn Ean / Ru 1 02 + H2 = H + H20 2.16 x 108 1.51 1726. 2 H + 02 = O + OH 1.91 x 1014 0. 8273. 3 O + H2 = H +OH 5.06 x 104 2.67 3166. 4 H + HO2 = H2 + 02 2.50 x 1013 0. 349. 5 H + HO2 = 20H 1.50 x 1014 0. 505. 6 O + HO2= OH +02 2.00 x 1013 0. 0. 7 OH + HO2 = H20 + 02 2.00 x 1013 0. 0. 8 H + 02 + M = HO2 + M 8.00 x 1017 .8 0. 9 H + OH + M = H20 + M 8.62 x 1021 2. 0. 10 2H + M = H2 + M 7.30 x 1017 1. 0. 11 H+O+M=OH+M 2.60 x 1016 .6 0. 12 20 + M = 02 + M 1.14 x 1017 1. 0. 13 20H = O + H20 1.50 x 109 1.14 0. 14 20H + M = H202 +M 4.73 x 1011 1. 3206. 15 OH + H202 = H20 + H02 7.00 x 1012 0. 722. 16 O + H202 = OH + HO2 2.80 x 1013 0 3220. 17 H + H202= H2 + HO2 1.70 x 1012 0 1900. 18 H + H202= H20+ 02 1.00 x 1013 0 1800. 19 2HO2 = H202 + 02 2.00 x 1012 0 0. Units are moles, seconds, centimeters, and Kelvin that the slower rates do impact both the computed shock location and the temperature field behind the shock. Their results are discussed further in the context of the results from the current study. A grid of resolution 57 x 41 over the nose region is used and the results are reported to be very similar to those by Yungster et al. (1989). Wilson and MacCormack (1990) compute the steady combustion cases at M=5.11 and 6.46 using the 32step reaction mechanism given in Table 2.7. The mechanism, as presented by Wilson and MacCormack (1990), is in a slightly different order and there is a typographical error in the preexponential constant for reaction 21. The constant An for reaction 10 is different than that used by Jachimowski (1988). Wilson (personal communication, 1997) reports this alteration has been made based on the rate information available in other sources such as the work of Warnatz (1984). In the current study, the rate for step 10 as defined by Wilson and MacCormack (1990) is used to facilitate comparisons. A fullyimplicit fluxsplitting technique based on the work by MacCormack (1985) and Table 2.7. Hydrogenoxygen reaction mechanism from Jachimowski (1988). n Reaction =An_ bn Ean / Ru Units are moles, seconds, centimeters, and Kelvin Thirdbody efficiencies: (6) H20 =6.0, (7) H20= 6.0, H2 : 2.0, (8) H20= 5.0, (9) H20 = 16.0, H2 = 2.0, (19) H20 = 15.0 Candler (1988) is used along with StegerWarming flux splitting (Steger and Warming, 1981). The governing equations also include vibrational nonequilibrium effects; however, Wilson (personal communication, 1996) reports these effects to be negligible. The study evaluates increasing the induction delay by changing the reaction rate for H + 02 OH + O, step 2 in Table 2.7. The rates (a) k2 = 2.60 x 1014 exp( 8459.64/T) (b) k2 = 1.91 x 1014 exp( 8277.35/T) (2.2) (c) kf = 1.2 x 107T91 exp( 8315.12/T) H2 + 02 = HO2 + H H + 02= OH + O O + H2 = OH + H OH+ H2 = H + H20 20H = O + H20 H + OH + M = H20 + M 2H + M = H2 + M H+O+M=OH+M H + 02+ M = HO2 + M 20+ M = 02+ M HO2 + H = 2 OH HO2 + H = H20 + O HO2 + 0 = 02 + OH HO2 + OH = H20 +02 2HO2 = H202 +02 H + H202 = H2 + HO2 O + H202 = OH + HO2 OH + H202 = H20 + HO2 H202 + M =20H + M 2N + M = N2 + M N + 02 + NO + O N + NO = N2 + O N + OH = NO + H H + NO + M = HNO + M H + HNO = NO + H2 O + HNO = NO + OH OH + HNO = NO + H20 HO2 + HNO = NO + H202 HO2 + NO = NO2 + OH H + NO2 = NO + OH O + NO2 = NO + 02 M + NO2 = NO +O + M L ____________________ I 1.0 x 1014 2.6 x 1014 1.8 x 1010 2.2 x 1013 6.3 x 1012 2.2 x 1022 6.4 x 1017 6.0 x 1016 2.1 x 1015 6.0 x 1013 1.4 x 1014 1.0 x 1013 1.5 x 1013 8.0 x 1012 2.0 x 1012 1.4 x 1012 1.4 x 1013 6.1 x 1012 1.2 x 1017 2.8 x 1017 6.4 x 109 1.6 x 1013 6.3 x 1011 5.4 x 1015 4.8 x 1012 5.0 x 10"11 3.6 x 1013 2.0 x 1012 3.4 x 1012 3.5 x 1014 1.0 x 1013 1.16 x 1016 0 0 1 0 0 2 1 .6 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. .75 1 0 .5 0 0 .5 0 0 0 0 0 0 28198.8 8459.64 4481.6 2593.28 548.87 0. 0. 0. 503.55 906.39 543.83 543.83 478.37 0. 0. 1812.78 3222.72 720.08 22911.53 0 3172.36 0 0 302.13 0 0 0 0 130.92 755.33 302.13 33234.3 are evaluated and it is determined choice (c), the rate suggested by Warnatz (1984) gives the best match for Lehr's experiments at M=5.11. This same rate is used to compute the M=6.46 case resulting in a good match of both computed shock and reaction front to the experiments. The M=6.46 case is computed on a 321 x 65 adapted mesh. In the current study, the 32step mechanism with the rates per Warnatz is labeled the modified 32step model. Wilson and MacCormack also evaluate altering the time of reaction by changing the rate constants for H + OH + M H20 + M, step 6 in Table 2.7. This alteration affects the heat release rate and the time required to reach the peak temperature. Reducing the original rate per Jachimowski (1988) by a factor of 4 and using choice (a) above for H + 02 ; OH + O produces no significant change in shock and reaction front location. Ahuja and Tiwari (1993) use an 18step mechanism, given in Table 2.5, to compute the M=6.46 case. The mechanism is reported to be based on the work of Jachimowski (1988) and matches the mechanism used in Drummond's work (1988) except for inconsistencies in steps 9 and 10 given by Ahuja and Tiwari (1993) which are believed to be simply typographical errors. The rate constants are not provided by Ahuja and Tiwari (1993) but are presented here based on the data available from Drummond (1988). The equations are solved using MacCormack's method and a LaxWendroff type scheme is used for the fluxes. The grid resolution used is 197x152 over a domain comprising a quarter sphere. The computed flowfield is reported to show a shock and reaction front but it is difficult to make a detailed comparison to experiments given the limited computational domain. Matsuo et al. (1993) use a 19step reaction mechanism based on 32step used by Wilson and MacCormack (1990). In the 19step mechanism the nitrogen reactions are neglected and as mentioned, the specific reactions used are steps 1 through 19 of Table 2.7. The governing equations are integrated explicitly, and the chemical reaction source term is treated in a linearly pointimplicit manner. Yee's nonMUSCL TVD upwind explicit scheme is used and computations have been carried out for the M=4.18 and 4.79 cases captured by Lehr (1972; personal communication, 1995). A grid of resolution 161x321 over a quarter sphere domain is used. For the M=4.79 case, computations show the frequency of the unsteadiness to be 725 kHZ where the reported experimental frequency is 712 kHZ. For the M=4.18 case, a frequency of 159.8 kHz is computed where the experiments of Lehr (1972) report a frequency of 150kHz. Also important in this study is their characterization of the chemical time scales which is discussed later. Sussman (1993) develops a source term modification scheme for use in combustion computations. Using a model problem, the error associated with the representation of the reaction source term is estimated. A scaling factor based on this error estimation is used in computing both a quasione dimensional shockinduced combustion scenario as well as the M=6.46 case captured by Lehr. This source term modification is discussed in detail in a later chapter and is evaluated in the current study. The 19step reaction mechanism presented in Table 2.6 (Oldenborg et al., 1990) is used. The fluid dynamic terms are solved explicitly using the HartenYee (Yee, 1989) TVD algorithm and the source terms are integrated using a point implicit method. A grid resolution of 44x49 is used. In the quasione dimensional case, the source term modification results in a reduction in the number of grid cells needed to resolve the induction zone. This same scheme when used in Lehr's problem does not give improved results without the addition of a grid adaptation technique. Wilson and Sussman (1993) use the 32step mechanism as in Wilson and MacCormack (1990) with the rate suggested by Warnatz (1984). A logarithmic form of the species conservation equations is used for which a description is provided later. A comparison between StegerWarming fluxvector splitting (Steger and Warming, 1981) based on work by MacCormack (1985) and Candler (1988) and HartenYee (Yee, 1989) upwind TVD scheme is made. The convective terms have been treated explicitly and source terms implicitly. Both steady and unsteady cases have been computed. For the steady case of M=6.46, a 52x52 grid is used with the logarithmic form of the equations. The computed shock and reaction fronts from both flux schemes compare well with experiments. For the unsteady case of M=4.79, a 375x161 grid is used but the standard form of the species conservation equations is solved. The frequency of the unsteadiness is calculated to be 530 kHz where experiments are reported to be 720 kHz. The reported experimental value is slightly different from the value of 712kHz reported by Matsuo et al. (1993). The computations are repeated using the original rate for H + 02 OH + O per Jachimowski and the frequency is computed to be 820 kHz. Matsuo et al. (1993) cite a private communication with Wilson and Sussman that after running the computation for a longer time, the results match the experimental frequency better than the 820 kHz cited. The findings that a different rate for the H + 02 OH + O reaction are needed for the unsteady and earlier reported steady cases is one area of concern regarding the available reaction mechanisms. Matsuo and Fujiwara (1993) compute the unsteady shockinduced combustion flow using the conditions from Lehr's experiments for M=4.18. The diameter of the nose is changed to generate different cases. A twostep mechanism using progression variables is used and the parameters are based on shocktube data for 2H2 + 02 + 7Ar. No direct comparison to Lehr's data is made but the computations are used to help explain the phenomena causing the instabilities. The solution scheme used is the same as in Matsuo et al. (1993). No details on the grid resolution are given. Sussman (1994) presents results of computations for both the low and high frequency cases studied by Lehr. The equations are solved using the technique of Wilson and Sussman (1993) including the logarithmic form of the species governing equations. All computations are performed on a 192x252 grid over a quarter sphere domain and the 19step mechanism of Table 2.6 (Oldenborg, et al., 1990) is used. Excellent agreement between the computed and experimentally measured frequencies is reported for all cases except at the lowest and highest projectile velocities. The discrepancy at the lowest frequency case is attributed to the fact that the energy release front is closest to the body at this velocity and either a finer grid in that region is needed or the neglect of the viscous effects is most prevalent here. The discrepancies at the high frequency end of the spectrum are suspected to be due to the difficulty in measuring the scales from the experimental results given their small amplitude. Also of interest in Sussman's study is the discussion of the mechanism determining the instabilities. Sussman reports finding an instability for projectile velocities of 2119 m/s where Lehr (1972) found no evidence of unsteadiness. Sussman finds that this transition point from an unsteady to a steady combustion scenario it directly related to the point at which the induction time of the reaction becomes smaller than the energy release time. This criterion is different from the projectiletodetonation velocity criterion mentioned earlier. The frequency of the instabilities is directly related to the induction time and given the fact that the computed values compare well with the experiments over the majority of the velocity regime, Sussman deduces that this aspect of the chemical process is being adequately modeled by the reaction mechanism being used. However, since the transition point from unsteady to steady combustion is over predicted, the energy release rate produced by the model is greater than that in the actual experiments. This conclusion that the ratio of energy release time, or reaction time, to the induction time is a key parameter is in agreement to the findings of Matsuo et al. (1993) and those of McVey and Toong (1971). Sussman also presents detailed computations of a onedimensional shockinduced combustion scenario and demonstrates the advantages of solving the logarithmic form of the species equations. Yungster and Radhakrishanan (1994) computed the unsteady cases captured by Lehr at M=4.18, 4.48, and 4.79. The 19step mechanism of Matsuo et al. (1993) is used with the original rates per Jachimowski. The governing equations are solved fully coupled using Yee's 2nd order TVD scheme (Yee, 1989) and a grid resolution of 220x220 over half the hemispherical nose region. Good agreement between computed frequencies of instabilities and those measured by Lehr are achieved. No computation of the steady case is performed. Shang et al. (1995) compute Lehr's M=6.46 case using the FDNS code (Chen et al., 1993), a pressure based finite difference code using the TVD scheme of Chakravarthy and Osher (1985). An operatorsplitting scheme is used where the chemical kinetics source terms are integrated implicitly with an ODE solver in a predictor step and then coupled with the fluid dynamics in a corrector step (Rhie et al., 1993). The M=6.46 case is computed on a uniform (87x97) and adapted grid (44x49) using a 7step mechanism, presented in Table 2.2, based on the 18step mechanism by Drummond et al. (1986) presented in Table 2.5. This mechanism is also documented in other work by Drummond (1988, 1990). Various interpolation techniques for defining the species concentrations in the predictor step are evaluated. Good comparisons with experimental results are achieved when either an adapted grid or a high order interpolation is used. These interpolation techniques are discussed in detail later and are evaluated in the current study. Matsuo and Fujii (1995) compute unsteady cases using the conditions from experiments of Lehr for M=4.79. The same solution scheme of Matsuo et al. (1993) and a 2step reaction model similar to that of Matsuo and Fujiwara (1993) are used. The 19step mechanism (Matsuo et al., 1993) is used to perform a time integration of a reduced equation for the species conservation with no convection or diffusion. Different activation energies for the exothermic reaction of the 2step mechanism are evaluated to generate different heat release rates or time of reaction. Results show the reaction time directly impacts the flowfield. Of the three activation energies evaluated, the one with a slow and the one with a medium heat release rate generate periodic unsteady results. A rapid release (Ea=0) produces a larger, more random disturbance. Matsuo and Fujii (1996a) perform computations using the experimental conditions of Lehr (1972) and Ruegg and Dorsey (1962). A grid resolution of 401x401 is used to represent a quarter sphere domain and the 19step mechanism used by Matsuo et al. (1993) is used. The pressure, temperature and velocities of the experiments are used but the diameter of the projectile is altered to generate steady, low frequency unsteady, and high frequency unsteady cases. A zerodimension time integration of the governing equations is used as in Matsuo and Fujii (1995) to determine correlation between a defined Damkohler parameter and the resulting combustion situation. The integration uses the conditions along the stagnation line and assumes internal energy and total density are constant. This results in a temporal profile depicted in Figure 2.3. The definition used for the Damk6hler parameter is S= tf/t = [d/a2]/[T2/(dT/dt)max (2.3) where d is the projectile diameter and the quantities with subscript 2 denote values behind the shock. Based on their results, a critical value of T = 80 is defined below which the unsteadiness is of a high frequency and above which it is a low frequency instability. No universal parameter to classify steady as well as unsteady cases is found. Such classification may require using the parameter defined in Matsuo and Fujii (1996a) as well as the critical parameter of Vp/VD as determined by Lehr (1972) where Vp is the projectile speed and VD is the detonation velocity of the mixture. Matsuo and Fujii (1996b) compute unsteady cases from experiments of Ruegg and Dorsey (1962). The 19step mechanism used by Matsuo et al. (1993) is employed and the governing equations are solved explicitly with the chemical reaction source term treated in initiation by shock dT/dt T2 Tind Tl1 t Figure 2.3. Depiction of the temporal profile from a time integration of the governing equations using the conditions behind the shock on the stagnation line. a linearly pointimplicit manner. Yee's nonMUSCL TVD upwind explicit scheme is used on a 601x601 grid over a quartersphere domain. Results from grids of 201x301, 401x401, and 801x801 showed similar flow features. No direct comparison to experimental results is presented but conclusions as to the mechanisms causing the instabilities are made. As is clear by the review here, the steady and unsteady shockinduced combustion scenarios captured by Lehr have been used in many computational studies. The current study computes the case at M=6.46 in a stoichiometric mixture of hydrogen and air. Therefore, those previous studies which have computed this same case are summarized in Table 2.8. Included in this summary are the particular fluid dynamics model and reaction model used as well as the grid resolution employed. Also, a depiction of the flow domain computed and resulting fronts are given. These images have been generated using the computed density field presented in each reference and orientating a line depicting the shock and energy release fronts based on this data. It is clear that those computations which do not cover a substantial portion of the projectile flow domain do not provide enough information to judge the accuracy of the computation. Therefore, the current study computes the entire length of the projectile as depicted in Figure 2.4 which is a larger domain than any of the previous computations. Also denoted in Figure 2.4 are the domains used in earlier computational studies in which the M=6.46 steady combustion case has been simulated. The utility of the larger domain is to be evident later. 2.3 Summary As detailed here, the shockinduced combustion problem experimentally studied by Lehr (1972) and computationally studied by several investigators is to be used as an initial problem in developing the needed models for simulating the munition concepts discussed earlier. As is evident in the review, both the phenomenology of the various models and the numerical aspects are key in determining the accuracy of a simulation for this class of problems. For example, the chemical time scales of induction and reaction time are found to be key. However, this trend in the physical and chemical processes is only detectable if Table 2.8. Summary of previous computational studies into steady shockinduced combustion. Lehr (1972) Experiments Yungster et al. (1989) Mechanism: 8step, Table 2.4 Grid : 42x44 Scheme : Fully implicit Lee and Deiwert (1989) Mechanism: 8step, Table 2.4 Grid: 57x44 Scheme : Fully implicit Wilson and MacCormack (1990) Mechanism : 32step, Table 2.7 Grid : 321x65 adapted Scheme : Fully implicit Ahuja and Tiwari (1993) Mechanism: 18step, Table 2.5 Grid: 197x152 Scheme: MacCormack Sussman (1993) Mechanism: 19step,Table 2.6 Grid : 44x49 Scheme : ExplicitImplicit Notes : Comparison improves with adapted grid. Wilson and Sussman (1993) Mechanism: 32step, Table 2.7 Grid : 52x52 Scheme : ExplicitImplicit Notes : Solve logarithmic form of equations. Shang et al. (1995) Mechanism: 7step, Table 2.2 Grid: 87x97 Scheme : ExplicitImplicit Notes : Comparison improves with extapolation scheme. (1) Extent of domain used by Yunster et al. (1989) Lee and Deiwert (1989) Ahuja and Tiwari (1993) (2) Extent of domain used by Wilson and MacCormack (1990) Sussman (1993) Wilson and Sussman (1993) Shang et al. (1995) (3) Extent of domain used in current study Figure 2.4. Computational domain used in previous investigations and throughout current study. an adequate spatial discretization and reaction model is used. In the current study, not only are the computations compared to experiments, but the change in the dependent variables as a result of alternative models is used to help quantify the roles of the individual models on the accuracy of the complete computational tool. This not only aids in determining the impact these aspects have on the computation but also helps quantify the phenomenological and numerical aspects involved in solving a class of problems such as the fuel cell scenario. This study uses the steady combustion case shown in Figure 2.2 to assess the aspects of modeling highspeed reacting flows. If the conditions for M=6.46 are used in the Damkihler parameter definition of Matsuo and Fujii (1996a), Eqn. (2.3), a value of 18 is calculated. This supports the use of this flow scenario to study the computation of reacting 36 flows where 2) = 0(1). As mentioned, the optical distortion in the image raises some concerns in a onetoone comparison between the computations and image. However, instead of reinterpreting the image, the locations denoted in Figure 2.2 are used as in the case of previous investigations. It should also be noted that the experimental data is not symmetric about the stagnation line which suggests the projectile had some slight angle of attack during its flight. Results for the steady shockinduced combustion case are discussed following the description of the computational models developed here. CHAPTER 3 GOVERNING EQUATIONS AND NUMERICAL TECHNIQUES In this chapter the governing equations and various numerical techniques needed to solve the current class of problems are detailed. After the presentation of the governing equations, the individual approaches are discussed. For each model there are various options as to the particular phenomenological and numerical form. Where appropriate these options are discussed. 3.1 Governing Equations The governing equations for the current problem are the axisymmetric form of the NavierStokes equations with appropriate source terms. The integral form of the governing equations is f Q dT + F.dS = dT (3.1) where Q is the dependent variables in volume Y, F is the flux of those variables across the surface (S) which defines the volume T and Q is any source term per unit volume for the dependent variables. The flux, F, is a nonlinear function of the dependent variables Q. The governing equations once discretized for a general curvilinear coordinate system are a J Q a J (F, Fy) 9 J (GI Gy) at + + J (H Hy) = J 2Q (3.2) where the dependent variables are Q = [Q, Qu, v,E,QXi ....QNs1] (3.3) The variables Q, u, v, and E are respectively the density, x component of velocity, y component of velocity, and energy. Also, ai is the mass fraction of ith species with the fluid being defined by NS total species. Note that the mass fraction of the NSth species is not NS explicitly modeled since the total density is included and the relationship p = Qpti holds. Another modeling option is to not solve an equation for total mass but to solve an additional species conservation equation for the NSth species and calculate the total mass using the summation given above. With the addition of the relationship between the individual species densities and the total mass, both implementations require attention to numerical roundoff issues. This is evident in previous investigations such as the work of Sussman (1993) and Wilson and Sussman (1993) in which it is noted that species densities can sometimes become negative due to numerical errors and that limiters must be used. In the curvilinear directions and r, the inviscid flux vectors are QU OV puU + SxP puV + YTxP QvU + yP QvV + ryP FI = U(QE + P) GI = V(QE + P) (3.4) 0aU galV QaNs_U PQNSIV where P is pressure, and the viscous flux vectors are 0 0 Jxtxx + lytxy Txtxx + TlyXxy Jxtxy + yTyy x'Txy + 1y'tyy FV = xx + MyIy G = lxnx + ryy (3.5) QalUl QaiVi eNS i NSl QaNsIVNs The third dimension is incorporated using the axisymmetric terms Ogv 0 Quv txy pv2 Tyy Too H = v(E + P) H (3.6) H y Qalv eQiv1 Q(NS I QONSlVNS where y is the distance from the line of symmetry and the vector of chemical reaction source terms is Q = [0,0,0,0,60,,...,oNS_ (3.7) The grid Jacobian J, an indication of the cell volume, and the contravariant velocities are defined as J = xyr) xy U = xu + yV (3.8) V = lTxu + T]yV where x y x X J y T) y, x (3.9) fix j 11y = J The components of the viscous flux vectors are defined as Iix = uTxx + vTxy qx Iy = utxy + VTyy qy Ty = 2 + + + (3.10) aay (ay x Y ee = 2 + (ay ax Y) where the bulk viscosity is represented using Stokes' hypothesis, = The heat flux components are NS NS q= k + h(aG^), qy= k T + _hi(aivi) (3.11) qx ax + ay i=l i=l where k is the thermal conductivity and ai i and aiVij are the diffusion velocities defined in the next section. 3.1.1 Molecular Diffusion A A The vector components Ui iV are the contravariant diffusion velocities for species i and are represented using Fick's law as described by Anderson (1989), and Libby and Williams (1994). The exact form of the modeled terms are A A A A Ui = xUi + yVi Vi = TIxUi + TlyVi A Oa aO (3.12) aui = D aiiv = Di ay The use of Fick's law is made possible by some simplifying assumptions as to the molecular diffusivity of individual species. A more exact representation of the diffusion of a single species i in a mixture composed of NS species is to calculate the diffusivity of the species i as a function of its diffusivity in each of the other species. One level from this in simplification is to define a background gas or the gas which comprises the majority of the mixture. This allows the diffusivity of each species in the predominant gas to be calculated and used for Di. For example, combustion problems involving hydrogen and air typically contain a large amount of nitrogen which can be defined as the background gas. The next level in simplification is to assume a single diffusion coefficient such that Di=D for all species, and that the binary diffusion can be related to the viscosity using the Schmidt number, Sc (3.13) QD This simplification, used in the current study, has proven reliable in many computational studies and is considered acceptable when all the species involved have approximately the same molecular weight. The only concern in the context of the present study is the hydrogen involved given its molecular weight is one order of magnitude less than the other species. However, this simplification is employed given the predominance of the other transport mechanisms for the species in the current flow scenario and considering the reduction in the computational requirement. Also, in the context of diffusion itself, the molecular diffusivity is secondary to the turbulent enhancement which is discussed later. 3.2 Equation of State and Temperature The equation of state used in the current computations is derived by assuming the ideal gas equation is valid for each individual species and has the form (Anderson, 1989) NS P = RuT ZMi (3.14) i=l where Ru is the Universal Gas Constant and Mi is the molecular weight of species i. There are flows for which an alternate equation of state is more applicable and from a computational model standpoint, such alternate equations simply require an alternate definition of terms for which the derivative of pressure with respect to the conserved variables are needed. These modeling aspects are presented in a later section. The temperature during the calculations must be extracted from the internal energy using the relationship NS e=E u = ah i=1 T (3.15) hi = hf0 + CpidT TR where TR is the reference temperature and hf the heat of formation of species i at TR. The specific heat, Cpi, of each species is a known function of temperature and is available from various references such as the JANNAF data (Stull and Prophet, 1971). The dependence of Cpi on temperature can be represented in various forms ranging from a constant up to a high order polynomial. In the shockinduced combustion computations presented here, a third order polynomial of the form Cpi = c1 + c2T + c3T2 + C4T3 (3.16) is used. The coefficients have been determined by fitting curves to the JANNAF data up to 6000K and are given in Table 3.1. Higher order polynomials have been used by previous Table 3.1. Coefficients used to calculate specific heats for the various species. n c1 C2 c3 C4 hf 02 6.72433 1.75132 x 103 3.35085 x 107 2.40485 x 1011 0. H2 6.18833 1.28193 x 10 1.59252 x 107 9.05493 x 1012 0. OH 6.96378 5.90649 x 104 2.75234 x 108 1.03921 x 1011 9.432 H20 6.9803 3.6993 x 103 6.55318 x 107 4.07621 x 10" 57.798 O 5.25211 2.0219 x 104 3.79786 x 108 0 59.559 H 4.968 0. 0. 0. 52.1 HO2 7.44933 7.3918 x 103 2.8673 x 106 3.955 x 10'1 5. H202 6.6606 1.49214 x 102 7.7563 x 10 1.3478 x 109 32.53 N2 6.51307 1.65586 x 10 3.74607 x 107 2.8942 x 10" 0. Cpn in cal / (mole K) hf in kcal / mole TR = 298 K investigators but here, the thirdorder has been found to provide a good match using a single set of constants across the temperature range. Given this representation of specific heats, an iterative NewtonRaphson procedure must be used to extract the temperature in each cell at each time step as detailed by Molvik and Merkle (1989). The procedure solves the equation em+ = em + e mAT (3.17) aT) where m denotes the iteration step. The initial guess for the temperature in each cell at a given time step is simply the temperature at the previous time step and convergence is typically achieved in only a few iterations if an overrelaxation scheme is used. Even though convergence is typically rather fast, it is clear to see that for a computation with a larger number of computational cells, this added aspect for reacting flows and for that matter high temperature flows can become costly. Therefore, for some applications a lower form for Cpi is worth using. In fact, if a linear representation for dependence of specific heat on temperature is used, the temperature can be calculated by solving a simple algebraic equation with no iteration. This representation enables the aspect of temperature dependence to be included with no added computational time. However, before using such a scheme, attention should be paid to the temperature range over which the computation is expected to vary to ensure the constants used best match this regime. 3.3 Thermal Conductivity and Viscosity For the current study, it is assumed that the thermal conductivity can be determined by specifying a Prandtl number where CtCp Pr = p (3.18) k This as the case with the diffusivity is a simplification which has proven applicable to the class of flows modeled here. Note, the Schmidt and Prandtl number are related through the Lewis number (Le = Sc/Pr) which is often near unity. In the current study Le= 1 is assumed and the Prandtl number is set to .72. We now have the binary diffusion coefficient and the thermal conductivity related to the viscosity of the mixture by a constant of order one. The viscosity of the mixture is modeled as a function of temperature and the species composition using Sutherland's law as described by White (1974), Vi T 5 Ti + Si oi T + Si (3.19) and Wilke's formula X j(ij j] = Nsij = + 8( 1)2()][( + 2 (3.20) j=1 where Xi is the molecular concentration of species i and Toi and Si are constants. The constants for use in Eqn. (3.19) are given in Table 3.2. One of the requirements for the fluid mechanics model is to limit the numerical errors in the models of quantities such as the Table 3.2. Coefficients used to calculate viscosity for the various species. diffusion of species and heat. Given the fact that both these quantities are related to the molecular viscosity with an order one constant, it is reasonable to conclude that the numerical errors in these two quantities are of the same order as numerical errors in molecular viscosity. The apparent errors in viscosity for typical computational models are not a function of the phenomenological model such as Eqn. (3.20) but are a result of numerical viscosity. Numerical viscosity is a direct function of the numerical scheme used in representing the governing equations. The particular scheme used in the current study is described next and the issue of numerical viscosity is addressed in that context. However, it is safe to assume that if the numerical viscosity is minimized, so will be the numerical errors in diffusion and conductivity. 3.4 Flux Vector Splitting Scheme For the solution of compressible flows, the most popular schemes for representing the inviscid fluxes are flux difference and flux vector splitting. These upwind schemes are attractive due to their ability to minimize numerical viscosity. A review of these various schemes for multispecies flows is provided in the work of Montagne et al. (1987), Glaister (1988), Cinnella (1989), Grossman and Cinnella (1990), Shuen et al. (1990) and Liou et al. (1990). Here the scheme used to define the inviscid numerical fluxes is the StegerWarming flux vector splitting (Steger and Warming, 1981) which has been extended to model multispecies flows by Liou et al. (1990). The flux vector splitting algorithm decomposes n o (Ns/m2) To (K) S (K) 02 1.919x 105 273 139 H2 8.411 x 106 273 97 HzO 1.120 x 10 350 1064 N2 1.663 x 105 273 107 F the inviscid fluxes into components based on the eigenvalues of the Jacobian AI = and aQ likewise for GI. The split fluxes take the form K* = kXK, + X:K2 + O:K3 (3.21) where Qu (u kxc) QV K _ 1 ( c2 1 (v kyc) yK h y 1 K2,3 y Q(ht ) (k22) Qal QaI L QaNS1 QaNS1 with the eigenvalues being Xk1 = (k Ikl) S= 13k c(3.23) 2 = Pk + clVkI k3 = Pk clVkl where Ok = kxu + kyV Ok = kxu + kyv k kY 1  (3.24) kx =kx ky k IVI =/k2x + k2 IVkl IVkl IVkl k +k( The above formulation gives K=FI when k=E and K=GI when k=ri. For the multispecies chemically reacting flow, c is the frozen speed of sound and y is the effective specific heat ratio. The calculation of c is detailed in the Appendix. In the finite volume formulation, the fluxes are evaluated at the cell faces and are either first or secondorder spatial representations. For the class of flows to be solved here, a secondorder scheme is desirable to limit the numerical viscosity. The flux at the face is a function of the states in the neighboring cells, (F),i+,j = F+(QL+Q + F(QR+ (3.25) If a firstorder spatial representation is used, then Q = Qi, Q ,i = Qi+,j .To achieve secondorder accuracy, a MUSCL approach described by Shuen (1992) is used in which cellcenter values are extrapolated to the interfaces. Also, to guard against the interpolation introducing any nonphysical extrema into the field in the region of large gradients, a limiter must be used. The formula for the neighboring states takes the form QL Qij + Q i+,j iQR, Q J (3.26) S',i+j = ij i+,J where the limiting function is S j = ,j minmod[AQ+i+I j,AQi+1,j 2 .( 3 .2 7 ) +j =i 2 minmod[AQ AQ ,j] with 2(Q i+ ,j i Q,j AQ+ij i+ j + ij( 2(Qi, Q(3.28) AQ'J fi,j + ilj Various limiters can be used with this implementation but here the minmod limiter is used where minmod [X, Y] = sign(X) max[0., min(IXI, Ysign(X))] (3.29) Note e ij is the celllength to provide weighting for nonuniform grid spacing. The same extrapolation procedure is carried out for the fluxes in r~ and can be performed on either the dependent or primitive variables. Previous investigations by Shuen (1992) have shown that using primitive variables gives better performance for flows with strong shocks and this is the method used here. Note, the secondorder accuracy is achieved by the assumption that the field quantities vary linearly across a cell. Higher order schemes are achievable if a higher order variation is assumed. However, in all cases, limiters must be used in conjunction with the higher order schemes to ensure against the introduction of any nonphysical extrema in regions of large gradients. The viscous fluxes are evaluated at the cell faces using central differencing and second order representations of the derivatives. 3.5 Chemical Reaction Source Term The particular reaction mechanisms used in the current study have been discussed in the earlier literature review of shockinduced combustion simulations. The production of each species is incorporated into the governing equations by way of the source terms in Eqn. (3.7). This source term is calculated using the finite rate chemistry model described here. The reaction between hydrogen and oxygen which drives the shockinduced combustion must be modeled by a specific number of reactions (NR) using the form NS kfi NS vi Xj <, j1j I i = NR (3.30) j=l kbi j=l where vi,' and vii" are the stoichiometric coefficients for species j in reaction i, and Xj is the molecular concentration of species j. The parameters kf, and kb, are, respectively, the forward and backward reaction rate constants for the ith reaction and can be expressed as a function of temperature with the Arrhenius expression ki = ATbexp Ea (3.31) where Ru is the universal gas constant, Ea, is the activation energy and Ai and bi constants for the ith reaction. The production rate of an individual species is calculated by summing the contribution due to each reaction NR NS NS j = M i i Vij ')x k XVn' kb, XV i (3.32) i=l n=l n=l The forward rate constants are available in various sources for most reactions of interest. The backward rate for a reaction is sometimes given explicitly using Arrhenius rates but is usually obtained using the equilibrium constant and the relationship kfl kei = . (3.33) The equilibrium constant can be expressed as a function of the stoichiometric coefficients and the Gibbs free energies (gn) of the species in the reaction through the formula NS ke,= [RuT] ,= exp RT (3.34) where NS Ag = (ni" Vni') gn (3.35) n=l The Gibbs free energy (gn) associated with each species is given as a function of temperature, gn = gn(T). Using the JANNAF data (Stull and Prophet, 1971), it is clear that a linear function of the form gn = gn,, + gn,T (3.36) is adequate. The values of the constants for the species to be used throughout the current work are provided in Table 3.3 as well as their molecular weights. From a solution perspective, the chemical reaction source terms are numerically stiff and require special attention which is discussed later. One note from the perspective of implementing a finiterate reaction model as described here is needed. The necessary routines can be coded up to handle a generic reaction mechanism. This requires multiple do loops to calculate the source term as well as the Jacobian of the source term (aQ2/Q) which is discussed later. This presents a problem when running such a code on a vector machine, or any machine for that matter, and literally orders of magnitude in computational time can be saved if the routine is coded up for a specific reaction mechanism. This is easily done if an analysis tool such as Mathmatica (Wolfram, 1988) is used to preprocess the needed equations and is well worth the time when large mechanisms or large grids are used. Table 3.3. Molecular weights and coefficients used to calculate the Gibbs free energy for various species. n Mn go gl 02 32. 0. 0. H2 2.016 0. 0 OH 17.08 9.1646 3.4 x 103 H20 18.016 59.29 1.35 x 103 O 16. 60.464 1.574 x 102 H 1.008 53.3078 1.396 x 102 HO2 33.008 4.195 1.194 x 102 H202 34.016 33.251 2.646 x 102 N2 28.016 0. 0. N 14.008 114.407 1.608 x 102 NO 30.008 21.468 2.94 x 103 HNO 31.016 23.2391 1.173 x 102 NO2 46.008 7.738 1.514 x 102 gn in kcal/mole 3.6 Special Treatment For Chemical Length Scales In the discussion of the fluxvector splitting scheme used in the current code, a detailed description is provided as to how to represent the numerical flux with a secondorder spatial model. Such a model reduces the numerical dissipation of field quantities and improves the accuracy of the scheme. The secondorder model is approximating the variation in variables across a cell as a linear function. As noted earlier, higher order schemes are achievable if a higher order variation is assumed. Second order schemes have proven to provide acceptable accuracy for simulations of compressible turbulent flows with little increase in computational requirements. However, this is due to the fact that the computational grid used in the simulations employ discretization of the appropriate size to capture the scales of the mean flow. With the addition of chemical reactions, the length scales which must be resolved can be orders of magnitude smaller. Of course one option is to construct the computational grid using a finer discretization. However, unless the location of key events such as the reaction front are known a priori this finer discretization must be used throughout a large portion of the flow field. An alternative is to employ a grid adaptation technique which is only advantageous if it is closely coupled with the fluid dynamics model. Both these options present an increase in the computational requirements. Other options of addressing the error associated with the small chemical length scales have been evaluated. Wilson and Sussman (1993) have used a transformed governing equation while Sussman (1993) and Shang et al. (1995) have tested source term modification techniques. These techniques take advantage of the known profiles of those species created in the chainbranching reactions. The variation of these species, such as OH, H, and 0, is more logarithmic in nature. The transformed governing equation used by Wilson and Sussman (1993) is discussed briefly but attention is paid to the schemes which are an alteration only to the source term. Typically the source term is calculated using the cell averaged value for temperature and mass fractions, a piecewise representation. The modified treatment calculates the source term using an altered value for the mass fractions, mij = 0(as)ij. The scheme used by Shang et al. (1995) differs from that of Sussman (1993) so both are described here. All the source term treatments are an effort to address the fact that the cell averaged source term is not necessarily equal to the source term calculated based on the cell averaged dependent variables. This is due to the highly nonlinear nature of the Arrhenius source term. 3.6.1 Transformed Governing Equations The modification option evaluated by Wilson and Sussman (1993) is a logarithmic transformation of the entire species conservation equation minus the diffusion term. The definition nn = ln() Qln(an) (3.37) is used with the total mass conservation equation to specify the transformed equation for species n, An+ a anuj= n (3.38) In practice these transformed equations are solved in conjunction with conservation equations for the individual elements to ensure the conservation of mass. Solutions by Wilson and Sussman (1993) for the M=6.46 case of Lehr using this set of transformed equations produce good results on somewhat coarse grids. Wilson and Sussman state this formulation is no less robust than the conventional form of the equations but do warn against the added distortion at contact discontinuities between different gases. This could be key in nonpremixed combustion scenarios. The computation time is cited to be 1020 % more per time step above the standard scheme. Also, of interest is the error trend when this scheme is used on coarse grids. If the shockinduced combustion flow scenario is considered, the flow structure is such that behind the shock, a small induction zone can exist which lies between the temperature rise due to the shock and the temperature rise resulting from the reaction front. If the standard form of the governing equations are solved, a coarse grid results in a computed flow field such that the shock and reaction fronts are merged and the induction zone is not resolved. Sussman and Wilson (1991) have shown that a coarse grid used in conjunction with the transformed governing equations results in an overprediction of the induction zone. In short, the computed reaction front is positioned too far down stream. Therefore, even though the transformed governing equation method has shown promise when applied to the shockinduced combustion scenario (Wilson and Sussman, 1993), this grid dependence should be considered when analyzing the results. 3.6.2 Logarithmic Interpolation The scheme of Shang et al. (1995) uses a logarithmic interpolation weighted by the mass flux into the cell and can be explained using the image of Figure 3.1. In this case the mass flux (rm) into cell (i,j) is from cell (il,j) and (i,j+l). The scheme first estimates the mass fraction at the cell faces using Log (i,j = Log a1,'j) + Log(ai)] Logi+) = Log ) + Log ) (3.39) L W = ILogij) + Lg(aij)] Figure 3.1. Schematic of mass fraction interpolation scheme per Shang et al. (1995). and then calculates the mass fraction to be used in the source term evaluation by m. i_,Log (i + mi. ,Log(t + 2 J, ij + cs = Exp +. (3.40) m. .+ mn. ? j '+ where a represents each species mass fractions. The interpolations are a sum of all cell faces for which the mass flux is into the cell. In actual implementations, Log functions are singular when the mass fractions are zero. Therefore, after simple arithmetic, the interpolated value for the example case is simply r/ \m \ 2!m 1 +m 1 \ as = 2" ilj)' T,+l) \ ++ (m i j) (3.41) Shang et al. (1995) present both solutions of a onedimensional autoignition case and the M=6.46 case captured by Lehr (1972). In the onedimensional case, the modified source term treatment produced "exact" solutions using only about a fourth of the grid points required to achieve the same accuracy with the piecewise representation. This interpolation scheme also improves the comparison between the experimental data and the computed flow 0 ij1 field for the M=6.46 case. This is accomplished with no addition of grid cells or adaptation. These results are discussed further in the context of the results of the current study. 3.6.3 Error Estimation Based on OneDimensional Model Equations The source term modification proposed by Sussman (1993) is based on the analysis of a onedimensional model problem of the governing equations. The analysis is summarized here for completeness. Consider the model problem aQ aQ + u = kQ at ax (3.42) Q(, t) = Qo and when discretized using a firstorder upwind scheme is Q+l = Q f At Qn) + Atn (3.43) where the source term is Qn = Q(Q) = kQ (3.44) Note, in Sussman's (1993) analysis, u is set to 1. The steadystate solution can be approximated by Q(x) = Qoekx/u (3.45) and the source term by 2(Q) = kQ (3.46) where k=kf. The values of the discretized function and source term, using a constant spatial discretization, are Qi = Qoek(i1)Ax/u gi = kQek(i1)Ax/u (3.47) Therefore, the approximation of the spatial derivative of Q at point i, using firstorder upwinding, is Di = Q = Qoek( 1)Ax/u[ ekAx/u] 1 (3.48) Ax/u The error present in a steadystate solution of the model problem can be estimated by analyzing the ratio of this operator to the source term, D [1 ekAx/u] QTi kAx/u (3.49) DZi kAx/u Recall that k represents the reaction rate and 1/k is the characteristic time scale of the reaction. Therefore, kAx/u is the ratio of the actual chemical time scale to the fluid dynamic time scale, associated with the grid spacing, used in the computation. This is a cell Damk6hler number, D)cfd. The variation of the ratio in Eqn. (3.49) with kAx/u is given in Figure 3.2. Both the figure and Eqn. (3.49) demonstrate that as kAx/u oo the ratio 0.9  0.8  DI/W 0.7 0.6  0.5  0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 kAx/U Figure 3.2. Discretization error for Sussman's (1993) model problem. Di/Qi is asymptotically proportional to [kAx/u]1. This relationship has direct implications to the solution of chemically reacting flows and demonstrates the expected trend. In a given flow scenario, the velocity is set by the solution as is the rate k, through the field variables and the reaction mechanism used. The advantage of adapting the grid to reduce Ax in regions of large k is obvious given the fact that the error changes exponentially with kAx/u. However, the effectiveness of such adaptation is dependent on the magnitude of k and it is clear for large k, Ax must be very close to zero to reduce the error. We next consider alternative ways of reducing the error introduced through the discretization of the exponentially varying function. Sussman (1993) uses the information from the error analysis to propose a modification to the source term evaluation using a scaling factor f such that a, = f ij. In the context of the model problem, O(Qi) = Q(fQi) = kfQi (3.50) where f is a scaling factor which gives Di/Wi = 1 if f [ ekAx/u] f = (3.51) kAx/u In a general problem f is not known a priori because k is not a constant. Sussman (1993) estimates k using the exact solution of the model problem, k= u ln( (3.52) ul(Qn) Using this estimate, which is exact if Q is truly an exponential function in space, the scaling factor when using a firstorder upwind scheme is 1r f = n(r) r = Qn /Q (3.53) and for a secondorder upwind scheme, is 3 2r + ir2 f = 2 )2 r /Qp (3.54) f ln(r) r I (3.54) As presented by Sussman (1993), to achieve the exact solution of the model problem with this modified source term treatment and a firstorder scheme requires only 1/8th the number of grid points if a piecewise treatment of the term is used with a secondorder upwind scheme. In applying this modification to a quasionedimensional shockinduced combustion case for hydrogenair, the modification is applied to only those species produced by the chainbranching reactions. For application to twodimensional flows, it is assumed the exponential variation occurs in the direction of the flow and k is estimated by k =  In(ri) ln(r) (3.55) Ax Ay J where cos 0 = u sin = (3.56) (u2 + v2)2 (u2 + v2) This results in a modification factor for use with secondorder upwind schemes of Cos1 2ri + r2? + in ( 2r + r2? f Ay (3.57) csn(ri) ln(rj) where the parameters ri and rj are determined using upwinding, ( Jilj/Qacij, u > 0; fQi,j /Qi,j, v > 0; ri = [Oi+l,j/ai,j, u < 0 =i lai,j+l/eai,j' v < 0 (3.8) A separate factor is calculated for each species a produced by the chainbranching reactions. To ensure the modification factor is bounded, the limits ri < 1 rj < 1 are imposed. Sussman (1993) demonstrated this technique for the M=6.46 case of Lehr but found the scheme did not produce results comparable to the experimental data unless an adapted grid is also used. These results are discussed in more detail later in conjunction with the current study. Given the above discussion, it is of interest to compare the modification of Shang et al. (1995) in the same context. Beginning with Eqn. (3.41), and assuming the mass fractions in the neighboring cells enter the formula through a similar representation of k in the two dimensions as that used in Eqn. (3.52), we have ekAx/u Ci ,j ekAy/v + (3.59) "ij ,i,j and after simple arithmetic, Eqn. (3.41) becomes [m kAx/u m ikAy/v 2m ( .,+m) ) , as = e TJ e + 2 where v is used to represent the velocity in the second direction. For comparison purposes, consider the onedimensional case as done by Sussman, we have Ay=0 and mn = 0. Therefore, the logarithmic extrapolation results in a scaling factor of f = e2Ax/u (3.61) which at first looks quite different from that of Eqn (3.51) but as evident in Figure 3.3, when Scaling factor 2 2.5 3 3.5 4 kAx/u Figure 3.3. Comparison of the scaling factors from the logarithmic extrapolation and error estimation. plotted for a range of kAx/u is comparable to the scaling factor by Sussman (1993) for firstorder upwind. Also presented in Figure 3.3 is the scaling factor for the secondorder upwind scheme of Eqn. (3.54). 3.6.3.1 Global Order of Accuracy Analysis The scaling factor proposed by Sussman (1993) has been derived by analyzing the local error introduced by the piecewise representation of the source term. It is also enlightening to consider the global order of accuracy and the impact of this scaling on the (3.60) magnitude of the global error. The global order of accuracy is dependent on the accumulation of error at the locations upstream of the local computational stencil (Shyy, 1994) as depicted in Figure3.4. This is particularly of interest in the current reactive flow problems since one Q(xn) Qo0 Xno Figure 3.4. Depiction of region over which error accumulation affects the global error of accuracy. of the key flow field features is the variation in mass fractions. Consider the model problem, Eqn. (3.42), used by Sussman (1993) which has an exact solution for the steadystate case of the form Q(x) = Qoea (3.62) where a = k/u. Employing a firstorder upwind difference operator with the piecewise source term representation, the difference equation used to represent the model differential equation is q(x + h) q(x) h a q(x) (3.63) where h is the grid spacing and for the current analysis constant spacing throughout the domain is assumed. The exact solution of the difference equation is qn = (1 + ah) qn1 = (1 + ah)[(l + ah)qn_2] = (1 + ah)n qo (3.64) qn = (1 + ah)x"/h oq and to estimate the global error the solution of the difference equation is subtracted from the exact solution of the differential equation E(xn) = Q(xn) qn = Qoeax (1 + ah)X"/hqo (3.65) Since the initial conditions for the difference equation are specified, qo = Q0 and the leading order of the global error is dependent of eaxn (1 + ah)x"/h. The second term can be expanded as (1 + ah)xn/h = exp ILog(l + ah) = exp [(ah 2 + (h3)) a2X h )1 (3.66) = exp axn a2h + O(h2) a2Xnh = exp(axn) exp(a2h) exp(O(h2)) where it has been assumed ah < 1 in the expansion of the logarithmic function. In implementation this may not be true but here the global order of accuracy, and specifically how the error decreases as h is reduced, is being assessed so the expansion is applicable. Using the expansion ex = 1 + x + x2/2! + x3/3! +... Eqn. (3.66) becomes (liahx)n/h i1 an2x~h (1 + ah) /h = exp 1 a2 O h2) ea2 (3.67) = e2 h (eax" + O(h2) and the global order of accuracy is E(xn) = Qo h eax" O(h2)] (3.68) This results in a normalized error, when the piecewise source term is used, of E(xn) Enorm(xn) = ( 0O(h) (3.69) Q(xn) which is an indication of the error at locations far from the boundary. This error, an accumulation of errors across the region from where the boundary conditions are imposed to the point of interest, can be key in problems such as the current shockinduced combustion flows where features such as induction zones and reaction fronts are of interest. This is due to the fact that regions such as induction zones are areas over which quantities such as radicals are growing. Since factors such as the mass fraction of radicals directly impact the resulting chemical process, any error in their growth pattern introduce inconsistences between the simulated field and what is actually experienced. When the scaling factor proposed by Sussman (1993) is used, the actual difference equation being solved is somewhat different from that of Eqn. (3.63) since the scaling factor is a function of the dependent variables. The exact form, using a firstorder upwind scheme, is q(x) q(x + h) q(x) 1 q(x+h) h = a q(x) qx) (3.70) Log(x + h) which can be simplified to qnLog( qn1Log ( ) + ah(q qn1) = 0 (3.71) and further reduced to (qn qn) Log( + ah = 0 (3.72) which has two roots, qn = qn1 qn = eah+Log(qn) (3.73) Dismissing the first root due to the nature of the problem being solved, the exact solution of the difference equation is qn = eah qn1 = 1eah[e n2 (eah) nq (3.74) qn = (eah)n/h q and when used to calculate the global error results in E(xn) = Q(xn) qn = Qoeax, q0eax = 0 (3.75) given the fact the initial conditions for the difference equation are specified so qo = Q0. Therefore, the exact solution of the difference equation being solved when the scaling factor is introduced matches the exact solution of the original differential equation. This infers the accumulation error when the scaling factor is included is minimized and is further reason to use the scaling factor in the shockinduced combustion computations to be performed here. 3.6.3.2 Current Application of Scaling Factor Scheme In the current study the two scaling methods described above are compared by applying each to the solution of the shockinduced combustion flow. It is obvious by this point that the error is introduced in the solution of the problem through the representation of the spatial derivative. This is due to the fact that the mass fractions of the species in many cases vary exponentially over a region much smaller than the typical grid spacing. Therefore, the source term for a given cell calculated based on the cell average of the dependent variables is not equal to the cell average of the true source term. Given the use of the scaling schemes, it is worth assessing how close to the governing equations is the model problem used by Sussman (1993). The species conservation equation can be reduced to __ + U = c (3.76) at ax using the continuity equation. In the actual governing equations the source term 6) is a function of all species involved in the reactions. Therefore the factor k in Sussman's model problem (Sussman, 1993) is an effective rate of production of a per concentration of a, or the derivative of the source term with respect to the species ac. This information is available in the form of the source term Jacobian, given in the Appendix. Alternate estimations for k may prove useful but here we use the form as proposed by Sussman (1993) to facilitate comparison to these earlier results. In the documented scheme by Sussman (1993), the criterion for selecting the upwind direction in specifying the ratio of mass fractions, Eqn. (3.58), uses the Cartesian velocity components. It is not clear if this has been done simply for notation but would be very erroneous if Cartesian velocities are used in conjunction with a curvilinear grid. Therefore, here the contravariant velocities are used to determine the upwind directions. Also, the particular form of Eqn. (3.57) used by Sussman (1993) in the twodimensional computations determines the error estimate using a weighted sum from the two curvilinear directions. However, the Cartesian velocities are used in the weighting which may introduce some inconsistencies. Therefore, the current study evaluates the scaling by Shang et al. (1995) and Sussman (1993), with curvilinear upwinding, and also evaluates the scaling factor in the form mw 2ri + r) + r n 2rj + r2) rTw ln(ri) finnln(rj) (3.77) 2,j i~j + mw = mn = S.2 2 2 2 m_,j +Im m + m + ij + ij+ where the components are as defined earlier and the flow scenario is that of Figure 3.1. The schemes of Shang et al. (1995), Sussman (1993) and the current scaling factor, Eqn. (3.77), based on Sussman's work (1993) are evaluated here and the results are presented later. When these scaling are used, the effective reaction rate is not just the change in the mass fraction over the distance between neighboring cells but a velocity scaling is included as detailed in the derivation using the model problem. As is evident by the attention to both adaptive gridding techniques, transformed equations, and modified source term treatments, the additional scales introduced by the chemical reaction processes in the current problem set can not be ignored. Various studies have evaluated these three options for treating the problem and here the modified source term treatment is tested. This is done in part because if such a treatment proves worthwhile, it is easily integrated into the solution process. The scaling factor offers an attractive option to grid adaptation but the error estimation required is only as good as the estimate of the effective rate k. Given the behavior of the error as kAx/U increases, any errors in this estimate are heightened for large k. 3.6.4 Nonphysical Wave Speeds LeVeque and Yee (1990) also document the problems associated with highly nonlinear source terms. In their study a model problem is used to study the behavior of the linear advection equation with a source term. Nonphysical wave speeds are noted for large At. The source of the nonphysical speeds is found to be solely a function of the presence of the source term since the conservative numerical methods used correctly reproduce the wave speed if no source term is present. LeVeque and Yee's (1990) analysis is summarized here. Consider the governing equation of the current reactive flow problems. If no source term is present, the numerical scheme is conservative if the summation over the domain results in the cancellation of the defined fluxes. In a word, for a finite volume scheme, the defined fluxes at the cell faces must balance. The addition of a source term results in a quantity that does not cancel when summed over the domain and therefore, if this term is incorrectly modeled, the accuracy in characteristics such as wave speed are reduced. In the numerical model of both the model problem studied by LeVeque and Yee (1990) and the full set of governing equations to be solved for reacting flows, the source term is typically represented as a function of the cell averaged a, w(Ui). In the context of the integral form of the governing equations which is used to validate conservatism, this replaces the average value of the source term, ii, with a source term which is calculated base on an average a. The discrepancy between these two is small for a smoothly varying a but in the case of reacting flows where sharp variations exist, the difference can be large. The conclusions of LeVeque and Yee (1990) on how to treat the occurrence are much the same as those already presented. First is the option of refining the grid which is not appealing as mentioned earlier. A second option proposed by LeVeque and Yee (1990) is to better represent the integral of the source term across a given cell. This requires a better model for the variation across the cell which is in the same line as the logarithmic variation approximations evaluated here. 3.7 Numerical Solution Methods Many options exist for the integration and numerical solution of the governing equations. These options typically fall in the categories of explicit or implicit schemes. Problems including chemical reaction source terms as here normally require an implicit method for the treatment of the stiff source terms. This can be achieved in two ways. The first option is to solve the governing equations in a fully coupled manner which requires the derivatives of both the fluid dynamic and source terms with respect to the dependent variable vector, Q. These derivatives, Jacobian matrices, can be quite large in the case of a multispecies flow. The second option is a splitting method in which the fluid dynamic terms are treated explicitly and only the source terms are treated in an implicit manner. Both of these schemes are to be used in the current study so a description of each is provided here. 3.7.1 FullyImplicit Formulation As detailed in the description of the flux scheme, a secondorder spatial accurate model is used. Likewise, it is advantageous to use a secondorder temporal scheme as well. A detailed review of the fundamentals of achieving such a scheme is available in Hirsch (1990a, 1990b) and the current implementation follows that of Shuen (1992). A concise review is provided here. The governing equation for the current problem, Eqn. (3.2), is of the form dQ d = M(Q, t) (3.78) where M is a nonlinear function of Q and represents the differential space operators and source terms of Eqn. (3.2). A onestep scheme for solving the system of equations is Qn+l Qn = 8t [{Mn+' + (1 )Mn] (3.79) and when 1 = 1/2 the scheme is second order in time. Since the whole purpose of solving the set of equations is to determine Qn+ the question is first how to represent Mn +1 = M(Qn +1). One option is to use a predictor and corrector approach which results in an explicit scheme. The second option is to linearize the nonlinear term M +1 which results in a fully implicit scheme. The term is linearized as Mn+ = Mn + ) t + O(6t2) QJ\ n 6t = Mn + at) at + 0(8t2) (3.80) = M + Nn 6Q" + O(t2) where 6Qn = Qn+1 Qn and the Jacobian, N = aM/aQ, has been introduced. Introducing the linearized terms the equation becomes 6Q = 6t [~(Mn + Nn"Q) + (1 )M"] = 6t [(Nn8Q) + Mn] (3.81) Using At = 6t/J and the definitions F = J(F, + Fv) G = J(G + G) (3.82) (3.82) H = J(H, + Hv) S = JQ the governing equation, Eqn. (3.2) becomes I + 'At(8A + B + H' S')1 AQ= At (tF + G+ HS) (3.83) where 6N and 68 represent the spatial discretization in the two curvilinear directions, A = aF/aQ, B = 8G/cQ, H' = aH/aQ, and S' = aS/aQ. Note, has been set to 1/2 for second order accuracy. When implemented, the Jacobians of the viscous terms are neglected which has proven to save computational time with no negative impact on the solution process. Introducing the flux vector splitting described previously, Eqn. (3.83) can be symbolically expressed as [D +L + U] AQ = At R (3.84) where the components of the implicit operator are D1 + (A+t A +B+ B +H' ) D = I + IAt A j A j + B j Bi,j + Hi,j S i,j L At (Ai1j+ Bij1) (3.85) U = At Ai+l.j + Bi,j+l) and the explicit operator is R = At (Fi+ F.,. + G. ,L G.. + Hi )Si (3.86) Due to the linearization, all the components of the implicit and explicit operator are evaluated at time level n. As detailed by Janus (1989), when solving unsteady flows, many times subiterations are used in the solution process to minimize the linearization error. The Jacobian matrices of the split fluxes are exactly A = 8F: /Q but are approximated here as in previous studies (Shuen, 1992) using A: = TAT' (3.87) The matrix A* is a diagonal matrix composed of the eigenvalues X* given in Eqn. (3.23). The matrix T is composed of the right eigenvectors of the matrix A, in columns, and T the corresponding left eignevectors in rows. The eigenvectors are provided in the Appendix. The exact Jacobians are available in Janus et al. (1992). Here the system of equations represented in Eqn. (3.84) is solved using a lowerupper (LU) approximate factorization scheme described by Shuen (1992) and Shuen and Yoon (1989) and based on the scheme developed by Jameson and Yoon (1987). The lefthand side of Eqn. (3.84) is approximately factored into the product of two operators [D + L]D'[D + U] AQ = At R (3.88) and when implemented can be solved in the sequence [D + L]AQ* = At R, [D + L]AQ = DAQ* (3.89) Qn+1 = Q" + AQ Various stability advantages of this solution scheme are detailed by Shuen (1992). 3.7.2 ExplicitImplicit Splitting Scheme As mentioned, options for the solution method for the current problem range from fully implicit to fully explicit schemes. Implicit schemes are essential when solving problems with strong source terms as is the case here. However, explicit methods are attractive due to their reduced computational requirements. Also, there is some question as to the advantage of fully implicit schemes for this class of problems. For highspeed reacting flows, small temporal steps are used to either capture the unsteady nature of the flow or in light of the small time scales of the chemical processes. The need for an implicit scheme for the treatment of the source terms and the efficiency of an explicit formulation can be met using the Strang splitting method described by Le Veque and Yee (1990). The scheme computes the dependent variables (Q) at the new time level (n+l) by performing various operations on the variables at the current time level (n). The splitting can be represented by Qn+I = SQ(At/2) SF(At) Su(At/2) Qn (3.90) where SQ represents the numerical solution operator for the source term and SF the numerical solution operator for the fluid dynamics conservation equations. The Strang splitting has two features which make it very attractive for the current problem. First the splitting maintains secondorder accuracy which is needed to reduce the numerical dissipation of the important phenomena. The second feature is the fact that the splitting allows proven solution schemes for each aspect, the chemistry and fluid dynamics, to be easily integrated. Such a splitting may raise some concerns about capturing the coupling between the fluid dynamics and the chemical processes. However, as demonstrated in the combustion instability simulation work of Shyy et al (1993), this approach yields a very satisfactory treatment for the source terms. As mentioned, in the splitting technique the fluid dynamics aspects of the problem are modeled using explicit schemes. To maintain secondorder accuracy, the fluid dynamics operator must be secondorder and here a predictioncorrection scheme is used in the form given by Hirsch (1990b). To clarify the steps in the solution process, the sequence of equations are as follows, S,(At/2) : I Q' AQ* =At Q* = Q + AQ* Sp(At) : AQ" = 8[aF + aG)' + H] Q** = Q* + AQ* (3.91) AQ*** = At[aF(2) + G(2) + H] Q*** = Q* + AQ*** Sn(At/2) : I Q' AQ****= At*** Qn+l = Q*** +AQ** where the superscript *'s denote what value of Q is used to calculate the fluxes, source terms, and Jacobian of the source term, Q'(Q"i,j) = aQ(Qnij)/Qni, j. The superscripts (1) and (2) denote the spatial order of the numerical fluxes. The term a F = Fi+. F. and S 1+ .1 1 2'j likewise for 9,G. The scheme used to define the inviscid numerical fluxes is the same as that used in the fullyimplicit method. However, here the Jacobians of the fluid dynamic terms do not have to be calculated which reduces computational time as does the fact that a banded system of equations does not have to be solved. 3.7.3 ODE Solution Issues Even when implicit integration is used difficulties can arise when implementing the solution process for problems with large reaction rates. Here some of these issues are discussed in the context of solving the ODE associated with the source term operator in the splitting solution method. One difficulty is due to the fact that the Jacobian matrix Q' can become illconditioned which can introduce errors when the equations are solved. The Jacobian matrix can be changed somewhat to help alleviate this problem. Consider the two step reaction H2 + O2 20H (1) H2 + 20H 3 2H20 (2) (3.92) which is used later in this study and has been presented in Table 2.1. The source term for each species, when expanded, is )H = MHj kf,(XH)(Xo + kb(XoH)2 kXH)(XOH)2 + kb XHo)2 o, = Mo{ kf, (X,)(Xo) + kb(XoH)2 OH = 2MoH kf,(XH2)(Xo) k OH kb(XH(XoOH)) + kb(XHo)} (3.93) HO = MH20 2k,(2XH)(XoH)2 2kb2(XHo)2 and after some simple arithmetic those for H2 and OH can be written in the form MH6. MH [b0 2MOH [o2] OH (3.94) 2"OH = M 2[ H20 02 HO The Jacobian matrix of the source term can be altered to reflect the relationships of Eqn. (3.94). This implementation has been used for reaction mechanism composed of two and eight steps and has proven to increase the robustness of the computational code in the case of the twostep mechanism. This is due to the fact that the reaction rate constants for the twostep mechanism are orders of magnitude larger than those for the higher order mechanisms and result in a Jacobian matrix with a much larger condition number. For the higher order mechanisms, the Jacobian matrix using the standard form is well conditioned. Even if the Jacobian matrix is well conditioned and a solution is obtained, the resulting state of Q can be nonphysical. This is typically manifested in the mass fractions of various species being less than zero. These nonphysical solutions are a result, once again, of the fact that quantities such as mass fraction vary in a logarithmic fashion. Because the solution process utilizes a linearization for the derivative of Q, these overshoots are produced. In the solution of problems such as the reacting flows currently studied, such nonphysical solutions must be guarded against. This can be done by either taking smaller time steps, which is not desirable, or by employing a scaling technique such as the method demonstrated by Ramshaw and Chang (1995). Here a technique is used which scales AQ if the solution of the ODE in Eqn (3.91) would result in a nonphysical Q. In the current code, the solutions from the ODEs of Eqn. (3.91) are evaluated to see if Q* or Qn+ 1 contain species mass fractions such that ai < 0 for all i. If so, a single appropriate scaling factor is calculated and applied to all AQ's such that 0 < ai < 1 for all i. The change in Q due to the fluid dynamic operator presents no problem. This is due to the fact that the At is calculated based on the CFL < 1 condition necessary for explicit schemes. When applying such a scaling, caution must be used to ensure the actual equation being solved is not altered. In the case of steadystate solutions, typically the scaling is needed only in the early stages of the evolution of the flow field. If the scaling factor converges to 1 long before the convergence of the governing equations to the steadystate solution, then there are no concerns that the equations have been altered. In the case of unsteady problems, an iterative process can be used to solve the ODE to ensure nonphysical solutions are not encountered. However, typically when solving unsteady problems, the time step is set small to resolve the evolution of the flow field and such problems are not encountered. It should be noted that these overshoots are different than the roundoff errors which can appear as mentioned earlier and as addressed in previous investigations (Sussman, 1993 and Wilson and Sussman, 1993). 3.8 Summary Here the various models needed to successfully simulate the problems of interest have been described. As is clear, there are numerous phenomenological and numerical models available. Those selected to be used in the current development offer the latest in methodology and have a wide range of applicability. This is important in the context of the current study given the fact that a multitude of scales are intertwined. Those simplification that have been made have been done so after analyzing the phenomenology of the process and based on previous investigations. From the numerical aspect, one of the key items to be highlighted is the source term treatment. This is due to the fact that previous investigations have shown the accuracy is highly dependent on the grid resolution. The source term treatments offer a viable alternative to the computationally intensive grid adaptation schemes. The primary chemical aspect to be focused on is the impact of the reaction mechanism. The model described here is used to compute the shockinduced combustion problem with various mechanisms in a effort to correlate mechanism characteristics to simulation accuracy. The impact of turbulence in the reactive flow scenarios presented here is also quantified. Both the modeling of the turbulence and the coupling of these effects with the chemical processes are detailed in the next chapter. The turbulence effects are known to be manifested in the governing equations through the viscous aspects. It is clear from the description of the governing equations that these aspects impact field quantities but also the chemical process through the dispersion of species and heat. CHAPTER 4 TURBULENCE MODEL The turbulence aspects of the problems of interest like the other physical and chemical aspects require a phenomenological model. As eluded to earlier, turbulence impacts many of the aspects already discussed as well as the particular models that have been described. Therefore, this aspect of the problem is dedicated its own chapter to emphasize both the phenomenological importance of turbulence but also the attention which must be paid to its modeling. A turbulent flow, like all continuum fluid problems, is governed by the NavierStokes equations. But to model these equations on the length and time scales to resolve all pertinent physical phenomena is not a practical method given today's computational platforms. Simulations which do resolve all the scales ranging from the small to the large eddies are conducted and are labeled Direct Numerical Simulations (DNS). A review of DNS is provided by Galperin and Orszag (1993). To highlight the computational resources required to conduct a DNS let us estimate the number of grid points that would be required for the fuel cell scenario. As mentioned, the phenomena of turbulence is known to be driven by a multitude of scales, particularly the large scale eddies which contain what is call the turbulent kinetic energy and the small scale eddies which work to dissipate this energy. The dynamics of turbulence consist of various transformation processes which transfer this energy between the eddies. To capture all the important phenomena, the governing equations must be discretized spatially such that all scales are resolved. This infers that the computational grid used must contain points spaced on the order of the size of the small eddies. The characteristic length scale of the small eddies is the Kolmogoroff scale (1) and the ratio of this scale to the characteristic length scale of the large eddies (le) is related to the large eddy Reynolds number by Yl/le = (Re ) where Rel = . Therefore, to capture all scales of interest, the grid spacing must be on the order of I] = le(v3/u31e3 ). Given the fact that these small eddies can have an arbitrary orientation in the flow, this same spacing must be used in all three physical dimensions. An estimate of le is needed and for the fuel cell scenario, the large scale eddies are on the order of the thickness of the various shear layers which mix the fuel and oxidizer. These shear layers are on the order of the size of the projectile so this length scale can be estimated as le = O(d) where d is the projectile diameter. Using the Reynolds number of the mean flow (Re = ud/v) we can estimate the Kolmogoroff scale to be 1 = 0 d(Re) ]. Recalling that the discretization must be on this order in all directions and considering the computational domain, at the least, must be on the order of the projectile diameter, the number of points to carry out a DNS study of this problem is on the order of 9 (Re)4. Using quantities for the viscosity, length, and velocities comparable to those in the problems to be solved here, we have Re = O(1.0x107) which would require on the order of 5x1015 spatial points on which the governing equations must be solved. This is clearly not practical for typical engineering analysis or much less easy given the current computational platforms. However, DNS is pursued for moderate Reynold's numbers and simple geometries in hopes of better understanding the dynamics of turbulence thereby advancing optional modeling techniques. These optional techniques are required to complement the solved set of governing equations in an attempt to reproduce the dynamics of those scales known to not be adequately resolved. The temporal scale of the small eddies must also be considered if a DNS is to be performed. This time scale can be thought of as the turnover time of the small eddies. Like the spatial requirements impact the number of points at which the governing equations must be solved and in turn the amount of computer memory needed, this requirement on the temporal discretization effects the required computational time as does the small time scales of the chemical reactions. The next level from DNS in the context of grid requirements is Large Eddy Simulation (LES) in which grid refinements are used such that the large scales of turbulence are resolved but a model is needed for the small scale phenomena. LES is also characterized by the use of spatial filters on the governing equations. This terminology is used because the computational grid has to be constructed such that the spatial discretization is on the order of the inertial subrange, or that region which lies between the large and small eddies and over which the primary physical process is the transfer of energy. Since the primary role of the small eddies is known to be one of energy dissipation, this is the phenomena which must be modeled. A review of LES is provided by Erlebacher et al. (1992) and Germano (1992). In the context of solving a full scale engineering problem, the term large eddy is somewhat misleading in that the size of the large scales in turbulence can still be orders of magnitude smaller than the domain over which the solution must be carried out. Therefore, the LES requirement to resolve only the large scales is still computationally restrictive. However, LES is currently being pursued aggressively and has been used for many problems of interest. LES will continue to become more viable as computational platforms become faster and cheaper. Based on the description of DNS and LES, it is clear that an optional modeling technique is needed for problems composed of large Reynolds numbers and large computational domains. The problem in modeling turbulent flows is that the highly nonlinear dynamics, represented by the NavierStokes equations, occurring at these unresolved scales produce phenomena when observed from a larger frame of reference appear random and nonsystematic. The model of choice for problems such as the reactive munition and shockinduced combustion scenarios remains a class of engineering models labeled moment methods which are solved in conjunction with a set of averaged governing equations. There are many options as to both how to average the equations as well as the form of the engineering models. More details on various engineering models is provided by Wilcox (1992) and a review in the context of reacting flows is provided by Libby and Williams (1994). Here the work of Krishnamurty (1996) is followed to describe the implementation of one particular model in an attempt to quantify the role of turbulence in the reacting flow problems of interest. 4.1 Averaging Techniques and Averaged Variables As mentioned, when a turbulent flow is observed, an apparent randomness exists. Conceding the fact that all scales are not to be resolved, it is assumed that a solution on the resolvable scales can provide a statistical average of the exact flow field. In LES a spatial averaging, or filtering, procedure is used. Here, however, the averaging is cast as either a time or ensemble average. A time average of a quantity with instantaneous value (d is to +T = lim j dt (4.1) T.00 T where T is the period over which the instantaneous values are sampled. This results in a formula for the instantaneous value at a given time of the form + = + (4.2) where (' is the deviation from the mean and by definition 4' = 0. When solving unsteady flows, it is more applicable to employ a methodology of ensemble averaging. The ensemble average of a quantity c( is an average over a number of realizations of the problem and mathematically is N = lim Z n (4.3) Nco N n=1 where again (' is the deviation from the mean and by definition 0' = 0. A realization is simply the same problem restarted from the same initial conditions with the quantity 4 being measured at the same point in space and at the same time after the start of the problem. The mental picture of this type of averaging is an excellent way to frame the problem of turbulence. Consider the task of collecting an ensemble average of say the velocity at a specific point in the flow field about a model in a wind tunnel at a specific time after the tunnel is started. Assuming the measurement can be made precisely at the same point in time and space, the variation from measurement to measurement is a function of the nonlinearity of the governing equations. Small perturbations can be generated in any quantity whether it is caused by the startup process, a slight movement of the model, a small contaminant on the tunnel wall, or a miniscule change in the environment. These small changes can grow within the framework of the NavierStokes equations and effect the field quantities at various points in space and time. And as mentioned several times now, to capture this evolution with a computational model is no easy task. Continuing with the premise that when the NavierStokes equations are solved on a grid known to introduce some statistical averaging, the equations which are actually solved need to be cast into a form accounting for this averaging process. The first step in doing so is defining the averaged variables based on the averaging techniques discussed earlier. When solving turbulent flows, averaged variables in the form of Eqn. (4.2) are typically defined as Reynolds averaged variables. With the added desire to model compressible flows, a second averaged variable definition is introduced. This is in light of the fact that unlike incompressible flows, the added aspect of compressibility can introduce deviations, or fluctuations, in density. By introducing a mass weighted average, denoted as a Favre average, the resulting equations are somewhat simpler. Using this form, the instantaneous value is given by 4)= 4) + )" (4.4) where S(4.5) which results in the relationships = = 0. (4.6) Note, the tilde and double primes are used to distinguish Favre averaged variables from Reynolds averaged ones. One reason this approach is taken in modeling compressible turbulent flows is that with the density weighting introduced, the final form of the governing equations is closely analogous to the incompressible form and it is in the context of incompressible turbulence that most models have been developed. Following the typical convention, a mix of Reynolds and Favre averaged variables are used. Specifically, Reynolds averaging is used for the density and pressure and Favre averaging for the velocity components and scalars such as temperature and species mass fractions. The forms are Q = +Q' p P + p' ui = uf + u," (4.7) T = T+ T" ai = di + a"ill The next section presents the form of the governing equations once these averaged variables are introduced. 4.2 Averaged Equations A review of the averaging process is provided by Wilcox (1992). Here the final form of the equations are presented. This final form is obtained by introducing the defined average variables of Eqn. (4.7) into the governing equations of Chapter 3 and making use of the relationships given above. Note, for convenience, the Cartesian form of the equations with tensor notation is used in the discussion. The continuity equation takes the form i() + )= 0 (4.8) and the momentum (ugiJ) + (ax u) = + x i u"u ") (4.9) t i xi axj ij where i v vi' + i I +X1 X 6 axx axi axk (4.10) [au." u." [ "1 (4.10) Tj"=v a + + + x 6,uk" j ax ax " ixk and the term u, "uj" is yet to be defined. The energy equation, with the introduction of total enthalpy, H = E + and specific enthalpy, h = e + becomes )+ ax( ) = "i + Tij Qu"uj") + ui"ti + ui Tij (4.11) qj quj"h" Quj'(,u,"u,") Recall for reacting flows, the heat flux is also a function of species diffusion velocities as given in Eqn. (3.11). Also, as a result of the averaging process, the the total energy, E, now includes the turbulent kinetic energy which is defined later. The final equation is the species conservation equation and when averaged for a single species n is (n) + nj c ax ean"uj = (4.12) where the source term is denoted as an average and the particular representation is discussed later. The first terms to be addressed are those that are dependent on the product of the molecular viscosity and the fluctuating terms. This includes Tij" in the momentum equation, ui "ij and qj" in the energy equation, and Sc I in the species conservation equation. These terms are small in magnitude when compared to the terms that are a function of the product between molecular viscosity and the mean field quantities. Therefore, these terms are ignored in the modeling. The remaining terms are significant and must be modeled. 4.3 Turbulence Closure The first term in need of modeling, or closure, is Qu "uj" which is a symmetric tensor with six independent components. Analysis shows that this term behaves much like the molecular stress term and is labeled the Reynoldsstress tensor. Boussinesq first proposed modeling this term in an analogous fashion, dependent on the gradients of the mean field. This requires an effective viscosity dependent on the turbulence, and the estimation of this turbulent viscosity varies based on the model employed. The first level in sophistication is the mixinglength model proposed by Prandtl (Tennekes and Lumley, 1972) which uses the gradient of the mean velocity field and a yet undetermined mixing length to define the viscosity. In a shear layer flow, the form for the turbulent viscosity is tt = Qmix (4.13) where u is the streamwise component of velocity and y is normal to the flow direction. This model performs well for some flows in which the mixing length, Imix, has been empirically determined. However, this model offers limited potential for complex flow fields in which a single length scale does not hold. A more appropriate model for such flows is one which incorporates the dynamics of turbulence in the estimation of turbulent viscosity. This requires the definition of a velocity scale and length scale. Various models have been used and a review is provided by Wilcox (1992). Here, one particular model is discussed and integrated with the other models required for the highspeed reactive flow problems. The model used here is the ke model, a twoequation model with its foundation in the work of Jones and Launder (1973). The exact equations used to govern these quantities are presented next but first the closure of the various terms in the averaged equations are discussed. The quantity k is the turbulent kinetic energy mentioned earlier and e is the dissipation of k and both are defined as k ui"u," 3ui," u,"' k = E = v aX (4.14) Q axi ax Using k and E to define the velocity and length scale of turbulence, the turbulent viscosity becomes kt = Cg (4.15) where C. is a constant found to be .09. The Reynolds stress term is modeled as r[aiu. 8a'xi [ik 1 Quu = t + axL k ij 2 (4.16) with the last term introduced to maintain the proper trace of the Reynolds stress tensor. The turbulent bulk viscosity is represented as is the molecular, Xt =  t. The next term in need of closure is Quj"h" which is modeled, using the Reynolds' analogy (Hinze, 1959), as CpIt aT( Quj"h" Prt axj (4.17) where Prt is the turbulent Prandtl number. Using this closure, a single turbulent Prandtl number, assumes the turbulent transport of heat and momentum are dependent on a single time scale, i.e. the ratio between the turbulent heat and momentum flux are the same throughout the flow field. This assumption is known to not hold in many flow scenarios and an optional modeling concept such as the work of Nagano and Shimada (1996) is to define a second velocity and length scale associated with the heat flux. This allows the thermal eddy diffusivity to be calculated explicitly throughout the flowfield. For the current study which is an initial investigation into the coupling of the turbulence and chemical reacting aspects of the flow, the closure of Eqn. (4.17) is deemed acceptable. In most flows, the turbulent Prandtl number is close to unity and that is the value used in the current study. The remaining term in the energy equation, Eqn. (4.11), is ui'"Tij" Quj"(ui"'ui") which represents the diffusion of energy due to turbulent fluctuations and is modeled as U7 t ak "" "(2 u,"u") = k (4.18) Ui Ci Qj u kaxj where ok is a Prandtl number associated with the turbulent kinetic energy and is set to unity. The one term in the species conservation equation needing closure is QCn"uj". Using the gradient treatment already introduced for the turbulent heat flux, the term is modeled as Qcg"u"  It aan (4.19) Sct axj where the turbulent Schmidt number, Sct, has been introduced. As in the case of the turbulent Prandtl number, assuming a single turbulent Schmidt number prescribes that the scales associated with the turbulent transport of species are the same throughout the flow relative to those of the momentum. Also, similar to the modeling option for the turbulent heat flux, additional equations can be solved for the correlations of Eqn. (4.19) but an additional equation is needed for each species which can be computationally expensive. Here, Sct=1 which also implies the turbulent Lewis number is 1. All the modeled terms are now expressed as a function of k and E so the next issue is the modeling of these two quantities. 4.4 kE Model The equation for k is derived by manipulating the momentum equation. A description is given in Krishnamurty (1996) but here only the final form of the equation is given, = p u1  (xk) + \li)jk . Quif"uj" ax uap aui t axI ax i ax axu " _,, (4.20) + a ui "ij Qu tu"ui"i) p'j] + p aX where some terms needing closure have already been defined. The definitions for the production of turbulent kinetic energy and the rate of dissipation of this same energy are such that ii Pk = Qui uj a uax (2 ^E = ~axp" Note, the term p'uj" is not explicitly modeled and is included in the closure of Eqn. (4.18). Bu." ap The terms p' and u" are, respectively, the pressure dilatation term and the enthalpic production term. A model has been proposed by Krishnamurty (1996) for these terms and has been evaluated for nonreacting blunt body flows with shocks. The work of Krishnamurty (1996) shows these terms may prove important in flows where compressibility effects are significant. The shockinduced combustion flow studied here is a highspeed, compressible flow and it may prove necessary to model these terms. However, for this initial investigation, the standard form of the kE model is used. This is due in part to establish a benchmark using the baseline model but also in light of the fact that the derivation of the models proposed by Krishnamurty (1996) need to be investigated in the context of multispecies reacting flows. There have also been several modifications proposed for addressing known compressibility effects and how they impact factors such as the dissipation rate of k. These modifications are also discussed in detail by Krishnamurty (1996). Here, the final form of the k equation used is S(k)t +a (ijk) = Pk +E + V ( + ~) (4.22) The governing equation for E is not as cleanly derived from the NavierStokes equations but is based on the assumption the time rate of change and the convection of the dissipation (E) is balanced by a production, dissipation and diffusion term in the form a + ijE) = PE E + Te + B (4.23) Unlike the k equation which is derived from manipulation of the NavierStokes equations, the terms in the E equation have been derived based on dimensional analysis and physical reasoning. The terms in the standard model are PE = E, = E2 Te t + (4.24) where Oj is a Prandtl number associated with the E and is set to 1.3. The coefficients in the production and dissipation terms can be modified to account for the imbalance between the production and dissipation of turbulent kinetic energy. These modeling option are discussed by Chen and Kim (1987) for C.1 and Thakur et al. (1996) for CE2. Here, the standard values of CE1 = 1.44 C,2 = 1.92 are used. The term BE represents the baroclinic torque term which arises due to gradients in density and pressure, () x (Vp) which do not appear in incompressible problems. A model for this term, which may prove important to flows such as the current shockinduced combustion scenario, is proposed by Krishnamurty (1996). Here, though the standard form of the equation, a + a ) = (Ck p CE2) + J + (4.25) is solved. Now that the equations for k and E have been presented, the next section specifies the needed boundary conditions. 4.5 Wall Function The kE model used here requires a wall treatment due to the low Reynolds number associated with this region. Here, the standard wall function (Wilcox, 1992) is used. This function is based on the known correlation between the nondimensional velocity near the wall and the nondimensional distance from the wall. This velocity is specified as U + = Up/UT where UT = wal/ and the nondimensional height is y + = (QUT yp)/t where yp is the distance normal to the wall and [t is the molecular viscosity. In the viscous sublayer, the region nearest the wall, the relationship U+ = y+ (4.26) holds. Outside this region, the relationship between the two nondimensional quantities is U = ln(y+) + B (4.27) where x is the von Kdrmdn constant set to .41 and B=5.1. Therefore, this region is labeled the log layer. Where the curves defined by Eqn. (4.26) and Eqn. (4.27) meet is defined yc, the critical value of the nondimensional height. In a finite volume scheme, such as the one developed here, the flow properties are defined at the cell centers. Therefore, when implementing the wall function, the y + corresponding to each cell which bounds the wall must be calculated based on the halfheight of the cell. Using this value along with velocity and viscosity, the boundary values for k and E are calculated, if y + < y+, using (Krishnamurty, 1996) U2 2 1 3 wall = wall wa (4.28) c+ w where l iY cp (C R Q wail YP S' ci = x Re = (4.29) 1 + I, It Re, If y + > yc+ then the boundary conditions are specified as U2 U3 wall Ewall = (4.30) Now that the equations for k and E have been specified as well as the boundary conditions for both quantities at the wall, the final information needed is the boundary condition for k and E at the inflow. This is typically specified in terms of a turbulent intensity, the magnitude of the velocity fluctuations relative to the mean value. Here a value of 3% is used base on data for similar flow scenarios used in earlier studies (Krishnamurty, 1996). The governing equations for k and E are solved in conjunction with the NavierStokes equations. These equations provide the necessary information with which to define pt. With the addition of the turbulent viscosity, an estimate of the enhanced diffusion of species and heat is obtained. This alone can impact the mean flow field and the impact this has on the shockinduced combustion scenario is studied. In addition to these effects, the turbulence aspects of the flow can also directly impact the reaction source term. This is the focus of the next section. 4.6 Turbulence Effects On The Reaction Source Term The chemical reaction source term, On, in Eqn. (4.12) is a function of temperature and species mass fractions. The chemical source term is a function of the mean values of these quantities and any deviations. The effect of turbulence on the source term is first represented by modeling the impact turbulence has on the mean temperature and mass fraction field. This is done in the framework of the twoequation type model by applying the averaging and closures described above. The next step in modeling the impact of turbulence on the source term is to account for the effect of the fluctuations. This is needed given the highly nonlinear dependence of the source term on the dependent variables. The analysis of Sagara and Tsug6 (1970) shows that the fluctuation in velocity does not effect the mean reaction rate, the effect of density fluctuations is relatively minor, and the primary effect of turbulence on the source term is manifested through the temperature fluctuations. They find the mean reaction rate in a turbulent field to be o = (1 + A)~exp + (1 A)Texp f l t) (4.31) where A = 6T/T and F = E/RT. The term 6T = T'T" is the rootmeansquare of the temperature fluctuation. This effective source term can be written in the form S= [(1 + A)w(T + 6T) + (1 A) (T 6T) (4.32) after simple arithmetic. Given the highly nonlinear nature of o, the resulting reaction source term can be quite larger than simply the order of magnitude of the temperature fluctuation, especially for large activation energies. This model has been used by Mizobuchi et al. (1997) to simulate a scramjet combustor scenario. In their computations the factors (1 A)2 are neglected but here, these are included. Their results show an improvement in the computed combustion efficiency when compared with experimental data. To this point in the modeling of the various terms in the averaged equations, the presence of temperature fluctuations has been solely in the context of the correlation of these fluctuations with those of velocity. As mentioned in the context of the closure of Eqn. (4.17), one modeling option is to solve additional equations for the temperature fluctuations scales. If such a model is used then the estimate of the rootmeansquare of the temperature deviation is explicitly provided. Since we have chosen to specify a single turbulent Prandtl number, these deviations are to be estimated based on the deviations in the velocity field. This is considered appropriate given the goal of the current study is to determine if any estimate for the effect of turbulence on the temperature field impacts the computed flow field in the problems of interest. Using the strong Reynolds analogy (Gaviglio et al., 1977), the ratio of the rootmeansquare of the temperature fluctuations to the mean value is related to the velocity field as 8T 2 S= (y 1)M2 (4.33) where U denotes the velocity magnitude and k is used to estimate the rootmeansquare of the velocity fluctuations. 4.7 Summary The model employed to represent the turbulent aspects of the flows has been described. The twoequation model has been selected due to its utility in solving fullscale engineering problems. The various modifications to the kE for compressibility effects and the imbalance between the production and dissipation of the turbulent kinetic energy have been presented. In the initial computations, the standard form of the equations is used to produce a benchmark to which further computations using these modifications can be compared. Turbulence is known to enhance the dispersion of species and heat and these aspects in turn affect the chemical processes. An additional effect explicitly accounted for in the present effort is the impact of temperature fluctuations. This is incorporated through 87 a modified source term dependent on the rootmeansquare of the temperature fluctuations. The impact turbulence has on the problems of interest is detailed in the context of results using these models presented in the next chapters. CHAPTER 5 SHOCKINDUCED COMBUSTION RESULTS AND DISCUSSIONS Through the review of the earlier computational studies into the shockinduced combustion scenario and model development discussion, it is clear that many issues related to both the numerical and phenomenological aspects of the models used can affect the outcome of the simulation. For example, it has been quantified in the discussion of the source term treatments that the numerical error due to the presence of the reaction source term is dependent on the cell Damkohler number, D)cfd = kAx/u. Using this parameter to frame the computational investigation into highspeed reacting flows, the phenomenological aspects of importance is the reaction rate, k, and the local flow speed, u, while the numerical aspect is the mesh spacing, Ax. To date, the numerical aspect has been the primary focus of investigations into the shockinduced combustion flow scenario. In the context of the shockinduced combustion problem, various techniques for addressing this numerical aspect have been evaluated including grid adaptation (Wilson and MacCormack, 1990, Sussman, 1993, Shang et al., 1995), source term treatments (Sussman, 1993, Shang et al., 1995), and solving a logarithmic form of the governing equations (Sussman and Wilson, 1991, Wilson and Sussman, 1993). In these earlier studies a single reaction mechanism is selected and used to simulate the M=6.46 shockinduced combustion case captured by Lehr (1972). A main focus of these earlier computations has been to quantify the effectiveness of particular numerical treatments. Wilson and MacCormack (1990) do compare the effect of altering the rate of the chainbranching reaction H + 02 OH + O and note the impact the alteration has on the induction time. However, none of the studies perform a comprehensive comparison over a wide range of reaction mechanisms which is one focus of the present investigation. A subset of those numerical techniques developed in these earlier studies is evaluated, particularly the source term treatments including the logarithmic interpolation proposed by Shang et al. (1995) and the scaling factor scheme proposed by Sussman (1993). Also evaluated is the explicitimplicit splitting scheme used to integrate the governing equations. Given the variety of reaction mechanisms that has been proposed for the hydrogenoxygen reaction process, the models used here have been selected based on their use in earlier computations of the shockinduced combustion scenario. These models are used to derive key parameters useful in quantifying reaction mechanisms and shown to be directly related to the computed flow fields. This establishes comparison parameters useful in future studies into the aspects of computing highspeed reacting flows and the coupling of the numerical and physical aspects. To assess the impact the numerical and physical aspects have on the accuracy of such computations, a starting point must be selected. The steady combustion scenario at M=6.46, captured by Lehr (1972) is chosen as the case to study. Based on the success of earlier computations, neither wall nor turbulence seems to have order one effects in the overall flow field and a logical starting point is to perform the computations with an inviscid model using piecewise representation of the source term. This is done in part to compare the current computations with those of earlier investigations. Based on the review of previous computational studies of Lehr's experiments (1972) it has been decided to first consider the 7step mechanism on a moderate grid and with no special treatment of the source term. The present computational grid is composed of 80x 120 cells with 120 being in the direction normal to the body. The current grid is shown in Figure 5.1. To simplify the issue related to the grid distribution, uniform spacing is used in the normal direction and the tangential spacing along the body is approximately uniform. The 7step hydrogenoxygen mechanism used in the baseline computation is given in Table 2.2 and is from the work of Drummond et al. (1986). The fully implicit solution Figure 5.1. Baseline 80 x 120 grid used for the M=6.46 case. scheme is used for the computations presented in this section. Also, for all shockinduced combustion computations, a thirdorder representation for the dependence of specific heats on temperature is used. The density contour from this baseline computation is shown in Figure 5.2 along with the reference points from the experimental shadowgraph presented earlier. The computation is carried out for only half the field but a mirror image is also presented to facilitate comparison to the experiments. The first noticeable feature of the comparison is the asymmetry in the experimental data, especially in the axial locations near the baseline of the projectile. It should be reiterated that there is also some optical distortion in the shadowgraph and the reference points have been altered in no way. They have simply been placed such that the experimental shock location on the stagnation line corresponds to * Experimental Shock Front o Experimental EnergyRelease Front Figure 5.2. Solution on the 80 x 120 grid for M=6.46 using the 7step mechanism and the piecewise constant representation of the reaction source term. that of the computations. Also, the data has been scaled 1:1 in length and height until the diameter of the experimental projectile matches that of the computational geometry. As is evident, this results in the reference points not extending as far to the base as in the experimental image of Figure 2.2. The data from the experimental shadowgraph can be altered to better match the length of the projectile in the computations but this requires additional assumptions on the optical distortion which are not made here. The comparison between the computed shock and reaction fronts are quite good in the nose region but in the region further downstream the reaction front is somewhat closer to the shock than the experimental data infers. This is one reason, unlike many of the previous studies, the current study models the complete length of the projectile to provide a more complete comparison. This thinning trend is also seen in the computations of Shang et al. (1995) in which the same reaction mechanism is used on a comparable grid. It is noted that all computations presented here have been converged to the point where the average density residual is reduced six orders of magnitude. Before studying other modeling aspects, the dependence of the solution on the grid resolution used in the computation must be measured. The 80x120 grid has been chosen based on the results of Shang et al. (1995) and the results from this grid as well as three alternate grids are shown in Figure 5.3. For comparison purposes, only the experimental data from the upper half of the shadowgraph in Figure 2.2 is used. To make a more precise comparison, in Figure 5.4 the variation in density along a line normal to the cylindrical body is presented. The x location of the lateral variation is one radius downstream of the hemispherical nosecylindrical body junction. All lengths are nondimensionalized by the projectile diameter. A similar comparison along the stagnation line results in lines that are essentially indistinguishable. The solutions on the coarser grids of 40x90 and 40x120 do alter the shock and reaction front location slightly, but even then the relative location of the reaction front to the shock are similar. Therefore, the thinning at the larger x must be due to factors other than the grid. Although Figure 5.4 demonstrates that as the grid is refined, the solution seems to asymptotically converge, caution must be exercised when determining grid dependency for highspeed reacting flows. This is due to the role of the chemical reaction source term which is dependent on the distribution of species and temperature. In the context of solving the multidimensional problem, this source term is calculated in each computational cell, or control volume. The discussion on source term treatments has shown that error introduced through the representation of the source term in each cell is proportional to kAx/u, the cell Damk6hler number, where k denotes the reaction rate and Ax the grid spacing. Figure 3.3 
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 EABKOPZZY_1T1PDQ INGEST_TIME 20131010T03:29:35Z PACKAGE AA00017651_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 