An Eulerian one-dimensional turbulence model: application to turbulent and multiphase reacting flows

This dissertation presents the development and validation of a variant of the One Dimensional Turbulence model (ODT) in an Eulerian reference frame. The ODT model solves unfiltered governing equations in one spatial dimension with a stochastic model for turbulence. The stand-alone ODT model implemented for this work resolves the full range of length and time scales associated with the flow, in 1D, with detailed chemistry, thermodynamics and transport in the gas phase. The model is first applied to a planar nonpremixed turbulent jet flame and results from the model prediction are compared with DNS data. Results indicate that the model accurately reproduces the DNS data set. Turbulence-chemistry interactions, including trends for extinction and reignition, are captured by the model. Differences observed between model prediction and data are the result of early excess extinction observed in the model. The reasons for the early extinction are discussed within the model context. A parameter sensitivity is also done for the current model. Simulations are performed over a range of jet Reynolds numbers for reacting and nonreacting configurations. Results from the simulations are compared with DNS and experimental data for reacting and nonreacting cases, respectively. Based on the identified sensitivity an empirical correlation is proposed and conclusions are drawn about the parameter estimation. The model is also applied to a planar premixed turbulent jet flame and results from the ODT simulations are compared with DNS data. Two different Da cases are considered in the study and comparisons between the model and DNS data in physical space are shown. Results indicate that the model qualitatively reproduces the DNS data set. Mixing is well captured by the model and the quantitative differences observed between model and data for thermochemistry are due to the curvature effects in the data. The reasons for the differences observed are discussed within the model context.

AN EULERIAN ONE-DIMENSIONAL TURBULENCE MODEL: APPLICATION TO TURBULENT AND MULTIPHASE REACTING FLOWS by Naveen Kumar Punati A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2012 Copyright c Naveen Kumar Punati 2012 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of has been approved by the following supervisory committee members: , Chair Date Approved , Member Date Approved , Member Date Approved , Member Date Approved , Member Date Approved and by , Chair of the Department of and by Charles A. Wight, Dean of The Graduate School. Naveen Kumar Punati James C. Sutherland 01/04/2012 Philip J. Smith 01/04/2012 Alan R. Kerstein 01/04/2012 Sean T. Smith 01/04/2012 Rodney C. Schmidt 01/04/2012 JoAnn Lighty Chemical Engineering ABSTRACT This dissertation presents the development and validation of a variant of the One Di-mensional Turbulence model (ODT) in an Eulerian reference frame. The ODT model solves unfiltered governing equations in one spatial dimension with a stochastic model for turbu-lence. The stand-alone ODT model implemented for this work resolves the full range of length and time scales associated with the flow, in 1D, with detailed chemistry, thermody-namics and transport in the gas phase. The model is first applied to a planar nonpremixed turbulent jet flame and results from the model prediction are compared with DNS data. Results indicate that the model accurately reproduces the DNS data set. Turbulence-chemistry interactions, including trends for extinction and reignition, are captured by the model. Differences observed between model prediction and data are the result of early excess extinction observed in the model. The reasons for the early extinction are discussed within the model context. A parameter sensitivity is also done for the current model. Simulations are performed over a range of jet Reynolds numbers for reacting and nonreacting configurations. Results from the simulations are compared with DNS and experimental data for reacting and nonreacting cases, respectively. Based on the identified sensitivity an empirical correlation is proposed and conclusions are drawn about the parameter estimation. The model is also applied to a planar premixed turbulent jet flame and results from the ODT simulations are compared with DNS data. Two different Da cases are considered in the study and comparisons between the model and DNS data in physical space are shown. Results indicate that the model qualitatively reproduces the DNS data set. Mixing is well captured by the model and the quantitative differences observed between model and data for thermochemistry are due to the curvature effects in the data. The reasons for the differences observed are discussed within the model context. The model is then extended to simulate a coal gasification process. A Lagrangian track-ing model is implemented for the particles, which are two-way coupled with the gas phase in the mass, momentum, and energy balance equations. A novel modeling technique is im-plemented for the particle-eddy interaction. For the coal particles, models are implemented for moisture vaporization, devolatilization of the raw coal and oxidation of the residual char. For this work, we consider the Chemical Percolation Devolatilization (CPD) model, which provides production rates of various gas-phase species during the devolatilization process. First a nonreacting particle laden jet simulation is performed and the results are compared with available experimental data. Results indicate that model qualitatively captures the par-ticle size influence on the dispersion behavior. For the coal gasification, simulation results are presented in the near field region of the jet. The model indicates that particle size has a significant influence on the initial heat up, vaporization and devolatilization processes. iv For my parents Subhadra and Basavaiah, and my brother Praveen, who have supported me in every possible way CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Gas Phase Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Coal Combustion and Gasification . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 One-Dimensional Turbulence Model . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.6 Outline of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 MODEL FORMULATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Governing Equations for ODT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Eddy Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3 NONPREMIXED TURBULENT JET FLAME . . . . . . . . . . . . . . . . . . . . 36 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3 Computational Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 PARAMETER SENSITIVITY ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . 48 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Turbulent Planar Jet Flame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3 Nonreacting Turbulent Planar Jet . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5 PREMIXED TURBULENT JET FLAME . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3 Computational Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6 TURBULENT PARTICLE LADEN JETS . . . . . . . . . . . . . . . . . . . . . . . . 92 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.3 Particle-Eddy Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.4 Computational Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7 COAL COMBUSTION AND GASIFICATION . . . . . . . . . . . . . . . . . . . . 115 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.3 Computational Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 7.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 8 CONCLUSIONS AND RECOMMENDATIONS . . . . . . . . . . . . . . . . . . . 129 8.1 Novel Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 8.2 Recommendations for Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . 130 A GOVERNING EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 B MODEL VERIFICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 C COAL MODELS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 vii ACKNOWLEDGMENTS I express my sincere and heartfelt gratitude to my mentor Dr. James C. Sutherland for his patience, guidance and motivation throughout my dissertation work. I cannot remember a time when James did not stop whatever he was doing in order to spend time answering questions and discussing issues with me. He has always been a constant source of motivation for me and is one of the very few people whom I look up to in my life. I have been very fortunate to be able to work with Alan Kerstein during my internship at Sandia National Laboratories, Livermore in summer 2009. Alan is the inventor of the ODT model and always been accessible to clarify my questions. The several discussions with him have been a rich learning experience. His passion for science and enthusiasm to help people truly inspires me. I also wish to thank my committee members, Professor Philip J. Smith, Dr. Sean T. Smith and Dr. Rodney Schmidt for taking the time to absorb this dense subject area and for providing direction and insight into the issues contained herein. I would like to thank the Department of Chemical Engineering at the University of Utah for supporting my graduate study. Thanks also to the administrative and support staff at the University of Utah for their help with issues through these years of graduate study. I would also like to acknowledge my colleagues, Michael Hradisky, Charles Reid, Julien Pedel, Tony Saad, Jeremey Thornock, Anchal Jatale and Babak Goshayeshi for all the invaluable discussions they had with me over different subjects. A special thanks to Yuxin Wu (Martin) for his invaluable suggestions for my dissertation work. I would like to thank my parents and brother for their constant encouragement and support. They are my pillars of strength and passion. I would also like to thank all my friends for their love and support throughout the years. I also acknowledge Department of Energy (DOE, Award Number FC26-08NT0005015) for providing financial grant for my dissertation work. CHAPTER 1 INTRODUCTION The conversion of chemical energy to sensible energy (heat) via a combustion process in a turbulent flow environment is necessary to meet ever-increasing energy demands. Combus-tion devices of practical interest include internal combustion engines, stationary and aircraft gas-turbine combustors, and industrial burners. The number of combustion systems used in the power generation and transportation industries are growing rapidly. This induces pol-lution and environmental problems to become critical factors in our societies. Combustion systems need to be operated such that the combustion reactions are brought to completion with a minimum of pollutants being formed. An accurate prediction of the essential physical and chemical properties of the combusting systems is important to achieve the two main objectives, optimization of combustion efficiency and the reduction of pollutants. In fact, turbulent combustion systems involve many phenomena and processes, such as turbulence, mixing, mass and heat transfer, radiation, and multiphase flow phenomena, which strongly interact. Their relative role depends on both the configuration and operating conditions. Turbulent combustion systems are often discussed in terms of the characteristic time scales required for mixing and reaction. If the mixing time scale (τm) is much higher the chemical time scale (τc) the assumption of fast chemistry (local chemical equilibrium) can be made. It is an assumption which introduces an important simplification, since it eliminates many parameters, those associated with chemical kinetics, from the analysis. This global comparison of time scales, however, may not be sufficient in turbulent flows where local diffusion time scales vary considerably. The fast chemistry assumption is then locally not valid and nonequilibrium effects must be taken into account. If the average diffusion time scale approaches the order of magnitude of the chemical time scale, local quenching will occur. A further reduction of diffusion time scales then leads to lift-off and even blow- 2 off of the entire turbulent flame. But already in globally stable flames the variation of diffusion time scales may interact selectively with the different chemical processes occurring in the system. Their chemical time scales may be quite different. For instance, the time required for combustion and generation of heat is much smaller than the time required for the formation of pollutants such as NOx and soot. Controlling of the time scales within the nonequilibrium range plays an important role in meeting the opposing requirements of fuel burnout, stability and low pollutant emission. "There is an axiom in physics which states that the simplest solution is usually the correct solution" [96]. In that vein, scientists often search for a simple, practical theory which will yield quantitative results for realistic problems in a relatively short time. For an excess of one hundred years physicists, mathematicians, and engineers have been searching for just that. However, a simple quantitative theory of turbulent combustion has not been identified [116]. No one can tell the future, but the likelihood of such a theory seems very distant in time. The scale-up from laboratory scale to industrial equipment is often a major problem, and it is often done by relying on experimental data and experience. The prediction of efficient operating conditions often based on empirical correlations, which tends to be unreliable. In the absence of a simple qualitative theory, the emphasis is on computing properties of turbulent combustion flows on computers. This work is intended to develop a numerical model, which will capture enough of the essential physics to predict quantitatively the turbulent combustion systems yet to be simple enough to be able to solve problems of practical interest. This chapter is organized as follows. First a brief introduction to computational fluid dynamics, gas phase combustion, turbulence models and coal combustion is given followed by the description of modeling technique adopted for the present work. This chapter concludes with an outline of the dissertation. 1.1 Computational Fluid Dynamics The philosophical study and development of the whole discipline of fluid dynamics is evolving with time. In the seventeenth century, the foundations of experimental fluid 3 dynamics were laid in France and England. The eighteenth and nineteenth century saw the gradual development of theoretical fluid dynamics. For most of the twentieth century the study and the practice of fluid dynamics involved the use of pure theory on the one hand and pure experiment on the other hand. The advent of the high speed digital computer combined with the development of accurate numerical algorithms for solving physical problems on these computers has revolutionized the way we study and practice fluid dynamics today. It has introduced a fundamentally important new third approach in fluid dynamics- the approach of computational fluid dynamics (CFD) [2]. In the present era CFD is an equal partner with pure theory and experiment in the analysis and solution of fluid dynamics problems and will continue to play this role indefinitely, for as long as out advanced human civilization exists. CFD is simply a new approach-but nothing more than that. It nicely complements the other approaches, pure theory and pure experiment, but will never replace either of these approaches. Fluid flow and related phenomena can be described by partial differential equations, which cannot be solved analytically except in special cases. In CFD, to obtain an approxi-mate solution numerically, a discretization method is used which approximates the differen-tial equations by a system of algebraic equations, which are then solved on a computer. The approximations are applied to small domains in space and/or time so the numerical solution provides results at discrete locations in space and time. When the governing equations are known accurately solution of any desired accuracy can be achieved in principle. However, for many phenomena (e.g., turbulence, combustion, and multiphase flow) the exact equations are either not available or numerical solution is not feasible [33]. This makes introduction of models a necessity. Even if we solve equations exactly, the solution might not be a correct representation of reality. In order to validate the models, we have to rely on experimental data. 1.2 Gas Phase Combustion Webster's dictionary defines combustion as "rapid oxidation generating heat, or both light and heat; also, slow oxidation accompanied by relatively little heat and no light." 4 Combustion is very complex and understanding the underlying chemical processes is essential in building more efficient systems. In many combustion processes chemical reaction controls the rate of combustion, and, in essentially all combustion processes, chemical rates determine pollutant formation and destruction. Ignition and flame extinction are intimately related to chemical processes. The study of the elementary reactions and their rates, chemical kinetics, is a specialized field of physical chemistry. Much progress has been made in understanding the combustion because the chemists have been able to define the detailed chemical pathways (for simple fuels) leading from reactants to products, and to measure or calculate their associated rates [1, 65]. With this knowledge, combustion scientists and engineers are able to construct computer models that simulate reacting systems. Combustion can be categorized into two different regimes based on mixedness of the reactants, i.e., premixed and nonpremixed [116]. If one looks at the complete range of the systems wherein turbulence and chemistry interact, one will find that many of the so called "mixing-sensitive" systems involve liquids or gas-phase reactions with modest density changes. For these systems, a key feature that distinguishes them from classical combusting systems is that the reaction rates are fast regardless of the temperature (e.g., acid-base chemistry). In contrast, much of the dynamical behavior of typical combusting systems is controlled by the fact that the reactants do not react at ambient temperatures. Combustion thus can be carried out in either premixed or nonpremixed modes, while mixing-sensitive reactions can only be carried out in nonpremixed mode [34]. 1.2.1 Premixed Combustion In a premixed flame, the fuel and oxidizer are mixed at the molecular level prior to the occurrence of any significant reaction. Fresh gases, fuel mixed with oxygen, and combustion products are separated by a thin reaction zone [119]. A strong temperature gradient exists between the fresh and burnt gases. In premixed flames the flame propagates towards the fresh gases. Because of the temperature gradient and the corresponding thermal fluxes, fresh gases are preheated and then start to burn. The most striking features of the premixed flames are counter-gradient diffusion and the large production of turbulent energy within 5 the flame [80]. Both these phenomena result from the large density difference between reactants and products and from the pressure field due to volume expansion. In applications, because of the explosion hazard, premixing is generally avoided. Nevertheless, there are several important applications of turbulent premixed combustion; the principal one is the (homogeneously charged) spark-ignition engine. Other examples are reheat systems in jet engines, industrial tunnel burners, and gaseous explosions in a turbulent atmosphere. 1.2.2 Nonpremixed Combustion In a nonpremixed flame, the reactants are initially separated, and reactions occur only at the interface between the fuel and oxidizer, where mixing and reaction both take place. Contrary to the premixed flame, in nonpremixed flames fuel and oxidizer are on both sides of a reaction zone where the heat is released. The burning rate is controlled by the molecular diffusion of the reactants toward the reaction zone (diffusion is the rate-controlling step); that is why nonpremixed flames are also referred as diffusion flames. The term diffusion applies strictly to the molecular diffusion of chemical species, i.e., fuel molecules diffuse towards the flame from one direction while oxidizer molecules diffuse toward the flame from the opposite direction. The turbulent convection mixes the fuel and air together on a macroscopic basis, whereas molecular mixing at the small scales then completes the mixing so that the chemical reactions can take place. Diffusion flames are mainly mixing controlled and the thickness of a diffusion flame is not a constant, but depends on the local flow properties. An example of a diffusion flame is a simple candle. 1.2.2.1 Extinction Efficient mixing is critical in nonpremixed combustion, because molecular mixing of reactants is necessary to allow chemical reaction. High molecular mixing rates characteristic of turbulent flows enhance reaction rates, and thereby improve combustion efficiency. How-ever, the interaction of finite-rate chemistry with excessive mixing rates (τm < τc) can lead to local extinction. Following this, local regions of fuel and oxidizer can mix and coexist without significant reaction, and may later reignite. 6 Extinction may lead to increased harmful emissions, and if pervasive, to flame desta-bilization or blowout. Understanding the intimate coupling between turbulence, molecular mixing, and finite-rate reaction is paramount to predicting the behavior of nonpremixed combustion processes. 1.3 Turbulence The governing equations describing the fluid flow phenomena can be found in many of the sources [9, 81]. The complex behavior of the Navier-Stokes equations has two general categories, one in which the viscous forces are extremely large, called laminar flow, and one in which the inertial forces are extremely large, called turbulence. The physics Nobel prize laureate Richard P. Feynman referred to turbulence as one of the last unsolved problems in physics. Many scientists have devoted their lives to studying turbulence and as such the volume of work on the subject is quite extensive. Turbulent flows are highly unsteady, three-dimensional and contain a great deal of vorticity. Stretching of vortices is one of the principal mechanisms by which the intensity of turbulence is increased [81]. Turbulence increases the rate at which conserved quantities are stirred. Stirring is a process in which parcels of fluid with different concentrations of at least one of the conserved properties are brought into contact. The actual mixing is accomplished by diffusion. Nonetheless, this behavior is often called diffusive. Turbulence brings fluids of differing momentum content into contact. The reduction of the velocity gradients due to the action of viscosity reduces the kinetic energy of the flow; in other words mixing is a dissipative process. The lost energy is irreversibly converted into internal energy of the fluid. It has been shown that turbulent flows contain coherent structures: repeatable and essentially deterministic events that are responsible for a large part of mixing [81]. However, the random part of turbulent flows causes these events to differ from each other in size strength, and time interval between occurrences, making study of them very difficult. Turbulent flows fluctuate on a broad range of length and time scales. The nonlinear terms in the Navier-Stokes equation and the pressure term make turbu-lent fluid flows difficult to solve even on the fastest computers [81]. With this in mind people 7 have made simplifications commensurate with what their needs and available tools (such as computer speed) were. Different people have different ways of categorizing ways of modeling turbulence. Based on the computational cost associated turbulence models can generally be classified into three groups: direct numerical simulation (DNS) models, Reynolds-averaged Navier Stokes (RANS) equation models, and large eddy simulation (LES) models [81]. 1.3.1 Turbulence Models As described in Section 1.1, the strategy in CFD is to approximate the continuous character of the flow and fluid properties with a discrete set of data. These data are dis-tributed across a computational domain and are mapped to spatial and temporal locations of the physical problem of interest. When enough spatial and temporal points are used to capture the smallest motions of the flow it is said that we are performing a DNS of the Navier-Stokes equations. DNS also utilizes high-order numerical methods to marginalize the impact of numerical error on the simulation results and also minimizes modeling error. Therefore DNS is a standard to which other turbulence models can be compared. In order to obtain the representation of the instantaneous velocity as a function of position (3-D) and time DNS must resolve all scales to the smallest (Kolomogorov) length scale [81]. Although the computational cost of such a calculation restricts DNS to small Reynolds numbers and simple geometries, considerable work has been done with respect to incorporating complex chemical kinetics and studying flame turbulence interaction [21, 22] using this method. LES utilizes a ‘filtered' velocity field to obtain the flow simulation. The LES strategy is to resolve scales far enough below the flow-dependent energy-containing scales so that the unresolved motions are within the inertial subrange, whose properties are presumed to be universal [94]. The fundamental questions about the conceptual foundations of LES, and about the methodologies and protocols used in its application are discussed by Pope in [82]. In LES the grid is not fine enough to capture all the energy containing motions, hence a substantial portion of the energy is in the Subgrid Scale (SGS). Models are required for the SGS and a LES simulation is a priori more dependent on the SGS modeling than a DNS. 8 RANS is the oldest and probably the most widely used of the methods for modeling industrial scale problems. Once the Navier-Stokes equations are Reynolds averaged there appear in the equations more terms than there are constitutive equations. Information about the turbulent fluctuations is lost in the averaging process, leading to the classic "closure problem" in turbulence. For detailed discussion on methods for solving the RANS closure problem please refer to [81]. RANS has been the CFD strategy for engineering applications for the last 30 years and to this day continues to be the most common way to solve problems with complex geometry. 1.4 Coal Combustion and Gasification Coal as an energy carrier plays an important role in the energy market and is a difficult fossil fuel to consume efficiently and cleanly. Compared with other fossil resources, coal has much greater reserve and involves lower costs, and so, is expected to remain an essential energy resource into the 21st century. One principal user of coal is the power plant, where pulverized coal combustion has become the generally accepted combustion system because of its excellent capacity to increase power production [7]. The combustion of pulverized coal is a complex process, involving coupled effects among heat and mass transfer, fluid mechanics, and chemical kinetics. Coal gasification offers a versatile and clean method for converting coal into gaseous fuel. In entrained flow gasifier, coal or coal slurry particles are usually injected into the furnace with pure oxygen at a high speed. Usually gasification process is carried out at high pressures and temperatures. The elevated pressure and high temperature in the gasifier guarantees a high carbon conversion in a short residence time. Under these conditions, the coal is broken apart into a gaseous mixture of CO and H2, which compose syngas fuel, the primary product of coal gasification, along with other products, such as CO2 and H2O. In addition to producing combustible gaseous fuel, coal gasifiers are also more efficient than traditional coal-fired boilers, both in thermal conversion of energy and in power cycle design. The gasification process in an entrained flow coal gasifier is very complex. A series of physical and chemical processes happen on the coal slurry particles, such as evaporation of 9 water, pyrolysis of coal, and heterogeneous coal char reactions. At the same time, there are strong coupling effects among turbulent fluctuation, chemical reactions, and heat transfer to particles. Especially, temperature and local velocities have a strong influence on these coupling effects, which makes the controlling mechanism and turbulent fluctuation effects change at different regions of an entrained flow gasifier. The coal combustion/gasification systems are turbulent and multiphase (solid-gas cou-pling) in nature. Both homogeneous and heterogeneous reactions occur in these systems. Modeling of such a complicated system needs accurate modeling of subprocesses. However nonlinear coupling occurring across a multitude of length and time scales in these systems poses a formidable challenge for accurate modeling. Even with the modern day computers, resolving the entire physics of the problem remains unfeasible. 1.5 One-Dimensional Turbulence Model Section 1.3.1 described the different turbulence models employed for the study of turbulent and multiphase combustion phenomena. In 1999, Kerstein developed the One-dimensional turbulence model [58]. As the name suggests in ODT the domain is restricted to 1D and the one dimensional line represents a line of sight through a three dimensional turbulent flow field. ODT simulates the evolution of fluid flow by completely resolving all the spatial and temporal scales along this line. In a loose sense, ODT is a one-dimensional surrogate for DNS. Being 1D, however, it does not suffer from the "curse of dimensionality" which makes DNS intractable for even modestly turbulent flow [69]. The distinctive feature of the model is the representation of turbulent advection by a postulated stochastic process rather than an evolution equation and the key attribute is computationally affordable res-olution of viscous scales in fully developed turbulence. However ODT is applicable only to flows that are homogeneous in at least one spatial coordinate. Many flows of fundamental interest and practical importance are of this type. ODT is an outgrowth of the linear-eddy model (LEM). In LEM, flow properties are specified empirically by assigning parameters governing the random event sequence. There is no provision for feedback of local flow properties to the random process governing subsequent 10 events. In contrast, ODT is formulated to capture this feedback with minimal empiricism. In this regard, ODT is both a turbulence model and a methodology for fully resolved simulation of mixing, chemical reaction, and related scalar processes in turbulence. The latter capability is a key feature distinguishing ODT from conventional turbulence models (LES and RANS) that require the incorporation of mixing submodels in order to treat scalar processes. The distinguishing features of ODT are its scope, simplicity, minimal empiricism, and capability to incorporate complex molecular processes (variable transport properties, chem-ical reactions, etc.) without introducing additional approximations [58]. Because ODT is a fully resolved simulation, various statistical quantities can be extracted that are not pro-vided by conventional closure methods. Being low dimensional the model is an inexpensive tool and it can be applied to problems of practical interest. The ODT model implemented for this work resolves full range of length and time scales associated the flow in 1D with detailed chemistry, thermodynamics and transport in the gas phase. More details of the ODT model used for the current study are given in the subsequent chapters. 1.6 Outline of the Dissertation The dissertation consists of formulation and validation of a new variant of the ODT model when applied to different set of problems. Chapter 2 presents a treatment to cast the various ODT formulations in a unified manner and a new variant in an Eulerian reference frame is described. Several derivations relevant to equations given in Chapter 2 are covered in Appendix A.1. This chapter establishes a mathematically sound basis for the various ODT formulations that will allow more clarity in comparing various approaches and will also allow a clear distinction between the equations being solved and the numerical method/algorithm used to solve the equations. Chapter 3 demonstrates the new ODT model's capability in reproducing the statistics for a nonpremixed reacting jet flame. The main focus of this chapter is to identify whether ODT model can capture significant finite rate chemistry effects like extinction and reignition. The results from the ODT simulation are compared with DNS data. 11 The main focus of the Chapter 4 is to identify the sensitivity of the model parameters. Simulations are performed over a range of jet Reynolds numbers for reacting and nonreacting configurations. Results from the simulations are compared with DNS and experimental data for reacting and nonreacting cases respectively. Based on the identified sensitivity an empirical correlation is proposed and conclusions are drawn about the parameter estimation. Chapter 5 demonstrates the model's ability to predict the important statistics for pre-mixed reacting jet flames. Results from the model are compared with DNS data. Mean profiles of velocities and temperature along with minor species are presented. Important statistics of premixed jet flames like flame surface density and surface area ratio are also compared with the data. In Chapter 6 the model is extended to simulate particle laden jets. For dispersed phase (solid particles), governing equations are derived in Lagrangian reference frame and two-way coupling, on momentum, between continuous and dispersed phases is implemented. A new particle-eddy interaction model is implemented to accommodate the eddy effects on dispersed phase. Turbulent particle laden jet simulations are performed and results are compared with available experimental data. Chapter 7 is an extension to Chapter 6, which mainly addresses the coal combustion and gasification process. Models describing the coal physics are implemented and the two-way coupling is extended for mass and energy. Qualitative assessment has been done for ODT coal gasification simulations. This dissertation concludes with a discussion on the findings from this study and some recommendations for future work in Chapter 8. CHAPTER 2 MODEL FORMULATION This chapter appears in much the same form as it is published in the technical report by Sutherland et al. [113]. Dr. Sutherland is the lead author of the report and I mainly contributed to the eddy events section. 2.1 Introduction The One-Dimensional Turbulence (ODT) model represents, conceptually, a line of sight through a turbulent flow field. First proposed by Kerstein [58], ODT is an outgrowth of the Linear Eddy Model [55-57, 71], but includes the solution of the local velocity field to determine the rate and location of eddy occurrence. Although ODT (and its predecessor LEM) has been implemented as a subgrid scale model in LES and RANS (see, e.g., [20, 36, 69, 70, 72, 88, 89]), much of its application has been as a stand-alone model. In stand-alone applications, ODT is applicable in situations where there is a direction of predominant large-scale gradients such as shear-driven flow (channels, jets), buoyancy-driven flow (plumes), etc. The one-dimensional domain is aligned perpendicular to the direction of primary gradients (e.g., across the shear layer), thereby resolving the primary driving force for turbulence. Of primary importance in ODT modeling is resolution of the streamwise (x-direction) velocity component (perpendicular to the direction of the ODT domain), as this velocity component captures the shear that results in the turbulent cascade. Indeed, early ODT formulations considered only the streamwise component of velocity. Later, the model was extended to include multiple components of velocity [59]. We refer to the ODT-aligned coordinate as the y-direction, and the streamwise coordinate as x throughout this chapter. This inherently assumes a Cartesian coordinate system, which is consistent with most ODT applications to date. However, ODT has been formulated in cylindrical coordinates as 13 well [64], where the ODT line is oriented in the radial direction, and the axial velocity component is the critical one that drives turbulence in the ODT model. ODT has been successfully applied as a stand-alone model for a variety of shear-dominated flows, both nonreacting [3, 35, 58, 59] and reacting [29, 50, 51, 64, 85]. The ODT model consists of two primary ingredients: • The governing equations written in terms of two independent variables: (t, y) for "temporal" ODT formulations and (x, y) for "spatial" ODT formulations. • Discrete "eddy events" that occur at various points in (t, y) or (x, y). In ODT, these eddy events are influenced by the local shear rate. Therefore, the majority of ODT formulations solve an equation to evolve the streamwise component of velocity. A notable exception is application of ODT to Rayleigh convection [124]. Stand-alone ODT models (the focus of the remainder of this chapter) have been formulated as temporally evolving, with (t, y) as independent variables, and spatially evolving, with (x, y) as independent variables. With each of these approaches, both Lagrangian and Eu-lerian variants can be used. Particularly in the case of variable-density flows, virtually all of the literature regarding ODT combines the numerical solution algorithm with the discus-sion of the governing equations so that it is not immediately clear what the actual governing equations being solved are. In some cases, the equations presented are not the equations being solved. In this chapter, we formulate the various stand-alone ODT approaches under a single umbrella and illustrate the differences between them. This is done without discussion of specific numerical algorithms, except in cases to illustrate nuances of implementations presented in the literature. We hope to establish a mathematically sound basis for the var-ious ODT formulations that will allow more clarity in comparing various approaches and will also allow a clear distinction between the equations being solved and the numerical method/algorithm used to solve the equations. The remainder of this chapter is organized as follows: Section 2.2 presents the governing equations solved for the ODT variants currently existing in the literature. The key modeling 14 concepts in ODT, the triplet map and kernel transformation, are then addressed in Section 2.3. Both of these sections are supplemented with material presented in the appendix. 2.2 Governing Equations for ODT In this section, we present several forms of the governing equations for use in ODT simulations. As shown in Appendix A.1, the governing equations for a single phase reacting system can be written in Lagrangian form as d dt ˆ V(t) ρψ dV = − ˆ V(t) Φψ · a dS + ˆ V(t) σψ dV, (2.1) ρ dψ dt = −∇ · Φψ + σψ, (2.2) where ψ is an intensive quantity, σψ is the net rate of production of ρψ, and Φψ is the non-convective flux of ρψ. Equations (2.1) and (2.2) are, respectively, the integral and dif-ferential forms of the Lagrangian evolution equations. In the Eulerian frame of reference, the corresponding equations are written in integral, strong differential, and weak differential forms as1 ˆ V(t) ∂ρψ ∂t dV + ˆ S(t) ρψv · a dS = − ˆ S(t) Φψ · a dS + ˆ V(t) σψ dV, (2.3) ∂ρψ ∂t + ∇ · ρψv = −∇ · Φψ + σψ, (2.4) ρ ∂ψ ∂t + ρv · ∇ψ = −∇ · Φψ + σψ, (2.5) where v is the mass-averaged velocity. Equation (2.2) is often most convenient for math-ematical manipulation, while (2.1) and (2.3) are more readily applied in a finite-volume numerical solution approach. Table 2.1 shows the forms of ψ, Φψ, and σψ for several com-mon forms of the governing equations. Appendix A.1 presents a derivation of the above equations and the terms in Table 2.1. 1Note that the weak differential form is obtained by applying chain rule to (2.4) and substituting the continuity equation (ψ = 1). 15 Table 2.1: Forms of terms in the governing equations (2.1)-(2.5). Here p is the pressure, τ = −μ Ä ∇v + (∇v)T ä + 2 3μ I∇·v is the stress tensor, v is the mass-averaged velocity, g is the gravitational vector, Yi is the mass fraction of species i, ji is the mass-diffusive flux of species i relative to a mass-averaged velocity, q = −λ∇T + ni =1 hiji is the heat flux, λ is the thermal conductivity, T is the temperature, and hi is the enthalpy of species i. Equation ψ Nonconvective Flux, Φψ Source Term, σψ Continuity 1 0 0 Momentum v pI + τ ρg Species Yi ji σi Total Internal Energy e0 pv − τ · v + q ρg · v Internal Energy e q −τ : ∇v − p∇ · v Enthalpy h q ∂p ∂t + v · ∇p − τ : ∇v ODT formulations can be broadly classified into two categories: • Temporal evolution, where (t, y) are chosen as the independent variables, and (through-out this chapter) y refers to the direction associated with the one-dimensional ODT domain. • Spatial evolution, where (x, y) are chosen as the independent variables and x refers to the streamwise direction. In each of these categories, there are both Eulerian and Lagrangian variants of ODT. Early ODT implementations were temporally evolving in a Lagrangian frame of reference. The spa-tially evolving ODT formulation was first introduced in the original description of ODT [58]. However, recent improvements to the model significantly broadened the range of flows and phenomena that the model can address [3, 64, 91]. Recently, Eulerian formulations for the temporal and spatial evolution have been demonstrated [84, 85]. The following sections consider each of these variants in turn. 16 2.2.1 Temporal ODT Evolution Equations We first consider temporal evolution. In this context, the ODT equations will describe evolution of various quantities on a line oriented in the y-direction and evolving in time. We consider two general forms of the governing equations: the Eulerian and Lagrangian forms. 2.2.1.1 Eulerian Temporal Form Retaining (t, y) as independent variables, (2.4) becomes ∂ρψ ∂t = −∂ρψv ∂y − ∂Φψ,y ∂y + σψ, (2.6) where Φψ,y = Φψ · y represents the component of Φψ in the y-direction, and v represents the local mass-averaged fluid velocity in the y-direction. Current approaches using the Eulerian form have solved the compressible form of these equations [84, 85]. Specifically, (2.6) (or its integral form equivalent) is solved as follows: • ψ = 1 is solved for ρ. • ψ = u is solved for the streamwise momentum (ρu) to provide the required information for the eddy selection (see Section 2.3). Note that the pressure gradient could be retained in this equation and imposed for pressure driven flow. • ψ = v is solved for the lateral momentum (ρv), which is mainly required for the continuity equation (ψ = 1). The pressure obtained from the equation of state is used to calculate the pressure gradient that appears in this equation. • ψ = e0 is solved for the total internal energy (ρe0). • ψ = Yi is solved for the species masses (ρYi). These equations are completed by an equation of state, p = p (ρ, T, Yi) . In a finite-volume context (as implemented in [84, 85]), the integral form of (2.6) is solved. While one could also solve the weak form of these equations (corresponding to equation (2.5)), there have been no implementations to date using the weak form. 17 2.2.1.2 Lagrangian Temporal Form In the Lagrangian frame of reference, choosing (t, y) as independent variables, (2.2) becomes ρ dψ dt = −∂Φψ,y ∂y + σψ. (2.7) Equation (2.7) is written in integral form as2 d dt ˆ ρψ dy = ˆ Å −∂Φψ,y ∂y + σψ ã dy. (2.8) Equations (2.7) and (2.8) are the forms most often used for temporally evolving ODT simu-lations. In the Lagrangian reference frame, the volume of a finite material element V(t), and its associated surface S(t), change with time according to the local mass-averaged velocity, v. To determine the locations of the cell centroids (and faces) ODEs may be solved for positions of cell centroids or faces by dy dt = v, (2.9) where v is the y-component of velocity. If we solve (2.7) for ψ = v then we have the lateral velocity component required for use in (2.9). Since v is the mass-averaged velocity, (2.9) describes the evolution of a surface defining a closed system for the mass, thereby enforcing continuity, d dt ˆ ρ dy = dm dt = 0. (2.10) Specifically, the limits of the integral in (2.8) are determined by (2.9), and (2.8) with ψ = 1 implies (2.10). Thus, by solving (2.9), we evolve the size of the control volume that, by definition, enforces continuity. Note that (2.10) need not be solved because it simply states that mass is constant.3 2See Section A.5 for details. 3Note that in the case of a multiphase system, where the continuity equation for one phase may have source terms due to interphase mass transfer, (2.9) is still the appropriate equation for enforcing continuity. 18 A full solution approach would solve: • Equation (2.9) for cell face positions to define the limits on the integral in (2.8). • Equation (2.10) need not be solved since its solution is simply that mass is constant. Density is obtained using this constant mass and the volume determined by the evo-lution of (2.9). • Equation (2.8) with ψ = u for the streamwise momentum to provide the required information for the eddy selection (see Section 2.3). Note that for pressure driven flow, the pressure gradient term can be specified accordingly. Otherwise, the streamwise pressure gradient is ignored. • Equation (2.8) with ψ = v for the lateral momentum. This is required for use in (2.9). • Equation (2.8) with ψ = e0, e, or h for energy conservation. • Equation (2.8) with ψ = Yi for species. • An equation of state p = p (ρ, T, Yi) . This is used in the lateral momentum equation (ψ = v). Early ODT formulations did not solve the y-component of velocity (v), and even among the ones that do (e.g., [59]), it is typically not used to supply the velocity for (2.9). Indeed, most ODT formulations to date use the y and z velocity components as repositories of kinetic energy rather than advective velocities. Thus, even if v is solved, rather than solving (2.9) to determine the limits for the integral in (2.8), (2.10) is used to describe the change in cell size. Specifically, (2.10) is discretized using a first-order time approximation to find (Δy)n+1 = ρn (Δy)n ρn+1 , (2.11) where the density is calculated from an equation of state (typically assuming constant pres-sure). This determines the new cell size, but does not specify position. In a time-split scheme, ρn+1 is evaluated from the equation of state, Y n+1 i , and Tn+1, where the pressure is typically assumed to be constant. Given the cell sizes at time n+1, the new cell positions However, the mass of the system will no longer be constant, and (2.10) (with the appropriate interphase exchange terms) would need to be evolved. 19 are determined by calculating = i (Δy)n i − i (Δy)n+1 i , adding /2 to the volume on each end of the domain, and then redistributing the control volumes with their new sizes over the domain length (which remains fixed). This approach has been employed in all variable-density temporal Lagrangian ODT simulations to date. Notably, it imposes a fixed domain size, whereas solution of (2.9) does not. To summarize, most current ODT temporal Lagrangian formulations solve (2.8) with ψ = u, ψ = e0 (or an equivalent energy variable), ψ = Yi, and (2.10) to obtain cell volumes that maintain continuity. However, as was shown in this section, an alternative would be to solve (2.7) for ψ = v and use this in (2.9) to obtain the positions of the Lagrangian cell centroids and faces. 2.2.1.3 Space-Time Mapping It is often useful to transform the time coordinate to an equivalent spatial coordinate. This can be done by solving an ODE for downstream position, and can be done in one of two ways: dx dt = u, (2.12) dx dt = u, (2.13) where u is the x (streamwise) component of the velocity. Equation (2.12) uses a suitably chosen average velocity (u) to determine the downstream position for the ODT domain whereas (2.13) uses the local velocity at each point on the ODT line and solves a position equation for each point. Figure 2.1 illustrates the difference between these approaches for a hypothetical constant (in time) u profile and u chosen in two different ways4: u = u∞ + ´ ρ (u − u∞)2 ´ dy ρ (u − u∞) dy , (2.14) u = ¨ u| u>uc ∂ , uc = α (max u − min u) (2.15) 4Note that these are only two of many reasonable choices for u. 20 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 y/δ x/δ τ/2 τ 2τ Fig. 2.1: Downstream position (x) as a function of lateral position (y) for various times. The solid line uses (2.13) (solid line), the dashed line uses (2.12) with (2.14), and the dotted line uses (2.15) with α = 0.05 (dotted line). Figure 2.1 clearly shows that the space-time mapping can be highly approximate, and must be used cautiously. 2.2.2 Spatially Evolving ODT Equations Because of the ambiguity in determining a downstream location (x) in the temporally evolving approach, it may be advantageous in some situations to formulate the governing equations so that (x, y) rather than (t, y) are the independent variables. Below we consider Eulerian and Lagrangian equation sets that use (x, y) as independent variables. 2.2.2.1 Eulerian Spatial Form The spatially evolving governing equations retain only (x, y) as independent variables in (2.4) to obtain ∂ρψu ∂x = −∂ρψv ∂y − ∂Φψ,x ∂x − ∂Φψ,y ∂y + σψ. (2.16) This is an elliptic equation, and is not readily amenable for use with the stochastic eddy events in ODT. If we neglect the term ∂Φψ,x ∂x , then we have 21 ∂ρψu ∂x = −∂ρψv ∂y − ∂Φψ,y ∂y + σψ, (2.17) which is an incompletely parabolic (convection-diffusion) equation set that may be solved using the method of lines for the streamwise fluxes, ρψu. Note that the continuity equation (ψ = 1), ∂ρu ∂x + ∂ρv ∂y = 0, (2.18) conserves mass flux rather than mass itself. The full set of equations to be solved is: (2.18), (2.17) with ψ = { u v Yi e0 }, and an equation of state. Note that we can obtain u = ρuu ρu , (2.19) ρ = ρu u = (ρu)2 ρuu , (2.20) ψ = ρψu ρu , ψ = { 1 u }. (2.21) An alternative solution strategy is to solve the weak form of (2.17), ∂ψ ∂x = − 1 ρu ï ρv ∂ψ ∂y + ∂Φψ,y ∂y − σψ ò , (2.22) together with an alternate form of (2.18) ∂ρ ∂x = −1 u ï ρ ∂u ∂x + ∂ρv ∂y ò . (2.23) The term ∂u ∂x in (2.23) may obtained from (2.22) with ψ = u. The full set of equations to be solved is: (2.23), (2.22) with ψ = { u v Yi e0 }, and an equation of state. As discussed in Section 2.2.1.1, the ∂p ∂y term comes from the equation of state while ∂p ∂x is only nonzero in the case where a pressure driven flow is considered, in which case a fixed value of ∂p ∂x is assigned. 22 2.2.2.2 Lagrangian Spatial Form Given that the independent variables for the spatial evolution equations are (x, y), we can write d dx = ∂ ∂x + dy dx ∂ ∂y = ∂ ∂x + v u ∂ ∂y . (2.24) Using (2.24), we can rewrite (2.22) in Lagrangian form as dψ dx = − 1 ρu ï ∂Φψ,y ∂y − σψ ò . (2.25) This applies to all ψ except ψ = 1 (continuity), since (2.25) is in weak form. The Lagrangian form of the continuity equation can be obtained by substituting (2.24) into (2.18) to obtain dρu dx = v u ∂ρu ∂y + ∂ρv ∂y . (2.26) In (2.25), dψ dx is interpreted as the local rate of change in ψ as it moves with velocity v. Equation (2.25) can also be written in integral form as d dx ˆ ρψu dy = ˆ Å −∂Φψ,y ∂y + σψ ã dy. (2.27) This is most easily shown by applying Leibniz' rule to (2.27) to find d dx ˆ y2(x) y1(x) ρψu dy = ρ2ψ2u2 dy2 dx − ρ1ψ1u1 dy1 dx + ˆ y2(x) y1(x) ∂ρψu ∂x dy = ρ2ψ2v2 − ρ1ψ1v1 + ˆ y2(x) y1(x) ∂ρψu ∂x dy = ˆ y2(x) y1(x) Å ∂ρψu ∂x + ∂ρψv ∂y ã dy, (2.28) where subscripts 1 and 2 indicate that the quantities are evaluated at y1 and y2, respectively, and we have used (2.29). Equation (2.28) shows that (2.17) and (2.27) are equivalent. By 23 virtue of the derivation of (2.25) from (2.17), we conclude that (2.25) and (2.27) are also equivalent. When (2.25) is solved, an equation for y is also required to determine the position of the Lagrangian system dy dx = v u , (2.29) where u and v are the local fluid velocities in the x and y directions, respectively. If solving the integral form of the Lagrangian evolution equations, the position is required to determine the limits on the integral in (2.27) for each discrete volume element. If solving the differential form of the equations (via, e.g., a finite difference method) then the position is required to evaluate the fluxes and their divergences. In both cases, the role of the velocity is to maintain the proper definition of the Lagrangian control volume as discussed in Section 2.2.1.2. Together with an equation of state, (2.25) and (2.29) form a complete set of equations. When solving these equations, primitive variables are obtained using (2.19)-(2.21). To date, ODT implementations using the spatial form of the governing equations have solved (2.25) - see, e.g., [91]. All of these formulations present equations (2.17) (Eulerian forms) as the governing equations to be solved, but the form of the governing equations actually solved in these formulations is (2.25) (Lagrangian forms).5 As discussed in Section 2.2.1.2, the v component of velocity was not solved in the early ODT formulations. Rather than solving (2.29), these formulations (and the ones cited above) obtain Lagrangian position (y) via (2.27) with ψ = 1, d dx ˆ ρu dy = 0, (2.30) so that a first-order time discretization results in (Δy)n+1 = (ρuΔy)n (ρu)n+1 , (2.31) 5References [91] only present the discrete form of the equations they are actually solving, but they are the discrete form of (2.25). 24 where n refers to the solution at streamwise position xn while n+1 refers to the solution at position xn+1. From the updated volume sizes, the local positions are obtained as described in Section 2.2.1.2. Alternatively (and equivalently), an equation for v could be solved (including the pres-sure term as shown in Table 2.1) and (2.29) could be solved for cell and face positions. However, as with the analogous approach in Section 2.2.1.2, this has not yet been demon-strated in ODT. Independent of which approach is taken to obtain the position evolution, the evolution streamwise mass flux, ρu, need not be solved since it remains constant (as is evident from (2.27) with ψ = 1). The exception is for multiphase flow where there may be a source term in the flux continuity equation, as discussed in Section 2.2.1.2 for the mass conservation analogue. 2.2.2.3 Time-Space Mapping Occasionally in a spatially evolving formulation we are interested in determining a "residence time," e.g., in order to advanced a chemical-kinetic mechanism [29,50]. In analogy to the discussion in Section 2.2.1.3 this can be obtained by solving one of dt dx = u −1, (2.32) dt dx = u −1. (2.33) Equation (2.32) accounts for the variation of residence time due to variation in u while (2.33) obtains a characteristic residence time for the domain assuming that it moves with some characteristic velocity u. As discussed in Section 2.2.1.3, the choice for u in (2.33) is somewhat arbitrary. 2.2.3 Summary This section has presented four general approaches for ODT formulations. These can be categorized as temporally developing and spatially developing equations, with Lagrangian and Eulerian variants of each. When solving the Eulerian equations (see Sections 2.2.1.1 25 and 2.2.2.1), the y-component of velocity advects fluid and serves to enforce continuity. On the other hand, when solving the Lagrangian equations (see Sections 2.2.1.2 and 2.2.2.2) continuity reduces to a statement that mass (temporal form) or mass flux (spatial form) remains constant. However, in the Lagrangian form, the position must be evolved. This can be done one of two ways: 1. Use the y-component of velocity to determine the system position by solving (2.9) (temporal) or (2.29) (spatial). No boundary conditions are imposed on this ODE for position, but the boundary conditions on v velocity directly influence the evolution of this equation. 2. Use an operator-splitting approach along with a discrete form of the continuity equa-tion (2.11) (temporal) or (2.31) (spatial). This also requires imposition of boundary conditions directly, as discussed in Sections 2.2.1.2 and 2.2.2.2. The boundary conditions mentioned in these two options are important as they directly affect entrainment, large-scale flapping, etc. In addition, initial conditions may be particularly important in the case of spatially developing flows because of the approximation discussed in Section 2.2.2 that eliminated the elliptical nature of the problem [64]. 2.3 Eddy Events Stand-alone modeling of turbulent flows using ODT requires a dominant direction of the flow (which we refer to as the y-direction) to be identified a priori. To mimic the three-dimensional nature of turbulence in one spatial dimension, a stochastic process is adopted whereby motions that accelerate mixing are modeled through a series of stochastic rearrangement events. These events may be interpreted as the model analogue of individual turbulent eddies which are referred to as "eddy events" or simply "eddies". Each eddy event modifies fields by applying an instantaneous transformation over some spatial interval (y0, y0 + ), where y0 represents the eddy starting location and is the eddy length. A complete definition of the model for an eddy event requires specification of: 26 1. A procedure for selecting the candidate eddy starting location (y0), length ( ), and the eddy rate distribution (which is a function of y0 and ). The selection procedure of y0 and is described in Section 2.3.1. 2. The transformation (mapping), which is the effect of an eddy on the solution variables. The following sections address these elements of the eddy model. 2.3.1 Eddy Starting Location and Length The selection of eddy event starting location and length are described in this section. At each integration step, eddy length ( ) and location (y0) are selected from randomly generated numbers and flow properties. From scaling analysis, the eddy length ( ) can be defined as = −2Lp ln 2Lpr ceddy + ln Ä−2Lp min ä (2.34) where r is a random number, Re is Reynolds number of the flow, Lp = exp ln(L)+ln(η) 2 is the most probable eddy length, L is the integral length scale, η = L Re0.75 is the Kolmogorov length scale, min = 6η is the minimum eddy length, max = L is the maximum eddy length, and ceddy is some arbitrary constant defined as follows, ceddy = 2Lp exp Ä−2Lp max ä − exp Ä−2Lp min ä. While computing the probability of eddy occurance, f(y0) and g( ) need to be specified. The commonly used probability density functions of y0 and are, f(y0) = ceddy 2 exp Å−2Lp ã (2.35) g( ) = 1 max − min . (2.36) 27 2.3.2 Transformations In ODT, each eddy is an instantaneous event and has no opportunity to interact directly with other eddies. Rather, the interaction is indirect, mediated by the flow evolution. An eddy event is represented by an instantaneous rearrangement in the form of a "triplet map." For a selected eddy event the triplet map of a function ψ(y) is ψ (f (y; y0, )) , with f(y; y0, ) given as f(y; y0, ) ≡ y0 + ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎩ 3 (y − y0) y0 ≤ y ≤ y0 + 1 3 2 − 3 (y − y0) y0 + 1 3 ≤ y ≤ y0 + 2 3 3 (y − y0) − 2 y0 + 2 3 ≤ y ≤ y0 + y − y0 otherwise . (2.37) The triplet map defined by (2.37) forms the heart of any ODT modeling approach, rep-resenting the effects of a three-dimensional eddy with a one-dimensional rearrangement. Triplet maps are qualitatively similar to turbulence in that they have the effect of increasing gradients by redistributing the fluid elements along the 1-D domain. The functional form chosen for the triplet mapping function is the simplest of a class of mappings that satisfy the physical requirements of measure preservation (the nonlocal analog of vanishing velocity divergence), continuity (no introduction of discontinuities by the mapping operation) and scale locality (at most order-unity changes in property gradients) over the eddy interval and also strengthen the local stretch rate just as turbulent fluctuations do [58]. The desired attribute of the triplet map is to provide a means of mimicking the increase in strain inten-sity, the decrease in strain length scale and the increase in mixing due to eddies in actual turbulent flow. This mapping rule assures that closest neighbors after the mapping event were no more than 3 cells apart before the mapping event. Hence the increased strain rate and shortening length scale is attained without undue introduction of discontinuities. While the triplet map itself is measure preserving, occasionally we wish to augment the transformation imposed by the eddy event to ensure conservation of other properties. For example, application of the triplet map to ρψ results in conservation of momentum, energy, 28 and mass. However, kinetic energy is not necessarily conserved. If an eddy occurs in the presence of a gravitational field, then there is an exchange of potential and kinetic energy that must be accounted for when the transformation is applied to the ρψ. To ensure conservation of kinetic energy when applied to the momentum fields, the triplet map can be augmented by a "kernel transformation," ciK(y), which ensures conservation of kinetic energy. Applying such a kernel transformation to velocity components rather than momentum components can lead to a violation of momentum conservation, so that a second kernel, biJ(y), must be added to repair momentum conservation in the situation where transformations are applied to ψ rather than ρψ6. In this context, we can write the effect of an eddy on a velocity/momentum field as ψi(y) = ψi (f(y; y0, )) + ciK(y) + biJ(y), (2.38) where K(y) = y − f(y; y0, ), (2.39) J(y) = |K(y)| . (2.40) The K(y) kernel enforces conservation of kinetic energy while the J(y) kernel enforces con-servation of momentum [3]. If ρψ rather than ψ, is transformed, then J is not required (or bi = 0) since momentum will be conserved by construction when (2.38) is applied to ρu, ρv, ρw. Generally speaking, selection of a kernel transformation is influenced by two primary considerations: 1. What variables are being transformed? Typically this will be one of ψ, ρψ, since it it typically most convenient to transform the solution variables. However, one could impose a transformation on any set of variables in general. 2. What constraints are placed on the transformation? Most frequently we seek to impose constraints on the momentum and kinetic energy when solving a temporal form of the 6Note that in cases where there are momentum sources from, e.g., a dispersed phase, a kernel ensuring momentum conservation must be applied even if ρψ is transformed. 29 equations and on the fluxes of momentum and kinetic energy when solving the spatial form of the equations. However, these are modeling considerations, not fundamental requirements. Additional constraints could be added as necessary. The derivation of the transformations for various choices of transformed variables and con-straints is presented in Appendix A.2. The coefficients ci and bi can be represented generally as ci = 1 2S Ñ −Pi + sgn(Pi) P2 i + α j TijP2 j é , (2.41) bi = Hci, (2.42) with specific forms for S, Pi, and H presented in Appendix A.2 and summarized in Table 2.2. Table 2.2 summarizes the transformations for designed for ψ and ρψ with constraints on conservation of kinetic energy and momentum and the fluxes of kinetic energy and momen-tum. In the original ODT model formulation, a single velocity component was considered along with a set of scalars [58]. The application of the model to buoyant stratified flows, where conversion between kinetic energy and potential energy was the key concern, moti-vated the development of a kernel transformation to enforce the kinetic energy conserva-tion [123]. Subsequently, a "vector" ODT formulation was considered, where an eddy event incorporated energy transfer between velocity components [59]. The rules governing the partitioning of kinetic energy among velocity components, which have been referred to as the "pressure-scrambling model," also incorporate an element of three-dimensionality into the 1D model. This model is derived in Appendix A.2 for various ODT formulations, and the key results are summarized in Table 2.2. 2.3.3 Eddy Selection The procedure to select an eddy event is described here in the context of temporally evolving flows. A similar analysis with appropriate scaling can be used for spatial flows. Unlike LEM, where frequency and eddy-size distribution of the events are based on a prede- 30 Table 2.2: Transformations for various quantities given the constraint of conservation of (a) momentum and kinetic energy, and (b) fluxes of momentum and kinetic energy. Derivations of the results shown here are presented in Appendix A.2. See equations (2.41) and (2.42) for how S, Tij, Pi, and H (whose defining equations are enumerated in the table) are used in obtaining the kernel coefficients bi and ci for use in (2.38). Kernel Coefficients for Conserving Quantities Transformed Quantity Transformations Momentum and KE Flux of Momentum & KE Tij H Pi S Tij H Pi S ψ ui(y) → ui[f(y; y0, )] + ciK(y) + biJ(y), ψi (y) → ψi [f(y; y0, )], ψ = ui (A.46) (A.51) (A.52) (A.53) (A.46) (A.71) (A.70) (A.72) ρψ (ρui) [y] → (ρui) [f(y; y0, )] + ciK(y) (ρψ)i [y] → (ρψ)i [f(y; y0, )], ψ = ui (A.46) N/A (A.63) (A.62) (A.46) N/A (A.85) (A.86) 31 fined kinetic energy spectrum [55], the eddy events are influenced by the local flow field in ODT. Similar to the dimensional relationship applied to turbulent eddies, for events defined in ODT, a relationship can be formulated between an eddy's size, time scale, and kinetic energy. Since ODT resolves one or more components of the velocity/momentum vector, the "turnover" time (τe) for an eddy can be calculated from the local kinetic energy and the length of a candidate eddy. 7 From τe, the eddy rate distribution (λ) that governs the eddy events is calculated from λ = C 2τe , (2.43) where C is a model constant often referred to as the "eddy rate constant." The models used to identify the turnover time (τe) are summarized in Table 2.3 for different ODT model variants. The quantity /τe is interpreted as an eddy velocity and ρ 3/τ2 e is interpreted as a measure of the kinetic energy of eddy motion. Based on the streamwise velocity (x-component), the kinetic energy will be computed and equated to eddy energy to formulate an expression for eddy velocity. 8 The model constant Z that appears Table 2.3 is a "viscous cutoff" parameter that provides a lower limit on the eddy size roughly analogous to the Kolmogorov scale. In principle, this is not necessary (and could be set to zero) since eddies smaller than the Kolmogorov and Batchelor scales will have a negligible effect on the physical evolution of the system.9 Because ρψ and ψ evolve continuously in time between eddy events, λ also evolves continuously in time. The unsteadiness of the eddy rate distribution is both a fundamental 7Note that τe can be interpreted as an eddy turnover time or the time between eddies (inverse of the eddy frequency). These two quantities are closely related, but there are situations where a clear distinction is important, such as in particle-laden flows where particle-eddy interaction is important. In these cases, τ −1 e is interpreted as an eddy frequency governing eddy sampling and the eddy turnover time is calculated using an adjustable constant of proportionality [122]. 8More recent formulations that employ the "vector" formulation and solve several components of velocity use the kinetic energy from all velocity components in determining the eddy velocity and time scale [59]. 9See [58] for exceptions. 32 Table 2.3: Models for the eddy time scale, τe, in terms of the eddy velocity, /τe. More details are found in Appendix A.2. ODT variant Eddy time scale, τe Model Pa-rameters Relevant Equations ψ, conserving KE and momentum Ä τe ä2 ∼ 4 2 27ρKK 8 27 Ä Qy + α j TyjQj ä −Z μ2e ρe α, C, Z See Section A.1.1 and equations (A.46), (A.50), (A.55), (A.90), (A.91). ρψ, conserving KE and momentum Ä τe ä2 ∼ ï Q2+α j TyjQj ρe −Z νe 2 ò α, C, Z See Section A.1.2 and equations (A.46), (A.64), (A.91), (A.92). ψ, conserving KE and momentum fluxes Ä τe ä2 ∼ ´ 2 y0+ y0 ρu dy î Q2 + α j TyjQj −Z ρu,KK νe 2 ó α, C, Z See Section A.2.1 and equations (A.46), (A.69), (A.74), (A.92). ρψ, conserving KE and momentum fluxes Ä τe ä2 ∼ ñ ´ 1 y0+ y0 ρu dy Ä Q2 + α j TyjQj ä −Z νe 2 ô α, C, Z See Section A.2.2 and equations (A.46), (A.88), (A.92). 33 property of the model and a key consideration in its numerical implementation. λ is used to compute the probability of the eddy occuring, pe = λΔt f(y0) g( ) , (2.44) where f(y0) and g( ) are the probability density functions for y0 and , respectively. The functional forms for f(y0) and g( ) can influence the computational cost of the simulation, but do not affect the results [58]. The commonly used functional forms of f(y0) and g( ) are given in Section 2.3.1. The probability (pe) is compared with a randomly selected number on the interval [0, 1]. If the random number is less than pe then the eddy will be implemented. For spatially evolving flows (see Section 2.2.2) Δt is replaced by Δx/ˆu, where ˆu = 1 ˆ y0+ y0 u dy (2.45) is the average velocity defined over the eddy interval. This results in a definition of the probability of an eddy occurring of pe = λΔx f(y0) g( )ˆu (2.46) for spatially evolving flows. 2.3.3.1 Large Eddy Suppression While the viscous cutoff parameter Z suppresses the least energetic eddies, we require a mechanism to prevent the occurrence of unphysically large eddies that result in unreal-istic behaviour. We briefly outline the common methods for large eddy suppression in the following subsections. 2.3.3.1.1 Eddy time scale method. Even though eddies are implemented as instanta-neous events, the turnover time associated with each eddy event can be calculated from the 34 scaling analysis as summarized in Table 2.3. The eddy turnover time can be compared with the simulation elapsed time (t), and eddy events are allowed only when t ≥ βτe, where β is a model parameter. This large eddy suppression mechanism is used in most of the temporal formulations [29, 50, 51, 85]. However, this approach can also be used in a spatially evolving simulation by using x ≥ β as the criteria for eliminating large eddies. 2.3.3.1.2 Median method. In this method, the eddy event rate (λ) for a given eddy event is evaluated two different ways and the smaller of the two results is used in evaluating the probability. One evaluation is by the expressions formulated in Section 2.3.3. The other evaluation replaces each velocity profile ui(y) by a profile that is linear in y, and evaluates of λ based on these linear proles. The slope of each profile is taken to be the median value | dui dy | within the eddy range [y0; y0 + ]. The procedure assigns a zero rate to any event for which each velocity profile is flat (zero slope) in more than half of the eddy range, thereby suppressing large eddies [59]. 2.3.3.1.3 Scale reduction method. The scale reduction method is the most common method for suppressing large eddies in spatially developing flows [3,64,91], although it could be applied to temporally evolving flows as well. It involves auxiliary eddy-rate computations for each of three equal subintervals of the eddy interval [y0, y0 + ]. For the selected eddy event, λ is evaluated as if the eddy interval were î y0 + (j − 1) 3, y0 + j 3 ó , for j = 1, 2, and 3, respectively. If any of these three candidate eddies are disallowed due to dominance of the viscous penalty, as described in Section 2.3.3, then the eddy is discarded. Otherwise it is unchanged from its value computed for the complete eddy interval [y0, y0 + ]. 2.4 Conclusions We intend that this chapter will serve as a reference for those interested in ODT as a modeling approach by providing a survey of the various ODT formulations along with a sound mathematical basis for the equations being solved. Most ODT formulations (particularly for variable density flows) have not clearly dis-tinguished the governing equations being solved from the numerical method employed to 35 solve them. The equations are often written in fully discrete form. This chapter attempts to clarify the equations being solved by the various ODT formulations and, in so doing, raise alternative solution techniques. Specifically, this chapter has formulated the governing equations for use in ODT simulations in several forms: • Temporally evolving Eulerian, • Temporally evolving Lagrangian, • Spatially evolving Eulerian, • Spatially evolving Lagrangian. In addition, the models for "eddy events" in ODT, including the transformations applied to the solution variables (with appropriate kernel transformations) and the eddy selection criteria, were discussed and compared for the various ODT formulations. Both the governing equations and the variable tranformations associated with the eddy events are presented in a general manner assuming variable density and a multicomponent reacting system. Simplifications can be made in the event where density or composition is constant. In such cases, the discussion here simplifies to many of the early forms for ODT presented in the literature. CHAPTER 3 NONPREMIXED TURBULENT JET FLAME This chapter appears in much the same form as it is published in the article by Punati et al. [86]. 3.1 Introduction Predictive methods based on fundamental principles to model turbulence-chemistry interactions are important in turbulent reacting flow simulations to improve combustion efficiency and to reduce emissions. The existence of a wide range of length and time scales in high Reynolds number flows makes Direct Numerical Simulations (DNS) computationally intractable [10]. To reduce the computational cost one generally averages or filters the governing equations to remove fine scales as in the Reynolds-Averaged Navier-Stokes (RANS) and Large eddy simulation (LES) approaches. These averaged equations are coupled with turbulent combustion models to address the nonlinear nature of chemical reactions occurring at molecular mixing scales (fine scales). Turbulent combustion models can be broadly categorized into moment methods and probability density function (PDF) approaches [77]. In moment methods, molecular trans-port is explicitly represented and a reduced parameter space approach is adopted for the solution of reacting scalars and their associated source terms. For PDF approaches, chemical source terms appear in closed form whereas mixing is implemented stochastically using a mixing model. The linear-eddy model, developed by Kerstein [55-57], is one such stochas-tic mixing model which has been used as an alternative strategy for closure in turbulent combustion [71, 72, 95] . The main objective of the present study is to perform stand-alone Eulerian ODT sim-ulations for a nonpremixed temporally developing planar syngas jet flame and to compare 37 the model prediction with DNS data [46, 48]. This work is the first time that a stand-alone ODT model has been compared directly with 3D DNS data for a reacting flow, and also demonstrates the richness of the data the model can produce. It also represents one of the first attempts to model the DNS data set, the only other approach to date being a combi-nation of LEM with LES [98]. In the present work, all simulation details, including mesh spacing, initial conditions, boundary conditions, and thermodynamic, chemical kinetic and transport models were matched with the DNS. The chapter is organized as follows. First we present the governing equations that are solved and then we evaluate the model's capability to reproduce finite-rate chemistry effects such as extinction and reignition. Flow entrainment effects are presented using axial statistics for velocity and mixture fraction. Conditional statistics of species and probability density functions of temperature and scalar dissipation are presented. Fig. 3.1: Schematic of the DNS configuration for the syngas jet flame (Case M) showing the logarithm of the scalar dissipation rate [47, 48]. 38 3.2 Model Formulation The transverse y-direction, which is the direction of the most significant gradients (Fig-ure 3.1), is considered here as the ODT domain. 3.2.1 Governing Equations The ODT model as formulated herein solves the following conservation equation set ∂ρ ∂t = −∂v ∂y , (3.1) ∂ρv ∂t = −∂ρvv ∂y − ∂τyy ∂y − ∂P ∂y , (3.2) ∂ρu ∂t = −∂ρvu ∂y − ∂τyx ∂y , (3.3) ∂ρe0 ∂t = −∂ρe0v ∂y − ∂pv ∂y − ∂τyyv ∂y − ∂q ∂y , (3.4) ∂ρYi ∂t = −∂ρYiv ∂y − ∂Ji ∂y + ωi, (3.5) where u and v refer to streamwise and lateral velocities, ρ is the density, p is the pressure, τ is the stress tensor, e0 is the total internal energy, q is the heat flux, Yi is the mass fraction of species i, Ji is the species mass diffusive flux, and ωi is the reaction rate. These equations are completed with the ideal gas equation of state and constitutive relationships for the diffusive fluxes: τyy = −4 3 μ ∂v ∂y , (3.6) τyx = −μ ∂u ∂y , (3.7) q = −λ ∂T ∂y + ns i=1 hiJi, (3.8) Ji = −ρYi Xi Dmix i ∂Xi ∂y . (3.9) where λ is the thermal conductivity, μ is the viscosity, T is the temperaure, hi is the enthalpy of species i, Dmix i is the mixture-averaged diffusivity for species i, and Xi is the mole fraction of species i. Code verification details are included in Appendix B.1. 39 3.2.2 Eddy Events As discussed in Section 2.3, turbulent motions that accelerate mixing are modeled through a series of stochastic rearrangement events. Continuously evolving gas phase is subjected to instantaneous rearrangement through triplet mapping (see Section 2.3.2 for additional details, specifically transformation is applied on conserved variables ρ, ρu, ρv, ρe0, ρYi). To suppress the large eddies eddy time scale method is implemented. 3.3 Computational Configuration DNS of three-dimensional (3D) temporal planar syngas jet flames with detailed chem-istry over a range of jet Reynolds numbers (Re) from 2510 to 9079 have been performed by Hawkes et al. [46, 48]. We consider a case with Re = 4478, which is addressed as Case M in the literature. The jet consists of a central fuel stream (50% CO, 10% H2 and 40% N2 by volume) surrounded by counter-flowing oxidizer streams comprised of 25% O2 and 75% N2. The stoichiometric mixture fraction is Zst = 0.42 and the steady extinction dis-sipation rate (based on a steady laminar flamelet solution using erfc distribution on χ) is χq = 2194 s−1 [48]. The fuel and oxidizer stream bulk velocities are U/2 and −U/2, re-spectively, with U = 194 m/s. The initial fuel stream thickness is H ≈ 0.96 mm and the characteristic jet time scale, computed using H/U, is tj ≈ 5 μs. Based on tj , a nondimen-sionalized time parameter is defined as τ = t/tj . The mixture fraction was computed from the local species compositions using Bilger's definition [8], and the scalar dissipation rate is calculated as χ = 2D Ä dZ dy ä2 with D = λ/(ρcp). Figure 3.1 illustrates the scalar dissipation rate field of the jet at τ = 40. The DNS data set exhibits significant finite-rate chemistry effects including extinction and reignition. Maximum extinction occurs at τ ≈ 20, while by τ ≈ 40 most of the flame has reignited. The ODT calculations consider a one-dimensional domain aligned with the y-direction in Figure 3.1. The initial conditions for all the variables transported in the ODT model are extracted directly from the DNS data. The detailed chemical mechanism considered in this study (consisting of 11 species and 21 reactions [48]), temperature and pressure-dependent 40 thermodynamic property evaluation, and the mixture averaged transport treatment are all consistent with the DNS simulations. The spatial and temporal resolution are likewise the same as used in the DNS simulation, with a spatial resolution of 15 μm and a time step of 2 ns. Simulations are run for 0.25 ms (50 tj) and results are analyzed over 400 ODT simulation realizations, which was enough to provide stationary statistics. ODT results at different times are compared with DNS statistics on xz planes in Figure 3.1. 3.4 Results and Discussion 3.4.1 Flow Entrainment The flow entrainment predicted by the ODT model is evaluated by comparing axial evolution of velocity and mixture fraction at different time intervals. The entrainment is sensitive to the choice of ODT parameters (particularly β and C), which have been tuned to match the spreading rate and decay of the velocity and mixture fraction. The parameter values used in this study are α = 0.5, C = 60, β = 1.0 and Z = 200. Figure 3.2 shows comparison of observed and model behavior for axial evolution of the mean streamwise velocity at τ = 6, 20 and 40. A similar comparison is shown for the mixture fraction in Figure 3.3. For both velocity and mixture fraction DNS mean values on the left half of the domain are a mirror image of the right half of the domain but for ODT, data at positive and negative y are not combined (i.e., spatial profiles are not symmetrized). For both velocity and mixture fraction, the decay and spreading rate are very well reproduced by the model, demonstrating the model's efficacy in capturing the flow entrainment. 41 −5 0 5 −100 −50 0 50 100 Y/H < u >[m/s] ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Fig. 3.2: Average streamwise velocity at τ = 6, 20 and 40. −5 0 5 −100 −50 0 50 100 Y/H < u >[m/s] ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Fig. 3.3: Average mixture fraction profiles at τ = 6, 20 and 40. 42 3.4.2 Conditional Statistics Figures 3.4 and 3.5 show the mean temperature and OH evolution, respectively, as a function of the mixture fraction at different time intervals. The steady laminar flamelet solution at the critical dissipation rate (χq = 2194 s−1) is also shown for reference. For the case simulated, mixing is initially rapid enough relative to reaction to cause local extinction, which is followed by reignition as mixing rates relax. The conditional mean for both tem-perature and species predicted by ODT is low compared to the DNS data atτ = 6, 20 over the entire mixture fraction range. The ODT model starts predicting local extinction earlier (at τ ≈ 6) than the DNS, as indicated by both temperature and OH species mean values dropping close to the extinction limit predicted by the laminar flamelet solution. As the simulation progresses (at τ = 20) the ODT model exhibits stronger extinction than the data as indicated by the low mean values. At τ = 40 both T|Z and YOH|Z are above the values predicted by the steady flamelet model at χq, indicating that reignition has occurred. Also note that the ODT values of both T|Z and YOH|Z are larger than the DNS. 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 1800 Z < T >[K] ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Flamelet Fig. 3.4: Conditional mean temperature, T|Z , at τ = 6, 20 and 40. The steady flamelet solution at χq is also shown for reference . 43 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 x 10−3 Z < YOH > ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Flamelet Fig. 3.5: Conditional mean temperature, YOH|Z , at τ = 6, 20 and 40. The steady flamelet solution at χq is also shown for reference. Figures 3.6 and 3.7 show evolution of the conditional probability density function (PDF) of T and log10 (χ/χq), respectively near Zst = 0.42 at three different time intervals (τ = 6, 20 and 40). In the early stages of the simulation (at τ = 6), the scalar dissipation PDF evolution shows a narrower distribution and is shifted toward higher values relative to the DNS data. These higher values of χ cause extinction in the early stages of the ODT simulations, resulting in a corresponding temperature PDF shift toward lower values with the most probable state near the steady flamelet extinction limit of T ≈ 1250 K. The higher χ predicted by ODT in the early stages of development is followed by a decrease in χ that is more rapid than exhibited by the DNS data. At τ = 20 the mixing rates are still high enough to cause extinction in the model, as indicated by the tails of the PDF in Figure 3.7. During the later stages of the simulation (τ = 40), mixing rates relax as indicated by the dissipation PDF shift toward lower values and the temperature PDF evolution starts shifting toward high values as reignition occurs. 44 600 800 1000 1200 1400 1600 0 0.05 0.1 0.15 0.2 T[K] at Z=Zst PDF ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Fig. 3.6: Conditional PDF of T|Zst at τ = 6, 20 and 40. −6 −4 −2 0 2 0 0.05 0.1 log10(χ/χ q)|Z=Zst PDF ODT−τ=6 DNS−τ=6 ODT−τ=20 DNS−τ=20 ODT−τ=40 DNS−τ=40 Fig. 3.7: Conditional PDF of log10 χ/χq| Zst at τ = 6, 20 and 40. 45 There are three different modes through which reignition can happen: autoignition, triple or edge flame propagation and turbulent flame folding [110]. The dominant reignition mechanism for the present case is turbulent flame folding [45, 48], where neighboring flame segments that are vigorously burning can provide a source for reignition. ODT can capture this mode of reignition because of triplet map action during an eddy event. The triplet map (2.37) instantly rearranges momentum and scalar fields enabling heat transfer from burning to nonburning regions. Since the domain is restricted to one dimension in ODT, the triple flame reignition mode (for which nonaligned gradients of mixture fraction and progress variable are needed [51]) cannot be addressed. Figure 3.8 shows comparison of χ|Z as a function of time for stoichiometric and fuel rich regions. The ODT model predicts higher χ|Z (exceeding χq), in the early stages, for both the regions compared to the DNS data, and as a result, early extinction occurs. In the fuel rich region χ|Z starts increasing as the simulation starts and exceeds χq as early as τ = 2 and starts decaying from τ = 3, whereas the corresponding times for the stoichiometric region are 4 and 6, respectively. The DNS data also exhibits regions in which χ exceeds χq consistent with the model results; however the maximum mean χ occurs later than in the model, resulting in earlier occurrence of extinction. In subsequent stages, the predicted χ|Z is lower than the DNS data. However, extinction continues until τ = 20 which is evident from Figure 3.6. In the later stages of the simulation, χ|Z predicted by ODT is lower in both regions compared to the DNS data, indicating a faster increase towards equilibrium resulting in higher mean temperature and OH species concentration. Early extinction observed in the model is the primary reason for the discrepancies observed between the ODT and DNS data. In the ODT model, the large eddies control the flow entrainment [29], which is well reproduced by the model (see Figures 3.2 and 3.3), whereas small eddies influence the small-scale mixing. Because of the high shear available in the initial stages, the eddy frequency is high. Implementation of an (instantaneous) eddy event further increases the strain rate within its interval, generating a turbulent cascade 46 0 5 10 15 20 25 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 τ < log10 χ > ODT−Z=0.42 DNS−Z=0.42 ODT−Z=0.55 DNS−Z=0.55 χ q Fig. 3.8: Evolution of log10 χ|Z with time (recall Zst = 0.42). The horizontal line indicates the steady extinction limit, χq. process (vortex stretching). For DNS, in the initial stages Kelvin-Helmholtz instabilities occur because of the high shear and significant velocity differences between fuel and oxidizer surface. Gradual growth in the size of coherent structures is observed for this case due to vortex pairing. Vortex pairing is inherently a multidimensional process that requires large structures at two different downstream locations to interact. The stand-alone ODT model cannot address this process because of its one-dimensional nature. This limitation of the model could help explain why early extinction occurs. As the simulation progresses, the strain rates become low and mixing rates relax. Once these rates relax the reignition takes place as described above. Overall the model exhibits stronger extinction and reignition characteristics compared to the DNS data. LES combined with ODT as a subgrid model may better capture both large-scale amalgamation as well as small-scale mixing processes representative of extinction in reactive jets. The present results may be compared with another recent effort to model this DNS flame that combined LEM with a three-dimensional LES and an artificial neural network approach 47 to accelerate chemistry computation [98]. The results obtained in the present work are of comparable quality to those obtained in [98], with a significantly reduced computational effort. Interestingly, both works exhibit the same key discrepancy with the DNS, i.e., the over-prediction of both the extinction and reignition processes. 3.5 Conclusions In this work, the ODT model is applied to a syngas jet flame, and direct comparison is made with DNS data. This study is first of its kind where a direct comparison has been made between ODT and 3D DNS data for a turbulent reacting flow. The present study focused on evaluating the model's ability to capture finite-rate-chemistry effects including extinction and reignition. A detailed comparison of jet spread rate as well as the thermochemistry for OH species has been presented. Results indicate that the ODT formulation can reproduce characteristics of the jet such as spread rate and entrainment. Additionally, the ODT model can qualitatively capture both extinction and reignition that are exhibited by the DNS data. The ODT calculations presented herein required approximately 2 hours per realization, and 400 realizations were used to provide well-converged statistics. Relative to DNS, ODT represents a very inexpensive modeling approach that can describe much of the physics present in the DNS, including PDF evolution, minor species evolution, finite-rate chemistry effects etc. Indeed, these results, together with the body of previous work in ODT of reacting flows [29,50,51], suggest that ODT can provide reasonably accurate predictions for turbulent combustion, even as a stand-alone model. CHAPTER 4 PARAMETER SENSITIVITY ANALYSIS 4.1 Introduction As described in Section 2.3, ODT model has number of adjustable parameters and simultaneous tuning of the parameters is needed to simulate a particular flow of interest. Presently, no theory or correlations exist to form the basis for parameter selection. In this chapter sensitivity analysis is performed to establish a common basis on which parameter values can be estimated. Two different configurations are chosen for the current study, turbulent nonreacting and reacting jets. 4.2 Turbulent Planar Jet Flame 4.2.1 Computational Configuration DNS of three-dimensional (3D) temporal planar syngas jet flames with detailed chem-istry over a range of jet Reynolds numbers (Rej) from 2510 to 9079 have been performed by Hawkes et al. [46-48]. Details of the DNS simulations (cases L, M and H) are summarized in Table 4.1. The number of grid points across slot width (D) is denoted by ND. Spatial (Δy) and temporal (Δt) resolutions are also given. The jet consists of a central fuel stream (50% CO, 10% H2 and 40% N2 by volume) surrounded by counter-flowing oxidizer streams comprised of 25% O2 and 75% N2. Table 4.1: Nonpremixed planar jet flame details [46]. Case D(mm) ND Rej U (m/s) tj = D U Δy Δt L 0.72 48 2510 144 5 μs 15 μm 2 ns M 0.96 64 4478 194 5 μs 15 μm 2 ns H 1.37 72 9079 276 5 μs 19 μm 2 ns 49 The fuel and oxidizer stream bulk velocities are U/2 and −U/2, respectively. The initial temperature of oxidizer stream is set to 500 K and pressure is set to 1 atm. For the chemical mechanism 11 species, H2, O2, O, OH, H2O, H, HO2, CO, CO2, HCO and N2, are considered with 21 reactions. For all of the cases considered here mixing is initially rapid enough relative to reaction to cause local extinction. Varying degrees of extinction are observed for the three cases. All three cases reignite following extinction, although at different rates. The dominant reignition mechanism for the cases considered here is turbulent folding [47]. The mixture fraction was computed from the local species compositions using Bilger's definition [8]. The scalar dissipation is computed based on the unity Lewis number (Le) assumption as χ = 2D Ä dZ dy ä2 with D = λ/ (ρcpLe). The stoichiometric mixture fraction is Zst = 0.42 and the steady extinction dissipation rate is χ = 2194 s−1. The transverse direction in Figure 3.1 is chosen as the ODT domain and initial condi-tions for all the variables are extracted directly from the DNS data. The detailed chemical mechanism considered in this study (consisting of 11 species and 21 reactions [48]), temper-ature and pressure-dependent thermodynamic property evaluation, and mixture-averaged diffusion are all consistent with the DNS simulations. The spatial and temporal resolution are likewise the same as used in the DNS simulation. The statistics from the DNS were extracted on xz planes in Figure 3.1. Figure 4.1 illustrates the temperature and χ/χq evolution with time in the form of contour plots for different DNS cases described in Table 4.1. These flames exhibit strong flame-turbulence interactions resulting in local extinction followed by reignition. The ex-tinction levels increase with increasing Rej . The maximum extinction is observed at τ = 20 for all three cases and is clearly evident from the low temperatures in Figure 4.1. The scalar dissipation rate χ is a quantity of critical importance in understanding and modeling non-premixed flames. As a measure of the local rate of molecular mixing, it plays a major role due to the intimate coupling of mixing and chemical reaction. It is heavily implicated in nearly every strategy to model turbulent combustion. 50 T 0 20 40 −4 −2 0 2 4 500 1000 1500 y D 0 20 40 −4 −2 0 2 4 500 1000 1500 τ 0 50 −4 −2 0 2 4 500 1000 1500 χ χq 0 20 40 −4 −2 0 2 4 0 0.5 1 0 20 40 −4 −2 0 2 4 0 0.5 1 τ 0 50 −4 −2 0 2 4 0.5 1 1.5 L M L M H H Fig. 4.1: DNS: Contour profiles of temperature (left) and scalar dissipation (right) for cases L, M and H. 51 The magnitude of χ/χq plays an important role in understanding the relative rates of mixing and chemical reactions. As Rej increases, χ/χq also increases and the maximum as seen in Figure 4.1 also peaks earlier for high Rej . Once the mixing rates relax all three flame reignite. However the onset of reignition event is different. Case L has largely reignited at τ = 40, whereas cases M and H are not fully reignited yet. Case H does not reignite until τ = 50. The greater degree of extinction at higher Rej , and therefore lower mean temperature and radical pool, account for the longer reignition times. Initial estimate on the model parameters indicated that the model behavior is quite sensitive to the choice of parameter C, and much less sensitive to the values α, β and Z. In Chapter 3, it is mentioned that the statistics are also sensitive to choice of β. However in the present study sensitivity analysis is restricted to the choice of model parameter C. The model parameter C determines the strength of the turbulence in ODT. Low values of C give a low rate of occurrence of eddies and consequently almost no eddies are imple-mented. In other words, when C is small enough, the flow is laminar. On the other hand, large values of C produce a lot of eddies, and thus the flow is very turbulent. The value of C is varied between 10 − 100. Strictly speaking 3 different C values are chosen, 10, 60 and 100. The other parameter values are α = 0.5, β = 1.0 and Z = 50 same for all the ODT simulations. The following procedure is employed to identify the influence of the C on the model behavior and also its dependency on the flow properties, • Check if a unique C value can reproduce the DNS statistics across Rej . • If a universally applicable C value is not identified, run simulations for cases L, M and H at different C values. • Identify the individual values for C that can produce qualitative agreement between the model and the data. • From the individual values identified for each of the L, M and H cases, deduce an empirical relationship between the parameter C and Rej . 52 • Check if the empirical correlation identified, when applied across Rej , can reproduce the DNS statistics, if not modify the coefficients of the correlation and identify the influence of them on the model outcome. 4.2.2 Results and Discussion Simulations are run for 0.23 ms, 0.25 ms and 0.3 ms for cases L, M and H, respectively. All the results reported here are analyzed over 400 ODT simulation realizations. The ODT calculations presented herein required approximately 2 h per realization. 4.2.2.1 Flow Entrainment (constant C) Figures 4.2, 4.3 and 4.4 show the comparison between the model behavior and DNS data of the mean streamwise velocity and mixture fraction evolution for cases L, M and H respectively. The following observations can be made when compared with DNS data: • For all the cases results indicate that changing C from 10 to 100 has little to no effect on the evolution of u and Z at τ = 6. • Case L: For C = 10, mixing is underpredicted at τ = 20 and 40, which is evident from the low spread and also slow decay of the u and Z . • Case L: For both C = 60 and C = 100, profiles match well at all the time intervals. • Case M: At τ = 20, mixing is underpredicted for C = 10 and overpredicted for both C = 60 and 100 (evident from fast decay in Z ). • Case M: The trend for over prediction in the mixing continues for C = 100 even at τ = 40, whereas mixing is in good agreement for C = 60. • Case H: For both C = 60 and 100 mixing is continuously over predicted (evident from fast decay in the u and Z and also low spreading). • Case H: At τ = 20, model shows very good agreement for C = 10. However underpre-dicts at τ = 40. 53 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.2: Case L: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). 54 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.3: Case M: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). 55 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.4: Case H: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). 56 4.2.2.2 Conditional Statistics (constant C) Figures 4.5, 4.6 and 4.7 show the comparison between the model behavior and DNS data of the mean temperature and hydroxyl radical, as a function of mixture fraction, for cases L, M and H, respectively. As described in Section 4.2.1 for all the cases simulated, mixing is initially rapid enough relative to reaction to cause local extinction, which is followed by reignition as mixing rates relax. The following observations can be made when compared with DNS data: • Case L: There is no difference in the way model behaves for C = 60 and 100. Both T|Z and OH|Z compare well with data, except at τ = 40 where strong reignition is observed. • Case L: At τ = 6, low temperature and OH values are reported for all values of C. However they are still above the extinction limit predicted by the laminar flamelet solution, indicating no extinction. • Case L: At τ = 20, for low value of C model predicts higher values of both T|Z and OH|Z . For medium and high values of C profiles fall below the extinction limit, indicating extinction. • Case L: At τ = 40, strong reignition is observed for all values of C. • Case M: At τ = 6, model indicates extinction for both C = 60 and 100. • Case M: At τ = 20, for all values of C the profiles fall below the extinction limit indicating extinction. The extinction event is stronger for C = 60 and 100. • Case M: At τ = 40, strong reignition is observed for all vales of C. Profiles of both T|Z and OH|Z are above the values predicted by the steady flamelet model at χq, indicating that reignition has occurred. • Case H: At τ = 6, model indicates extinction for all the values of C, i.e. profiles fall below the extinction limit. • Case H: At τ = 20, global extinction is observed for both C = 60 and 100. The temperature profiles are well below the extinction limit and the OH|Z is almost zero. For low value of C, the flame sustains but the values are still low compared to data. 57 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.5: Case L: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). The steady flamelet solution at χq is also shown for reference (red line). 58 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.6: Case M: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). The steady flamelet solution at χq is also shown for reference (red line). 59 0 0.2 0.4 0.6 0.8 1 500 1000 1500 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 500 1000 1500 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 500 1000 1500 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.7: Case H: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (C = 10), medium gray (C = 60), dark gray (C = 100). The steady flamelet solution at χq is also shown for reference (red line). 60 • Case H: At τ = 40, reignition is not observed for C = 60 and 100, where the profiles match the data for C = 10. 4.2.2.3 Empirical Correlation One of the main objectives of this study is to find if a unique value of C can reproduce the DNS statistics over a range of Rej . From the results discussed so far, in Sections 4.2.2.1, 4.2.2.2, it is clearly evident that C value needs to be changed based on Rej . Overall, for cases L, M and H, the selected model parameter values are 100, 60 and 10 respectively. These values are selected based on the following criteria. • The selected C value should reproduce the qualitative trends for extinction and reig-nition. For case H, only C = 10 reproduced both extinction and reignition. • If more than one value of C reproduces the qualitative behavior for both extinction and reignition, a value of C which quantitatively shows good agreement with data for both mixing and thermochemistry should be selected. The common notion of the ODT modeling community is that C value should be increased to increase the turbulence intensity (increases the number of eddies being implemented). However a reverse trend is observed in the current study (with increasing Rej , C decreases). Eddy events selection in ODT is a stochastic process and the acceptance/rejection of the candidate eddy depends on both the shear kinetic energy and the model parameter (see Section 2.3 for more details). For the cases described here the magnitude of gradients increase as the Rej increases. For case H, the gradients are so high that even a value of C = 10 reproduces the data. Based on the individual C values identified for all the cases the following empirical correlation is developed. Least squares regression is applied to identify the coefficients. Three different sets of coefficients are considered and included in Table 4.2 along with the residual values (R2). log C = logb + a log Rej , (4.1) 61 Table 4.2: Proposed empirical correlation (4.1) coefficients, a and b, values. Residual values are also given. set a log b R2 1 -1.49 15.11 0.834 2 -1.69 17.4 0.652 3 -1.69 16.9 0.768 When the correlation is considered for the simulations, C changes during the course of the simulation based on the local Rej which is computed from ujD/ν, where uj is the difference between maximum and minimum velocities and ν is the kinematic viscosity. The following sections focus on the comparison of model behavior ,with data, for different sets proposed in Table 4.2. 4.2.2.4 Flow Entrainment (variable C) Figures 4.8, 4.9 and 4.10 show the comparison between the model behavior and DNS data of the mean streamwise velocity and mixture fraction evolution for cases L, M and H, respectively. The following observations can be made when compared with DNS data: • Case L: At τ = 6, for all the different sets of coefficients considered here model behavior shows good agreement with the data. • Case L: At τ = 20 and 40, the model overpredicts the decay and spreading of both u and Z . • Case M: At τ = 6, the model behavior shows good agreement with the data. • Case M: At τ = 20 and 40, the model overpredicts the decay and spreading of both u and Z . However the set-1 is in close agreement with the data compared to other two cases. • Case H: At τ = 6, the model behavior shows good agreement with the data. • Case H: At τ = 20, model indicates the same behavior for set-1 and set-3. The spreading matches well but the velocity at the jet center is high. Both spreading and decay of the centerline velocity matches well for set-2. • Case H: At τ = 40, model indicates high mixing for set-2, whereas set-1 and set-3 behaves much the same way. 62 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 <u> [m/s] y/D −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −80 −60 −40 −20 0 20 40 60 80 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.8: Case L: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). 63 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −100 −50 0 50 100 <u> [m/s] y/D −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.9: Case M: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). 64 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (a) τ = 6 −4 −2 0 2 4 −100 −50 0 50 100 <u> [m/s] y/D −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 <Z> y/D (b) τ = 20 −4 −2 0 2 4 −100 −50 0 50 100 y/D <u> [m/s] −4 −2 0 2 4 0 0.2 0.4 0.6 0.8 1 y/D <Z> (c) τ = 40 Fig. 4.10: Case H: Mean velocity (left) and mixture fraction (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). 65 4.2.2.5 Conditional Statistics (variable C) Figures 4.11, 4.12 and 4.13 show the comparison between the model behavior and DNS data of the mean temperature and hydroxyl radical, as a function of mixture fraction, for cases L, M and H, respectively. The following observations can be made when compared with DNS data: • Case L: At all the time intervals the model behaves the same way for different sets of coefficients considered here. • Case L: At τ = 6, model reports low temperature and OH values. However they are still above the extinction limit. • Case L: At τ = 20, profiles fall below the extinction limit indicating extinction. Com-pared to data, higher extinction is reported in the simulation. • Case L: At τ = 40, strong reignition is observed in the model behavior. The values predicted by the model are higher compared to DNS data. • Case M: There is little to no difference in the way model behaves for different sets of coefficients considered here. • Case M: At τ = 6, model reports low temperature and OH values. However they are still above the extinction limit. • Case M: At τ = 20, for all values of C the profiles fall below the extinction limit indicating extinction. The extinction event is stronger in the model. • Case M: At τ = 40, reignition is observed in the model. Profiles of both T|Z and OH|Z are above the extinction limit indicating that reignition has occurred and also compare well with the data. • Case H: At τ = 6, model indicates early extinction for all the sets which is not observed in the data. • Case H: At τ = 20, profiles fall below the extinction limit for all the sets indicating extinction. For set-1 and set-3 the model shows good agreement with the data. • Case H: At τ = 40, reignition is not observed for set-2. The profiles for set-1 and set-3 match well with the data. 66 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 <T> [K] Z 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.11: Case L: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). The steady flamelet solution at χq is also shown for reference (red line). 67 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 <T> [K] Z 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.12: Case M: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). The steady flamelet solution at χq is also shown for reference (red line). 68 0 0.2 0.4 0.6 0.8 1 500 1000 1500 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 Z <Y OH > (a) τ = 6 0 0.2 0.4 0.6 0.8 1 500 1000 1500 <T> [K] Z 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 <Y OH > Z (b) τ = 20 0 0.2 0.4 0.6 0.8 1 500 1000 1500 Z <T> [K] 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 x 10−3 Z <Y OH > (c) τ = 40 Fig. 4.13: Case H: Conditional mean temperature (left) and hydroxyl radical (right) profiles for different C values at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), medium gray (set-2), dark gray (set-3). The steady flamelet solution at χq is also shown for reference (red line). 69 4.2.2.6 PDF Evolution From the results discussed so far it is clearly evident that proposed correlation can reproduce the DNS statistics. The choice of the coefficients influence the mixing and ther-mochemistry behavior. The extinction and reignition observed in the model is reproduced by the model. However it is important to look at the higher order statistics like probability density function evolution (PDF) before drawing major conclusions about the model. In this section PDF evolution of temperature and scalar dissipation are compared with the DNS data. The PDFs are generated for both a constant value of C and correlation based C based with set-1 coefficients. The constant values of C are 100, 60 and 10 for cases L, M and H, respectively. Figures 4.14, 4.15 and 4.16 describe the PDF evolution of both T and log10(χ/χq), con-ditioned on stoichiometric mixture fraction (Zst = 0.42), for cases L, M and H, respectively. The following observations can be made from the comparison: • For all the cases (L, M and H), qualitatively, model indicates the same behavior for both constant and varying C. • Case L: At τ = 6, scalar dissipation PDF evolution shifted towards higher values relative to the DNS data. These higher values of χ cause extinction in the early stages of the ODT simulations, resulting in a corresponding temperature PDF shift towards lower values. • Case L: At τ = 20, mixing rates are still high enough to cause extinction in the model, as indicated by the χ PDF. The temperature keeps dropping as indicated by the PDF. • Case L: At τ = 40, mixing rates relax and temperature PDF evolution starts shifting towards high values as reignition occurs. • Case M: At τ = 6, mixing rate is high as indicated by the χ PDF shift towards higher values. The corresponding temperature PDF shifts towards lower values. • Case M: At τ = 20, the mixing rates are comparable to DNS. However the rates are still high enough to cause extinction in the model. 70 1200 1300 1400 1500 1600 0 0.05 0.1 0.15 0.2 0.25 T | Z=Z st [K] PDF −2 −1 0 1 2 0 0.05 0.1 0.15 0.2 0.25 PDF log 10 (χ/χ q ) | Z=Z st (a) τ = 6 500 1000 1500 2000 0 0.02 0.04 0.06 0.08 0.1 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.02 0.04 0.06 0.08 0.1 PDF log 10 (χ/χ q ) | Z=Z st (b) τ = 20 500 1000 1500 2000 0 0.05 0.1 0.15 0.2 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.05 0.1 0.15 0.2 PDF log 10 (χ/χ q ) | Z=Z st (c) τ = 40 Fig. 4.14: Case L: Probability Density Function profiles of temperature (left) and log10(χ/χq) (right), conditioned on stoichiometric mixture fraction (Zst = 0.42), at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), dark gray (C = 100). 71 800 1000 1200 1400 1600 0 0.02 0.04 0.06 0.08 0.1 T | Z=Z st [K] PDF −4 −2 0 2 4 0 0.05 0.1 0.15 0.2 PDF log 10 (χ/χ q ) | Z=Z st (a) τ = 6 500 1000 1500 0 0.02 0.04 0.06 0.08 0.1 0.12 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.02 0.04 0.06 0.08 PDF log 10 (χ/χ q ) | Z=Z st (b) τ = 20 800 1000 1200 1400 1600 0 0.02 0.04 0.06 0.08 0.1 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.02 0.04 0.06 0.08 0.1 PDF log 10 (χ/χ q ) | Z=Z st (c) τ = 40 Fig. 4.15: Case M: Probability Density Function profiles of temperature (left) and log10(χ/χq) (right), conditioned on stoichiometric mixture fraction (Zst = 0.42), at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), dark gray (C = 60). 72 500 1000 1500 2000 0 0.02 0.04 0.06 0.08 0.1 0.12 T | Z=Z st [K] PDF −4 −2 0 2 4 0 0.02 0.04 0.06 0.08 0.1 0.12 PDF log 10 (χ/χ q ) | Z=Z st (a) τ = 6 500 1000 1500 0 0.02 0.04 0.06 0.08 0.1 0.12 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.02 0.04 0.06 0.08 PDF log 10 (χ/χ q ) | Z=Z st (b) τ = 20 500 1000 1500 2000 0 0.02 0.04 0.06 0.08 0.1 T | Z=Z st [K] PDF −6 −4 −2 0 2 0 0.02 0.04 0.06 0.08 0.1 PDF log 10 (χ/χ q ) | Z=Z st (c) τ = 40 Fig. 4.16: Case H: Probability Density Function profiles of temperature (left) and log10(χ/χq) (right), conditioned on stoichiometric mixture fraction (Zst = 0.42), at τ = 6, 20 and 40. Dashed line (DNS), solid line (ODT), light gray (set-1), dark gray (C = 10). 73 • Case M: At τ = 40, mixing rates relax as indicated by the dissipation PDF shift towards lower values and the temperature PDF evolution starts shifting towards high values as reignition occurs. • Case H: At τ = 6, mixing rates are high in the model compared to data as indicated by the χ profiles. The corresponding temperature values are low. • Case H: At τ = 20, mixing rates are comparable to the data and corresponding temperature PDF profiles also show good agreement. • Case H: At τ = 40, mixing rates relax and the dissipation PDF compares well with the data. Interestingly temperature PDF shows bimodal distribution whereas data indicate a monomodal distribution. Overall, the model exhibits stronger extinction and reignition characteristics compared to the DNS data. The early extinction observed in the model is not directly explicable from mean profiles of u and Z . However the χ PDF profiles indicate that the mixing rates are high in the model compared to data. The high mixing rates in the early stages causes extinction in the model, see profiles of T|Z and OH|Z at τ = 6 for cases M and H. The reasons for high mixing rates in the early stage are discussed in Chapter 3. As the flame evolves in the time, the decay of mixing rates allow T|Z and OH|Z to move towards their equilibrium values. 4.3 Nonreacting Turbulent Planar Jet In this section the proposed correlation is applied to a different nonreacting configura-tion. 4.3.1 Computational Configuration Temporally developing planar jet configuration with air, at room temperature and pres-sure as the fluid, is simulated to validate the model. The initial conditions for simulating planar jet are given in Table 4.3. The streamwise velocity at the inlet is specified using the following hyperbolic tangent function and is shown schematically in Figure 4.17. 74 Table 4.3: Initial conditions of different Rej cases defined for parameter sensitivity analysis. Rej D (m) u0 (m/s) u∞ (m/s) dt (s) dy (m) 2250 0.002 28 10 2e-7 100e-6 5000 0.003 37 10 2e-7 100e-6 9000 0.004 46 10 2e-7 100e-6 14000 0.005 55 10 2e-7 100e-6 27500 0.007 73 10 2e-7 100e-6 36000 0.008 82 10 2e-7 100e-6 Fig. 4.17: Schematic of tanh profile used to specify streamwise velocity profile. ϕ(y) = A 2 ï 1 + tanh Å y − L1 w ãò ß 1 − 1 2 ï 1 + tanh Å y − L2 w ãò™ , (4.2) where A is the amplitude of the change w is the width of the transition and L1 and L2 are the midpoints of the transition. 4.3.2 Experimental Data Following semiempirical relation developed by Gutmark [40], for center line velocity (uc) decay of spatially developing planar jet, is used to compare with simulation data, 75 uc = u∞ + u∗ 0.188(x−x0) b , (4.3) where u∗ = » u0(u0 − u∞), offset (x0) near the jet nozzle depends on the jet velocity con-ditions, for the current study a value of −4D is considered. To compare the temporally evolving flow simulation results with spatially evolving jet experimental data, the procedure described in Section 2.2.1.3 is followed. 4.3.3 Results and Discussion Figure 4.18 shows the comparison of mean centerline velocity evolution between sim-ulation and experimental data. Simulations for the low Rej range (2250-9000) performed only using the empirical correlation. For high Rej regime, simulations are performed using a constant value of C and also the corrleation based C. As it can be clearly seen from the comparison, for Rej ranging from 2250-9000, the model reproduces the experimental data. The proposed correlation for C underpredicts the velocity decay for high Rej (14000-36000), whereas a constant C value (C = 10) reproduces the experimental data. For the proposed correlation, in the limit of Rej →∞, C goes to zero. Model parameter, C determines the turbulence strength in the flow, if C value is very low fewer eddies will be implemented and the f