Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Effects of cold wall quenching on unburned hydrocarbon emissions from a natural gas HPDI engine Turcios, Marco Antonio 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2011_fall_turcios_marco.pdf [ 10.34MB ]
Metadata
JSON: 24-1.0072348.json
JSON-LD: 24-1.0072348-ld.json
RDF/XML (Pretty): 24-1.0072348-rdf.xml
RDF/JSON: 24-1.0072348-rdf.json
Turtle: 24-1.0072348-turtle.txt
N-Triples: 24-1.0072348-rdf-ntriples.txt
Original Record: 24-1.0072348-source.json
Full Text
24-1.0072348-fulltext.txt
Citation
24-1.0072348.ris

Full Text

Effects of Cold Wall Quenching on Unburned Hydrocarbon Emissions from a Natural Gas HPDI Engine by Marco Antonio Turcios BASc., The University of British Columbia, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2011 © Marco Antonio Turcios 2011 Abstract The quenching of hydrocarbon flames on cold surfaces is considered to be a potential source of unburned hydrocarbon (UHC) emissions from inter- nal combustion engines, but its contribution to emissions has been difficult to determine due to the strong coupling between physics, chemistry and flame/wall geometry. This is particularly problematic for high pressure di- rect injection (HPDI) engines where high pressures, inhomogeneous mixtures and complex piston geometry are present. In this work, a computational model is implemented to determine the distance at which hydrocarbon flames quench on cold walls during numerical simulation. This model accounts for variable pressure, temperature, gas mixture and the geometry conditions. The model presented in this work is an extension of the experimental work done by Boust et. al. with stoichiometric premixed flames at low pressures. The validation of this model for high pressure and diffusion flames is presented and shows that the correct trends in heat flux and order of magnitude of quench distance are observed. This model is further refined for engine simulation and enhanced by a two-zone mass diffusion model to account for post-quench oxidation of boundary fuel. A selection of engine cases are simulated for a variety of different con- ditions to determine the spatial distribution and temporal evolution of un- burned fuel cold surfaces. It was found that wall quenching on the piston contributes up to 50% of the total UHC during the combustion cycle, the ma- jority of which is oxidized during the expansion stroke; the final contribution is at most 10% but frequently near or less than 1%. As the injection pres- sure was increased, quenching on the piston surface became more extensive, through the quenching thickness itself decreased. UHC from wall quenching occurs more readily for higher load conditions due to the richer mixtures and incomplete mixing. Altered engine timing introduced coupled effects of changed flame/wall interaction and combustion characteristics. The data obtained from the model can be used to evaluate attempts to reduce UHC by changing combustion chamber geometry. ii Preface The majority of research contained in this report was undertaken by the au- thor. Chapter 2 and Chapter 4, Section 3.2 contain work that appears in the proceedings of the 2011 JSAE/SAE International Powertrains, Lubricants and Fuels meeting in Kyoto, Japan (JSAE 20119123, SAE 2011-01-1997). Drs. Jim Huang, Sandeep Munshi and Gordon McTaggart-Cowan at West- port Innovations were involved in the preparation and proofing of the final paper, with input from Dr. Carl Ollivier-Gooch of the University of British Columbia. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 The Problem of Unburned Hydrocarbons . . . . . . . . . . . 1 1.2 Goals and Overview . . . . . . . . . . . . . . . . . . . . . . . 2 2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Wall Quenching Terminology . . . . . . . . . . . . . . . . . . 3 2.2 Thermal Considerations . . . . . . . . . . . . . . . . . . . . . 4 2.2.1 Stoichiometric Laminar Premixed Flames . . . . . . . 4 2.2.2 Strained Diffusion Flames . . . . . . . . . . . . . . . 6 2.2.3 Flame Wall Interaction at a Distance . . . . . . . . . 10 2.2.4 Turbulence and Combustion Coupling . . . . . . . . . 10 2.3 Chemistry Considerations . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Chemical Mechanisms . . . . . . . . . . . . . . . . . . 12 2.3.2 Post-quench Diffusion . . . . . . . . . . . . . . . . . . 15 2.4 Modelling Considerations . . . . . . . . . . . . . . . . . . . . 16 iv Table of Contents 2.4.1 Assumptions and Simplifications . . . . . . . . . . . . 16 2.4.2 Phenomenological - Quench distance calculation . . . 16 2.4.3 Physical - Flame structure modification . . . . . . . . 19 2.4.4 Model Selection . . . . . . . . . . . . . . . . . . . . . 20 3 Model Development and Validation . . . . . . . . . . . . . . 21 3.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 High Pressure Premixed Flames . . . . . . . . . . . . . . . . 21 3.2.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.2 Flame Speed . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.3 Heat Flux . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.4 Quenching Distance . . . . . . . . . . . . . . . . . . . 29 3.3 Strained Diffusion Flames . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Initial and Boundary Conditions . . . . . . . . . . . . 32 3.3.2 Heat Flux Calculation . . . . . . . . . . . . . . . . . . 34 3.3.3 Flame Power Calculation . . . . . . . . . . . . . . . . 36 3.3.4 Simulation Results . . . . . . . . . . . . . . . . . . . . 38 4 Model in Engine Simulation . . . . . . . . . . . . . . . . . . . 41 4.1 Overview of Algorithm . . . . . . . . . . . . . . . . . . . . . 41 4.2 Validation under Engine Conditions . . . . . . . . . . . . . . 41 4.2.1 Order of Accuracy . . . . . . . . . . . . . . . . . . . . 44 4.3 Natural Gas HPDI . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3.1 Compression Ignition and HPDI . . . . . . . . . . . . 44 4.3.2 Westport HPDI . . . . . . . . . . . . . . . . . . . . . 45 4.4 Flame Power Variation . . . . . . . . . . . . . . . . . . . . . 46 4.4.1 Density . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.4.2 Flame Speed . . . . . . . . . . . . . . . . . . . . . . . 47 4.5 Peb and δl Calculation and Filtering . . . . . . . . . . . . . 59 4.6 Diffusion/Oxidation Model . . . . . . . . . . . . . . . . . . . 61 5 Engine Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1 Engine Conditions . . . . . . . . . . . . . . . . . . . . . . . . 67 5.2 Geometry and Mesh . . . . . . . . . . . . . . . . . . . . . . . 67 5.3 Engine Simulation Results . . . . . . . . . . . . . . . . . . . 69 5.3.1 Comparison with Experiments . . . . . . . . . . . . . 69 5.3.2 Numerical Quenching Results . . . . . . . . . . . . . 70 5.3.3 Mode 8 Results . . . . . . . . . . . . . . . . . . . . . 73 5.3.4 Mode 9 Cases . . . . . . . . . . . . . . . . . . . . . . 81 5.4 Engine Timing Changes . . . . . . . . . . . . . . . . . . . . . 97 v Table of Contents 5.4.1 Normal Injection Pressure . . . . . . . . . . . . . . . 97 5.4.2 High Injection Pressure . . . . . . . . . . . . . . . . . 110 5.5 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . 115 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.1 Future Work and Applications . . . . . . . . . . . . . . . . . 120 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 A Thermal Quench Model . . . . . . . . . . . . . . . . . . . . . . 128 B Numerical Model Details . . . . . . . . . . . . . . . . . . . . . 131 C Quench Model Code . . . . . . . . . . . . . . . . . . . . . . . . 140 D Refinements to Model . . . . . . . . . . . . . . . . . . . . . . . 146 vi List of Tables 3.1 Experimental conditions . . . . . . . . . . . . . . . . . . . . . 23 3.2 Peak heat fluxes, head on quenching . . . . . . . . . . . . . . 27 3.3 Peak heat fluxes, sidewall quenching . . . . . . . . . . . . . . 28 3.4 k −  wall function summary . . . . . . . . . . . . . . . . . . 36 3.5 Gülder laminar flame model . . . . . . . . . . . . . . . . . . . 38 4.1 Validation mode statistics . . . . . . . . . . . . . . . . . . . . 43 4.2 RNG k −  wall function summary . . . . . . . . . . . . . . . 44 5.1 Engine cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Phase breakdown . . . . . . . . . . . . . . . . . . . . . . . . . 69 vii List of Figures 2.1 Quenching quantities . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Quenching of strained diffusion flames . . . . . . . . . . . . . 9 2.3 Quenching of turbulent premixed flames . . . . . . . . . . . . 13 2.4 Species profiles through a premixed flame . . . . . . . . . . . 14 2.5 Post-quench diffusion of boundary fuel . . . . . . . . . . . . . 15 2.6 Unburned hydrocarbons from pre-chamber engine . . . . . . . 17 2.7 Thermal quenching model, head on quenching . . . . . . . . . 19 3.1 Computational domain . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Convergence of temperature . . . . . . . . . . . . . . . . . . . 24 3.3 Laminar flame speeds . . . . . . . . . . . . . . . . . . . . . . 26 3.4 Methane head on quench distance . . . . . . . . . . . . . . . 30 3.5 Methane sidewall quench distance results . . . . . . . . . . . 30 3.6 Heptane quench distance results . . . . . . . . . . . . . . . . 31 3.7 Initial conditions for strained diffusion flame . . . . . . . . . . 33 3.8 Quenching model results for strained diffusion flames . . . . . 39 3.9 Quench distance model . . . . . . . . . . . . . . . . . . . . . . 40 4.1 Engine condition map . . . . . . . . . . . . . . . . . . . . . . 42 4.2 Location of Peclet number sampling . . . . . . . . . . . . . . 43 4.3 Westport fuel system diagram . . . . . . . . . . . . . . . . . . 46 4.4 QΣ related contours . . . . . . . . . . . . . . . . . . . . . . . 49 4.5 Density variation at quenching . . . . . . . . . . . . . . . . . 51 4.6 φ− T analysis of Sl . . . . . . . . . . . . . . . . . . . . . . . 52 4.7 Local flame power and heat flux comparison . . . . . . . . . . 53 4.8 Peb from local variables . . . . . . . . . . . . . . . . . . . . . 55 4.9 Stoichiometric flame power and heat flux comparison . . . . . 56 4.10 Peb based on stoichiometric flame speed . . . . . . . . . . . . 58 4.11 Peb filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.12 Variation of Peclet number with location . . . . . . . . . . . . 62 4.13 Schematic of post-quench diffusion model . . . . . . . . . . . 64 4.14 Comparison of cell center and quench layer values of YCH4 . . 65 viii List of Figures 5.1 Computational domain . . . . . . . . . . . . . . . . . . . . . . 69 5.2 Fuel consumption and power comparison . . . . . . . . . . . . 71 5.3 Emissions comparison . . . . . . . . . . . . . . . . . . . . . . 72 5.4 Mode 8 combustion results . . . . . . . . . . . . . . . . . . . 74 5.5 Mode 8 Pinj = 250 bar flame propagation . . . . . . . . . . . 76 5.6 Mode 8 Pinj = 250 bar Peb distribution . . . . . . . . . . . . 77 5.7 Mode 8 Pinj = 250 bar global quenching . . . . . . . . . . . . 79 5.8 Mode 8 Pinj = 250 bar global quench distance . . . . . . . . . 80 5.9 Mode 8 Pinj = 600 bar flame propagation . . . . . . . . . . . 82 5.10 Mode 8 Pinj = 600 bar Peb distribution . . . . . . . . . . . . 83 5.11 Mode 8 Pinj = 600 bar global quenching . . . . . . . . . . . . 85 5.12 Mode 8 Pinj = 600 bar global quench distance . . . . . . . . . 85 5.13 Mode 9 combustion results . . . . . . . . . . . . . . . . . . . 86 5.14 Mode 9 Pinj = 150 bar flame propagation . . . . . . . . . . . 88 5.15 Mode 9 Pinj = 150 bar Peb distribution . . . . . . . . . . . . 89 5.16 Mode 9 Pinj = 150 bar global quenching . . . . . . . . . . . . 91 5.17 Mode 9 Pinj = 150 bar global quench distance . . . . . . . . . 91 5.18 Mode 9 Pinj = 250 bar flame propagation . . . . . . . . . . . 93 5.19 Mode 9 Pinj = 250 bar Peb distribution . . . . . . . . . . . . 94 5.20 Mode 9 Pinj = 250 bar global quenching . . . . . . . . . . . . 96 5.21 Mode 9 Pinj = 250 bar global quench distance . . . . . . . . . 96 5.22 Mode 8 altered timing combustion results . . . . . . . . . . . 98 5.23 Mode 8 Pinj = 250 bar τa flame propagation . . . . . . . . . . 99 5.24 Mode 8 Pinj = 250 bar τa Peb distribution . . . . . . . . . . . 100 5.25 Mode 8 Pinj = 250 bar τa global quenching . . . . . . . . . . 102 5.26 Mode 8 Pinj = 250 bar τa global quench distance . . . . . . . 103 5.27 Mode 8 Pinj = 250 bar τr flame propagation . . . . . . . . . . 105 5.28 Mode 8 Pinj = 250 bar τr Peb distribution . . . . . . . . . . . 106 5.29 Mode 8 Pinj = 250 bar τr global quenching . . . . . . . . . . 108 5.30 Mode 8 Pinj = 250 bar τr global quench distance . . . . . . . 108 5.31 Mode 8 high injection pressure, altered timing combustion results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.32 Mode 8 Pinj = 600 bar τr flame propagation . . . . . . . . . . 111 5.33 Mode 8 Pinj = 600 bar τr Peb distribution . . . . . . . . . . . 112 5.34 Mode 8 Pinj = 600 bar τr global quenching . . . . . . . . . . 114 5.35 Mode 8 Pinj = 600 τr bar global quench distance . . . . . . . 115 5.36 Unburned hydrocarbon comparison . . . . . . . . . . . . . . . 116 5.37 Quench model comparison . . . . . . . . . . . . . . . . . . . . 117 A.1 Simplified model of head-on quenching . . . . . . . . . . . . . 129 ix List of Figures D.1 Ypg1 and QW at quenching locations . . . . . . . . . . . . . . 148 x List of Programs 4.1 Algorithm for computing quench distance . . . . . . . . . . . 41 B.1 Code for solving Navier-Stokes equations . . . . . . . . . . . . 132 B.2 fvSchemes file . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 B.3 fvSolution file . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 C.1 qDist.H . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 C.1 updateBoundary.H . . . . . . . . . . . . . . . . . . . . . . . . 144 xi Nomenclature Note: ( ) used as a placeholder for sub- and superscripts. ( )B Quantity at boundary ( )st Quantity at stoichiometry ( )t Turbulent quantity ( )u Properties of fresh (unburnt) gases ( )W Quantity at wall α Thermal diffusivity ∆H Heat of reaction δl Thermal flame thickness δq Quenching distance  Turbulent kinetic energy dissipation rate κ Thermal conductivity µ Dynamic viscosity ν Kinematic viscosity φ Equivalence ratio ρ Gas density ζ1 Conserved scalar for methane, with φ = ζ1 ζ1,st D Mass diffusivity h Enthalpy k Turbulent kinetic energy xii Nomenclature Peb Burning Peclet number: δq α/Sl Pr Prandtl number: να Q Heat flux QΣ Flame power Re Reynolds number:ρULµ , where U and L are the representative velocity and length scale (respectively) of the flow Sl Laminar flame speed Sc Schmidt number: µρD tf Flame time: α S2l Yi Mass fraction of species i Ypg Reaction progress variable ATDC After top dead center CFD Computational fluid dynamics ODE Ordinary differential equation PDE Partial differential equation TDC Top dead center UHC Unburned hydrocarbon xiii Acknowledgments Much appreciation goes to Drs. Jim Huang, Sandeep Munshi and Phillip Hill for giving me the opportunity to work with such an interesting project at Westport Innovations. Thanks in particular to Drs. Jim Huang and Carl Ollivier-Gooch for being excellent supervisors and colleagues. This research was funded by NSERC IPS in cooperation with Westport Innovations. This research has also been aided by the use of computing resources provided by WestGrid and Compute/Calcul Canada. xiv Dedication I would like to personally thank my family for their support and encourage- ment, without them I could not have persevered and accomplished such an undertaking. xv Chapter 1 Introduction 1.1 The Problem of Unburned Hydrocarbons Internal combustion engines that use natural gas as fuel are seen as a viable alternative in several industries hoping to reduce the amount of engine emis- sions produced during power generation while maintaining (or increasing) the same level of efficiency and performance. Natural gas is abundant and currently available at a lower cost than other hydrocarbon fuels, and has already been implemented in several transportation sectors. The influence of engine emissions on global climate and human health are still matters of scientific debate and research[39], but government and regulatory bodies are already taking steps to set acceptable limits on greenhouse gas emissions[56]. These initiatives, particularly those concentrating on methane emissions, re- sult in the the requirement that very little unburned hydrocarbons (UHC) remain after combustion. One possible avenue for UHC reduction is through the use of catalytic converters, though the conversion efficiency rapidly decreases to low levels (less than 20 %) over timescales too short for engine use; there are also additional difficulties posed by the stability of methane, the high exhaust temperature required for efficient conversions and the sensitivity of catalysts to contaminants[32]. The addition of another sub-system increases the cost and complexity of the engine. Another method is to ensure that in-cylinder combustion is as complete as possible. The difficulty lies in separating the results of different inter-related physical phenomena during the combustion process for any fuel[4, 38]. The main ways in which UHC is thought to be produced in the turbulent diffusion flames encountered in high pressure di- rect injection engines are incomplete mixing of fuel and oxidizer, incomplete reaction of fuel due to low temperature, the quenching of flames on cold walls and the quenching of flames due to excessive strain. To optimize the design process of internal combustion engines with re- spect to UHC emissions, it is important to know the relative contributions of each source of UHC, and the influence that parameters such as geometry, timing and injection pressure may have on emissions. Numerical simulation 1 1.2. Goals and Overview is frequently used in the early stages of engine design to test a variety of different concepts and determine important trends that can then be verified through experimentation. Unfortunately, the length scale associated with wall quenching (∼ 100 µm) prohibits direct simulation in engine cases unless some sort of modeling is used. 1.2 Goals and Overview The objective of this study is to develop a model for wall quenching that can be applied to general combustion simulation to obtain an estimate of the amount of UHC that remains in the quenching layer during combustion in a high pressure direct injection (HPDI) natural gas engine. The model should be able to account for changes in temperature, pressure, fuel/oxidizer mixtures and geometric configurations. A survey of the relevant literature will identify the main characteristics of wall quenching and determine currently existing models, as well as their suit- ability for engine simulation. From this review, a thermal model is selected and evaluated to see if it is applicable to the high pressure flames present in the engines under consideration. An implementation of this model will be evaluated by comparing the results from the simulation of strained diffusion flames with recent numerical simulation, as the nature of flames in HPDI combustion is different than premixed quiescent stoichiometric flames. After the model has been shown to provide consistent information about the quenching distance for high pressure diffusion flames, it will be further refined by examining the results for an HPDI engine at high load. Particular attention will be paid to the inhomogeneity of the fuel-air mixture and its effects on model results. Refinements to the model will be made to calculate the quenching distance over all the cold walls of the combustion chamber, and a mass exchange model that accounts for the post-quench diffusion and oxidation of fuel trapped in the quench layer will also be implemented. Once testing of the model is complete, several test cases are presented that vary the injection pressure, timing and load to determine what effect they have on the amount of UHC present for a single cycle. The temporal and spatial evolution of UHC over the piston surface is also observed, as variation of the the aforementioned parameters will change the interaction between the flame and the piston. This information will be useful to en- gine designers wishing to minimize the amount of UHC from wall quenching without negatively affecting performance in other areas, and those contem- plating new piston geometries to improve combustion characteristics. 2 Chapter 2 Literature Review 2.1 Wall Quenching Terminology In early studies on the visualization of combustion inside internal combustion engines, it was noticed that there were dark regions (indicating low temper- ature) between the propagating flame and the cold walls of the chamber[14]. It was speculated that this region could not support combustion due to the low gas temperature preventing flame from propagating (because of enthalpy losses), and that this region was contributing significant amounts of UHC to exhaust gas emissions. Sampling valve measurements at first seemed to confirm this idea[1], but subsequent tests with an improved valve[26], and other experiments seemed to indicate that although combustion is indeed inhibited near cold walls, the majority of unburned fuel is subsequently oxidized[9, 43, 54]. This analysis applies mostly to the extended surfaces of the piston crown and the cylinder walls; there are still significant contribu- tions to UHC from crevice volumes that remain fairly constant over a range of equivalence ratios. To compare quenching under different conditions it is useful to non- dimensionalize quench distance δq ; the most commonly accepted way is by using the diffusive flame thickness[40], defined as δl = α/Sl . This gives a burning Peclet number Peb = δqSl α = Sl α/δq (2.1) where Sl is the laminar flame speed and α is the thermal diffusivity of the gas. In its second form in Equation 2.1, Peb gives a measure of how the propagation speed of the flame and thermal diffusion play key roles in determining quenching behaviour. To observe wall quenching under more controlled and easily observable conditions, studies have focused on single wall quenching, where a planar flame either impinges normally on a cold wall (know as head-on quenching) or propagates parallel to a cold wall (side-wall quenching). Both numer- ical and experimental studies have studied the effect cold walls have on 3 2.2. Thermal Considerations chemistry[20, 54] and heat transfer [9]. For a variety of conditions, it was found that the Peclet number was primarily dependent on geometry with values of ∼3 for head-on-, ∼8 for sidewall- and ∼50 for tube quenching[4, 41] (typically at a gas and wall temperature of 300 K and pressure of 1 bar). A higher Peb can be seen as both an indication that δq is larger than δl and that the speed of flame propagation has more influence on that particu- lar quenching configuration. Quenching on the combustion chamber wall is similar in behaviour to single wall quenching, whereas quenching in crevices is more closely associated with two-wall or tube quenching. Two quantities are used to quantify the heat transfer processes that occur during wall quenching. The first is the heat flux from the flame to the wall, which is calculated as QW = −k∇T · n̂, (2.2) where k is the thermal conductivity of near-wall gases, T is the temperature and n̂ is the outward normal at the wall. As ∇T will point away from the wall and back into the flame at quenching (since the wall temperature is assumed to be lower than the burning mixture), QW is positive when the hot flame approaches the cold wall. Related to the wall heat flux is a quantity known as the laminar flame power, with several different definitions. Two of the most common are QΣ = ρSlcp∆T (2.3) QΣ = ρSl∆HYfuel (2.4) where ρ is the density of the gas mixture, cp is the isobaric heat capacity, ∆T is the change in temperature across the flame front, ∆H is the heat of combustion and Yfuel is the mass fraction of fuel present in the unburned gases. In both forms, QΣ is used to represent the maximum heat release possible from the flame. Flame power is frequently used to normalize QW , Φ = QW QΣ . (2.5) This is not to be confused with the equivalence ratio φ. 2.2 Thermal Considerations 2.2.1 Stoichiometric Laminar Premixed Flames As mentioned earlier the contribution of wall quenching to UHC emissions was initially overestimated. Subsequently studies that improved understand- 4 2.2. Thermal Considerations ing of flame quenching in internal combustion engines have focused on sin- gle wall quenching, particularly on two geometric configurations: head-on quenching where the flame is parallel to the wall and is propagating normal to the wall, and sidewall quenching where the flame is perpendicular to the wall and propagating in a direction parallel to it. Studies have focused on premixed flames at lower temperatures and pressures since they are more easily controlled and observed, though there has been recent work at higher pressures[27]. Experimental and numerical investigations have quantified the size of the δq , and although the exact number fluctuates depending on the criteria used to define the quenching distance (chemical species concentration, smallest flame-wall distance, etc)[9, 10, 13, 54], the order of magnitude is typically less than 1 mm for methane mixtures, decreasing with increasing wall tem- perature and pressure, and increasing with lower equivalence ratios. Few studies have made comparison with other fuels, but Westbrook’s seminal numerical study compared quenching between methane and methanol and found the quench distance for methane to be double that of methanol[54]. As it is believed that hotter flames (due to fuel type, pressure or equiva- lence ratio) are thinner and thus able to propagate closer to a cold wall, this behaviour is consistent. Measurement of the quench distance was only one aspect of the phe- nomenon. Because chemical reactions that occur during combustion are highly temperature dependent, it is desirable to determine what happens as different species are consumed and produced. Although limited exper- imental work has been performed, numerical studies are the primary tool used to study kinetics. It was found by multiple researchers that unburned hydrocarbons near the wall diffuse back into the hot region of the flame and are oxidized quite effectively over timescales of order 1 ms[9, 20, 54]; more will be said on the subject when chemistry modelling is discussed. The extent to which UHC from wall quenching contributes to exhaust emissions is difficult to assess as there are many coupled effects. Amano and Okamoto attempted comparison by taking experimental data from a methane engine and using that as input to numerical simulations that cal- culated the quench distance and assumed all fuel in that zone was unburned. This assumption led to an overestimate of the quenching contribution[4]. A similar study of wall quenching in diesel engines attempted to isolate quench- ing by controlling the wall temperature and combustion phasing, with lim- ited success[38]. At the moment of quenching, the near-wall region experiences a large temperature gradient which results in elevated wall heat flux, as well as 5 2.2. Thermal Considerations rapid changes to temperature dependent quantities such as density, ther- mal diffusivity and pressure. It is these gradients that make wall quenching difficult to model directly, as the time and distance scales are exceedingly small. To illustrate, Figure 2.1(a) shows results of a numerical study that displays the normalized values of quenching parameters for two different wall temperatures; the horizontal time axis is scaled by the flame time tf , here defined as tf = Sl/δf . As stated in the source “[a]ll parameters are normalized using their laminar undisturbed values except the wall heat flux and the Peclet number. The wall heat flux is normalized by the laminar flame power [QΣ] and the Peclet number yields the flame-wall distance nor- malized by the laminar flame thickness, based on the maximum temperature gradient.”. The order of magnitude of QΣ , tf and δl are 0.5 MW/m 2, 1 ms and 0.1 mm respectively. Figure 2.1(b) shows the dimensional wall heat flux and timescales for quenching of a propane flame. 2.2.2 Strained Diffusion Flames There has been limited work performed to extended the study of quenching distance and heat flux under conditions other than homogeneous stoichio- metric mixtures. Delataillade et. al. performed a numerical study of the “head on quenching with strain” (or HOQS) configuration using single step chemistry[15]. A schematic of the configuration is shown in Figure 2.2(a); this is the most easily controlled configuration of diffusion flames. In Figure 2.2(b), the effect on a diffusion flame of increasing strain rates on quenching distance and peak heat flux are shown. The strain rate on the horizontal axis was scaled with tf = α/δq,st , the quench distance on the left hand axis with δl , and the peak heat flux on the right axis with QW = ρSlcp∆T . All quantities were measured at the stagnation point of the flow. As a point of reference, the peak value of Φ for the stoichiometric mixture was 0.13. This result shows that the peak heat flux for strained flames can in fact exceed the stoichiometric premixed heat flux, but that the quench distance does not appear be any smaller than the laminar flame thickness (the de- creasing trend could possibly extend for higher strain rates). This suggests that if a suitable quenching model for quiescent premixed flames can be found (that is, one that accounts for differences gas composition, pressure and temperature), then it might be suitably applied to strained diffusion flames. It also stands to reason that sub-sonic motion of cold walls would have no effect on peak heat flux or quenching distance, as these quantities should be the same if measured in either a stationary reference frame or one moving with the wall; in other words the strain could be imposed by either 6 2.2. Thermal Considerations (a) Time dependence of quenching parameters[43]. Reprinted from Combus- tion and Flame, Volume 108, Issue 3, P. Popp and B. Baum, Analysis of wall heat fluxes, reaction mechanisms, and unburnt hydrocarbons during the head-on quenching of a laminar methane flame, p.333, ©Combustion and Flame (1997), with permission from Elsevier Figure 2.1: Quenching quantities 7 2.2. Thermal Considerations (b) Heat flux magnitudes[16]. This article was published in the Symposium (International) on Combustion, Volume 24, Issue 1, O. Ezekoye, R. Greif and R. F. Sawyer, Increased surface temperature effects on wall heat transfer during unsteady flame quenching, p. 1467, ©The Combustion Institute (1992) with permission from Elsevier Figure 2.1: Quenching quantities (cont.) 8 2.2. Thermal Considerations Wall Cold  Oxidizer Flame Front Fuel (a) Head on quenching with strain (b) Quenching distance and heat flux Figure 2.2: Quenching of strained diffusion flames (adapted from [15]; axis labels have been added by the author). Reprinted from Proceedings of the Combustion Institute, Volume 29 Issue 1, A. De Lataillade, F. Dabireau, B. Cuenot, T. Poinsot, Flame/wall interaction and maximum wall heat fluxes in diffusion burners, pp.796,799, ©2002 The Combustion Institute with permission from Elsevier 9 2.2. Thermal Considerations acceleration of the flow of by moving the walls. 2.2.3 Flame Wall Interaction at a Distance The inhibition of combustion to heat losses is only one effect that cold walls have on flames. Because combustion couples chemistry, thermodynamics and fluid dynamics, it is important to consider what effects walls have on each of these different areas. In numerical simulation, the smallest turbulent length scales are rarely resolved due to computational expense. Instead, models are used that either primarily account for the time averaged effects of turbulence of the flow (Reynolds Averaged Navier Stokes or RANS models), or simulate only the largest turbulent eddies and use sub-grid models for the effects of smaller ones (Large Eddy Simulation or LES). As flow approaches the wall, the size of turbulent eddies decreases and the flow will become essentially laminar at some distance from the wall. This type of interaction becomes important when modelling these effects in numerical simulation. A common way of coupling the effect of chemistry and turbulence is through flamelet models; the local flame structure is as- sumed to be that of a laminar flame. The flame is broken up elements, and these ’flamelets’ are then transported along with the turbulent flow. Un- der these conditions, flames are strained (stretched) which increases their area, and in the case of diffusion flames enhances the mixing of fuel and oxidizer. The laminarization of a flame inhibits both of these processes and reduces the extent of the flame[40]. One approach to address the near-wall laminarization has been to incorporate a sink term in the transport equa- tion for flamelets at cells adjacent to the wall[41]. The strength of this sink term, determined by the time- and distance scales over which flames quench, is related to the quenching observed in laminar flames; most studies have found that there is agreement to within the order of magnitude between the laminar and turbulent peak heat flux and δq . This model was used in a qualitative study under engine-relevant conditions and found to produce results that mirrored physical behaviour and gave better agreement for local quantities[41]; no quantitative studies have been performed with this model, but the results have led to the development of more sophisticated models. 2.2.4 Turbulence and Combustion Coupling An additional effect walls have is their influence on flames aside from di- rect quenching. For non-premixed combustion (which is what occurs in high-pressure direct-injection engines), many studies have looked at altering 10 2.3. Chemistry Considerations reactor geometry to improve the mixing of fuel and oxidizer. In this situa- tion, the fluid dynamics portion of the problem is coupled to the chemistry. More directly, it has also been observed that flames propagating near a cold wall are influenced (deceleration of flame front, reduction of flame surface) even before quenching[8, 41]; the presence of walls affects the structure of turbulence at length scales larger than the quenching distance, but smaller than the characteristic length of the domain. In combustion codes that model turbulence, problems arise when the length scales of turbulence fall below a certain level, usually governed by the CFD mesh spacing. This happens near walls, which will cause unphysi- cal behaviour (for example, eddy-break-up turbulence models predicts that flames will travel faster near walls, which is not observed in experiments[40]). As such, even before the flames have reached the wall they are behaving un- like what is seen in reality. Modifying the flame code to account for near-wall turbulence is substantially more difficult, but such changes to the combus- tion model will improve not only predictions about UHC, but all other com- bustion properties as well. Early experimental results have shown that there is an abrupt increase of turbulent intensity near the wall, which affects wall shear and the fine scale variation of wall heat flux[11]. Figure 2.3(a) shows the time evolution of quenching parameters for turbulent combustion from both experimental (Figure 2.3(a)) and numerical (Figure 2.3(b))studies; it is clear that while turbulence affects the time dependence of quenching pa- rameters, peak values are less affected by turbulent variations. 2.3 Chemistry Considerations Up to this point, the literature review has focused on studies of quench- ing from a thermal perspective, as this allows us to obtain the large scale phenomena of quench distance, wall heat flux and flame propagation. When simulating combustion, the selection of chemistry mechanisms can have a sig- nificant impact on the accuracy of solutions as well as the run-time required. Chemical mechanisms can range from the most simple that only involve one reaction (fuel + oxidizer), to multi-species and reaction mechanisms that at- tempt to replicate the many intermediate steps required for combustion[48]. Comparison between one step- and detailed chemical kinetics appears to support the claim that major properties (flame speed, peak heat flux) can be reasonably predicted with a well posed one step mechanism[50]. With the advances in computer hardware there has been greater flexibility with this choice; depending on the application requirements a simpler chemical 11 2.3. Chemistry Considerations mechanism can be chosen to reduce the overall run-time, or a more complex mechanism can be employed to gain greater insight into emissions or similar phenomena. 2.3.1 Chemical Mechanisms There also exists the question of how strong the coupling is (or should be) between chemistry, fluid- and thermodynamics. As an example, two effects that need to be addressed are thermophoresis (or the Soret effect) and the Dufour effect. Thermophoresis occurs when a mixture of different species is subjected to a temperature gradient and each species undergoes thermodif- fusion differently; the Dufour effect occurs when a thermal gradient is estab- lished through species diffusion. Both of these effects occur over scales much smaller than CFD mesh length scales (as the relevant length scale is the λ, the gas mean free path length) and can be important to the phenomenon of quenching; the intermediate species required for combustion of hydrocar- bon fuels will diffuse from the high temperature region of the flame front to the cooler region near the wall where the high activation energy required for the intermediate reactions is not available. Figure 2.4 shows the species concentration and temperature profiles for a premixed laminar flame using the Trajectory Generated Lower Dimensional Manifold (TGLDM) method, showing both major and minor species, being compared with a detailed chemical mechanism with a total of 71 species and 379 reactions. More will be said on TGLDM later in this work; for now it is sufficient to know it is a method whereby complex chemical mechanisms can be simplified and solved for with far lower computational requirements[23]. From the figure is it evident that each of the species is responding dif- ferently to the exothermic reaction, though this is due to the progress of the reaction; these species gradients could lead to Soret diffusion. In the TGLDM model, differential species diffusion is ignored and the results com- pare favourably with detailed chemistry calculations. The influence of Soret and Dufour effects is typically ignored in most combustion simulation and flame propagation models; a study including these effects determined that the Soret effect could be ignored for quenching on cold walls (300 K < TW < 400 K), but becomes important at higher wall temperatures where radical concentrations increase[43], though a later study has suggested that the “lower” wall temperature could be as high as 600 K[5]. Numerous stud- ies have also shown that the wall temperature rises negligibly (< 10 K) throughout the quenching process[16, 25, 27, 53], so that cold walls can be considered inert throughout quenching. 12 2.3. Chemistry Considerations (a) Turbulent quenching[11]. Reprinted from Proceedings of the Com- bustion Institute, Volume 31, Bastien Boust and Julien Sotton and Marc Bellenoue, Unsteady heat transfer during the turbulent combustion of a lean premixed methane-air flame: Effect of pressure and gas dynamics, p.1414,©2007 The Combustion Institute with permission from Elsevier (b) Turbulent quenching (QW only)[41]. Reprinted from Pro- ceedings of Combustion and Flame, Volume 95 Issues 1-2, T.J. Poinsot and D.C. Haworth and G. Bruneaux, Flame/wall in- teraction and maximum wall heat fluxes in diffusion burners, p. 125, ©1993 The Combustion Institute with permission from Elsevier Figure 2.3: Quenching of turbulent premixed flames 13 2.3. Chemistry Considerations (a) TGLDM Major Species (b) TGLDM Minor Species Figure 2.4: Species profiles through a premixed flame[23]. Reprinted from Combustion Theory and Modelling, Volume 11 Issues 6, J. Huang and W.K. Bushe, Simulation of transient turbulent methane jet ignition and combus- tion under engine-relevant conditions using conditional source-term estima- tion with detailed chemistry, p. 984, ©1993 Combustion Theory and Mod- elling with permission from Taylor and Francis Ltd 14 2.3. Chemistry Considerations 2.3.2 Post-quench Diffusion As mentioned earlier, the initial studies that first suggested wall quenching as a significant source of UHC required revision after it was discovered that the correlations for two wall quenching weren’t applicable to the quenching configuration inside IC engines. It must also be noted that the amount of fuel trapped in the quenching layer does not remain constant after the initial quenching event, and that this can be attributed to species diffusion that occurs after the flame has quenched[3, 9, 43]. Because of the hot burnt mixture that remains near the wall after combustion, the fuel may diffuse out and be oxidized. This can be most easily seen in Figure 2.5 which shows the reduction in fuel concentration near the wall after quenching for two different wall temperatures predicted by Popp and Baum; the time and distance have both been scaled by the tf and δf as defined previously in Figure 2.1(a). Figure 2.5: Post-quench diffusion of boundary fuel[43]. Reprinted from Combustion and Flame, Volume 108, Issue 3, P. Popp and B. Baum, Anal- ysis of wall heat fluxes, reaction mechanisms, and unburnt hydrocarbons during the head-on quenching of a laminar methane flame, p.333, ©1997 Combustion and Flame with permission from Elsevier Popp and Baum’s results suggest that for both wall temperatures, the distribution of UHC within the quenching layer is self similar, and confirms the process is governed by species diffusion. 15 2.4. Modelling Considerations 2.4 Modelling Considerations 2.4.1 Assumptions and Simplifications Based on the examination of the literature and identification of important phenomena, the following assumptions will be made regarding flame wall interaction: • Quenching occurs because the low temperature walls inhibit chemical reactions from proceeding • Wall quenching typically occurs over length scales smaller than the smallest computational cell • Unburned hydrocarbons from quenched flames are reduced due to post-quench diffusion and oxidation • The length and time scales of direct quenching and flame wall inter- action in turbulent flow are of similar magnitude as in laminar flow • Moving walls have negligible effect on quench distance • Wall temperature stays constant during quenching A discussion of the two main types of models used to represent flame quenching will be presented. The basic assumptions, simplifications, advan- tages and disadvantages of each grouping will be summarized to yield an understanding of the rationale behind the final decision to proceed with the chosen model. 2.4.2 Phenomenological - Quench distance calculation One family of models dealing with wall effects looks solely at wall quench- ing. The idea is taken from the initial proposition that the unburned fuel in the quenching layer is responsible for UHC emissions. The size of the quench zone is predicted by modelling the heat transfer during quenching as quasi-steady conduction through a thin slab of gas[13, 40]. The thermal gradient from the hot reactants to the cold wall is the driving force behind this heat transfer, and the heat flux can be directly related to the quench distance. This is a fairly simple way of dealing with near-wall combustion, but the principles are well accepted by many researchers; it also provides a reasonable reflection of what occurs in reality. The disadvantages of this ap- proach include the problem of accounting for the diffusion of species within 16 2.4. Modelling Considerations the quenching zone after quenching occurs, and the calculation of tempera- tures and/or heat flux near the wall with such high gradients. However, so long as the fluid dynamics and heat transfer are adequately represented a reasonable estimate of the quench distance can be calculated. Amano and Okamoto used this type of approach in their analysis of UHC emissions from a pre-chamber compression ignition engine undergoing fully premixed combustion of methane fuel. They experimented with a variety of engine conditions and found that the amount of unburned hydrocarbons from wall quenching was mainly influenced by the global equivalence ratio. Their results are shown in Figure 2.6. Figure 2.6: Unburned hydrocarbons from pre-chamber engine[4]. ORG, HTR and HTR+FR refer to the geometry of the engine: Original, High Top Ring and High Top Ring + Fire Ring respectively. Reprinted with permission from SAE Publication 2001-01-3528 ©2001 SAE International. The model employed by Amano and Okamoto took the average wall tem- perature (453 K), pressure and equivalence ratio from experiments as input to the combustion simulation software PREMIX to calculate Sl , which was 17 2.4. Modelling Considerations then used in a quench distance correlation that assumed a Peb of 8 (this would mean predominantly sidewall quenching at atmospheric temperature and pressure)[4]. Thus this model did not address the spatial distribution of quench distances; this is not as large an issue for spark ignited engines where the majority of quenching occurs nearly simultaneously on all surfaces, but becomes more important for HPDI engines where the flame front’s interac- tion with combustion chamber walls is non-uniform. A more general model proposed by Boust et. al.[10] calculates the quenching distance as a function of peak heat flux and flame power. The concept is to perform a first law control volume analysis on the flame wall interaction. The wall heat flux is non-dimensionalized by QΣ, defined as QΣ = ρuSl∆HYfuel, (2.6) where ρu is the density of the fresh gases, Sl is the laminar flame velocity, ∆H is the heat of combustion for the overall chemical reaction, and Yfuel is the fuel mass fraction in the fresh gases; this is a measure of the heat release from the flame front. Flame power should increase with pressure, as flame speed decreases with pressure at a slower rate than density increases. The thermal model then balances the heat flux from the flame to the wall with the heat conduction through the quench layer. Details of the derivation can be found in Appendix A; the final result is Peb = 1 Φ − 1, . (2.7) Thus the Peclet number (and the quench distance) can be influenced by local gas conditions (through QΣ ) and heat transfer (through QW ). The study in which this model was proposed compared the measured quench distance from direct photography to that calculated via the model and the measured peak heat flux. The model was found to predict the quenching distance for both configurations and over a pressure range of 1 to ∼3 bar. This model has predicted quench distances for head-on- and sidewall quenching under stoichiometric and lean conditions and several initial pres- sures. The results of this model are shown in Figure 2.7; the symbol φ refers to the equivalence ratio of the mixture. The model predicts that even at a pressure of only 2 bar, the quench distance is reduced to a value of roughly 100 µm for a stoichiometric mix- ture. This is substantially less than that predicted by Cleary and Farrell’s polynomial wall temperature correlation (∼ 373 µm)[13]. There is also a strong dependence on equivalence ratio, with δq for an equivalence ratio of 0.7 being higher than the stoichiometric case over all pressures. This result suggests that a thermal model can predict the proper quench distance. 18 2.4. Modelling Considerations Figure 2.7: Thermal quenching model[10]. Reprinted from Combustion and Flame, Volume 149, Issue 3, B. Boust, J. Sotton, S.A. Labuda, M. Bellenoue, A thermal formulation for single-wall quenching of transient laminar flames, p.333, ©2007 Combustion and Flame with permission from Elsevier 2.4.3 Physical - Flame structure modification The other approach to modelling the effects of walls on combustion is mod- ify the flame structure and speed as the flame approaches the cold wall. The results from single wall quenching studies are still used, but instead of directly trying to determine a quenching distance or heat flux, the informa- tion is used to inform the relative sizes of the time and length scales over which walls influence flames. From this information source or sink terms are introduced into the transport equations, with some model parameters to tune their effects. Poinsot et. al. proposed a sink term in the flamelet transport equations [41]. The strength of this sink term is determined by the time- and distance scales over which flames quench. These scales were extracted from 2D DNS studies of turbulent wall quenching, and were found to be closely related to the quenching observed in laminar flames; studies have found that lam- inar and turbulent quenching have similar magnitudes of quench distance 19 2.4. Modelling Considerations and peak heat flux[11]. This model was used in a qualitative study under engine-relevant conditions and found to produce results that mirrored physi- cal behaviour and gave better agreement for local quantities; no quantitative studies were performed with this model; rather it was used as a launching point for future work. A more sophisticated method is to account for the destruction of flame elements close to the wall in a way that is globally applied but only fully active at the boundary cells. Bruneaux, Poinsot and Ferziger used DNS re- sults of turbulent channel combustion and proposed closure for all the terms in the flamelet transport equation while accounting for enthalpy loss at cold walls [12]. This approach has been used in engine simulations successfully, but only for that particular turbulent combustion model[7]. 2.4.4 Model Selection The current combustion model used in Westport’s simulation software em- ploys conditional source term estimation (CSE) with the Trajectory Gener- ated Lower Dimensional Manifold (TGLDM) to predict reaction rates and species concentrations[23]. Most amended models that account for flame- wall interaction are based on flamelet approaches, and as such would require adaptation to the CSE-TGLDM method (to the author’s knowledge there are no such models currently in existence). These models are the most chemically and physically accurate ways of accounting for flame quenching, and in principle are the most desirable. Difficulties arise in the initial development of the model, with the closure of transport equations and the introduction of source terms requiring DNS data of flame wall interaction to formulate closure relations. Though this type of model is attractive from an academic point of view, it would be more suited to an advanced research project. After considering both possibilities, it was decided to use the thermal model of Boust as the basis for the computational quenching model. The reasons behind this decision include the fact that the existing combustion model adequately predicts chemical species in the majority of the domain, and that it is only the amount of UHC that exists after quenching that is desired. 20 Chapter 3 Model Development and Validation 3.1 Methodology To use a detailed chemical mechanism, the Trajectory Generated Lower Di- mensional Manifold (TGLDM) method[42] was implemented with the Open- FOAM computational fluid dynamics (CFD) toolbox, version 1.6[35]. This method has been previously validated with various combustion mechanisms and conditions[23, 51]. TGLDM allows the description of complex chemistry by separating out the shortest timescales that govern chemical kinetics. This allows projecting reaction trajectories computed from detailed chemistry to a low-dimensional table. Retrieving the reaction rate from TGLDM be- tween steps of CFD calculation based on a progress variable Ypg reduces the computational time by nearly two orders of magnitude compared to directly solving a large system of coupled ordinary differential equations (ODEs; this excludes solving for the flow field which accomplished using conven- tional finite volume methods); details of the solver operation can be found in Appendix B. Using this approach yields the ability to run simulations far more quickly than using OpenFOAMs built in ODE solvers; as an example, at 1 bar a case with comparable mesh density takes 3-4 days to perform on a quad- core desktop computer using OpenFOAM built in chemistry solver, whereas TGLDM only requires approximately 6 hours on a single core machine. 3.2 High Pressure Premixed Flames To determine if the Boust thermal model could be adapted for use in engine simulation, it was important to verify that the model could be applied to the high pressure diffusion flames present during HPDI combustion. Breaking down the problem it was decided that high pressure stoichiometric flames would be an important first step for several reasons, which include the fol- 21 3.2. High Pressure Premixed Flames lowing: • To establish a baseline of δq values for high pressure combustion which were unavailable at the time of the study • To evaluate the performance of the TGLDM method for solving chem- ical schemes for predicting flame speed and peak heat flux • To assess the accuracy of the Boust model under controlled conditions before proceeding to high pressure diffusion flames For this purpose, a premixed stoichiometric flame was simulated imping- ing on an rectangular obstacle. The simulation domain is shown in Figure 3.1, and was considered two dimensional (the grid resolution in the figure is 8 times coarser in each direction than that used in simulation, and is shown to give an idea of the refinement used near walls). The stagwall bound- ary measured 3 mm and the overall solution domain measured 19 mm by 9 mm. This setup is partially based on the recent experimental work investi- gating quench distance, where combustion took place in a constant volume chamber[10]. As the flame initiated in that study was spherical, the dis- tance between the wall and the ignition source had to be larger to minimize the effects of flame stretch. Comparison between point and plane ignition sources showed no appreciable effect on parameters such as flame speed and wall heat flux, and thus the plane flame configuration was employed. Figure 3.1: Computational domain Points A and B are where the wall heat flux is probed for stagnation and sidewall quenching respectively. Point B was placed 4 mm away from the in- tersection the stagWall and sideWall boundaries, to ensure that the flame 22 3.2. High Pressure Premixed Flames had cleared the corner and stabilized to a sidewall configuration. The bound- ary conditions for the open sections of the domain were constant pressure and zero gradient for temperature and velocity, except inletHeat (where the flame is initiated by raising the temperature to 2000 K over 0.2 ms) and centerLine which is a symmetry condition. The species mass fraction had zero gradient conditions on all boundaries. Table 3.1(a) shows the experi- mental matrix that describes the temperature and pressure conditions; the stoichiometric gas mixtures were assumed to be composed of fuel and dry air, with the initial species mass fractions shown in Table 3.1(b). Table 3.1: Experimental conditions (a) Temperature and pressure TW [K] Pressure [bar] 300 1, 2, 20, 40 600 1, 2, 20, 40 (b) Stoichiometric mixture fractions Fuel Yfuel YO2 YN2 Methane 0.055 0.22 0.725 Heptane 0.062 0.218 0.72 3.2.1 Convergence The computations were performed on a structured rectangular grid, with the mesh density increasing normal to wall boundaries. The ratio of the first cell near the wall to the last cell of the expansion region was 1.4. To ensure numerical convergence, the temperature profile normal to the wall at the mo- ment of quenching was calculated on three different grids, with the density of cells increasing by a factor of 2 in each direction for every iteration; the total number of cells in each grid were 93.6K, 374K and 15M respectively. Two related quantities were observed for convergence: the temperature profile out to 2 mm from the wall and the peak heat flux for stagnation quenching. Both of these quantities were calculated for a methane-air mixture at 300 K, wall temperature of 300 K, and pressure of 1 bar; this was done so that the results could also be used for validation, since the majority of experimental data in the literature is obtained at lower pressures. Temperature Profile The temperature profile was re-sampled up to 2 mm normal to the wall at the locations marked A and B in Figure 3.1. A comparative plot of the temperature profiles for each of these grids is shown in Figure 3.2. 23 3.2. High Pressure Premixed Flames  0  500  1000  1500  2000  2500  0  0.0005  0.001  0.0015  0.002 T e m p e r a t u r e  [ K ] Distance From Wall [m] 93600 Cells 374400 Cells 1497600 Cells (a) Stagnation T profile  0  500  1000  1500  2000  2500  0  0.0005  0.001  0.0015  0.002 T e m p e r a t u r e  [ K ] Distance From Wall [m] 93600 Cells 374400 Cells 1497600 Cells (b) Sidewall T profile Figure 3.2: Convergence of temperature 24 3.2. High Pressure Premixed Flames It is obvious that the second grid captures all the necessary features of the temperature profile for the stagnation profile, while the sidewall profile is captured well even on the coarsest grid. This is consistent as the sidewall quench distances are typically about twice that of the stagnation case[54, 41, 8]. All the reported results for this report were calculated on the second grid. The order of convergence was calculated with these results was calculated to be 2.88, which indicates roughly third order accuracy; this will be taken as the order when calculating the error in the quenching distance, which is based on the temperature profile. Peak Heat Flux The wall heat flux was calculated using the formula QW = α ∂h ∂n ∣∣∣∣ W , (3.1) where ∂h∂n ∣∣ W is the normal derivative of enthalpy evaluated at the wall and α is the thermal diffusivity. This is slightly different than the usual definition in Equation 2.2 This is due to the particular formulation of the transport equations in the OpenFOAM framework. A convergence calculation similar to that of temperature profile was performed for the peak heat flux value, by calculat- ing the peak heat flux on grids of finer and finer resolution. The estimated numerical uncertainty was calculated to be 2.8%. The convergence for peak heat flux is much slower, with an order of 0.866; as this quantity is calculated at the boundary, we expect a lower order of convergence. The extrapolated relative error was calculated to be 2.8% on the medium grid; all peak heat flux values are taken to have the same relative error. 3.2.2 Flame Speed One of the first measures of simulation validity was the predicted laminar flame speeds. To accomplish this, the displacement of the reaction zone (where Ypg changes from 0 to 1 rapidly) was tracked as a function of time in post-processing. The location of the flame was defined as the area-weighted average of the x coordinate of all cells whose value of Ypg is within the range [0.6,0.8]. Flame speeds were calculated for all time steps and then averaged to obtain a single figure. The uncertainty was taken as the twice the standard deviation of the time average, as this was considered the greatest source of 25 3.2. High Pressure Premixed Flames  1  10  100 100 101 S l [ c m / s ] P [bar] TGLDM Results Gu[18] Rozenchan[41] Warnatz[46] (a) Methane  1  10  100 100 101 S l [ c m / s ] P [bar] TGLDM Results Extrapolated from Smallbone[42] (b) Heptane Figure 3.3: Laminar flame speeds 26 3.2. High Pressure Premixed Flames error. The results for methane and heptane are shown in Figures 3.3(a) and 3.3(b). These values mostly lie within experimental agreement with the cited literature. There are fewer studies on the flame speed for heptane, and even fewer at higher pressures. To provide some comparison, data was extrap- olated from recent experimental data at pressures from 0.5 to 2 bar and temperatures 298 and 300 K[47]. An empirical laminar flame speed correla- tion similar to that of Gülder[19] was used: Sl = S 0 l ( T T0 )a( P P0 )b , (3.2) with the exponents a and b calculated to be 1.54 and -0.43 respectively; S0l , T0 and P0 refer to the laminar flame speed, temperature and pressure obtained under reference conditions (typically 300 K and 1 bar, as is the case for this work). As such, the close agreement at high pressures for heptane should not be considered exact. All subsequent calculations involving Sl used the values computed in this section. 3.2.3 Heat Flux The transient heat flux was calculated as described in the validation section at the locations marked A and B for stagnation and sidewall quenching respectively. The calculated values of heat flux are shown in Tables 3.2, 3.3 for both fuels and quenching configurations. (a) Methane QW [MW/m 2] @ P [bar] TW [K] 1 2 20 40 300 0.315 0.472 1.198 2.210 600 0.298 0.569 1.830 2.609 (b) Heptane QW [MW/m 2] @ P [bar] TW [K] 1 2 20 40 300 0.425 0.667 1.687 1.717 600 0.430 0.720 0.648 0.533 Table 3.2: Peak heat fluxes, head on quenching 27 3.2. High Pressure Premixed Flames (a) Methane QW [MW/m 2] @ P [bar] TW [K] 1 2 20 40 300 0.191 0.309 0.893 1.227 600 0.211 0.343 0.978 1.477 (b) Heptane QW [MW/m 2] @ P [bar] TW [K] 1 2 20 40 300 0.258 0.438 1.084 1.555 600 0.288 0.469 0.353 0.327 Table 3.3: Peak heat fluxes, sidewall quenching The heat flux values calculated are at least 30-40% lower than those recorded by previous researchers for the pressures where data is available[10, 16, 43, 27]; the higher pressure cases show further deviation from the trend. This is important for evaluation of the quench distance model of Boust et. al., as it uses peak heat flux as an input parameter. It should be noted that the heat flux calculation does not account for gas compressibility; a recent paper by Rakopoulos et. al. has shown that neglecting compress- ibility results in erroneous QW [45], though to a degree not as severe as the present study. It is also clear that this trend is not observed for the sidewall quenching of heptane at 600 K wall temperature. In the future a better calculation for heat flux should be implemented; for the present study these figures will be used as they follow the trend of increasing heat flux at higher pressures (non-linearly with QW ∝ P a, a ' 0.5[27]) and higher temperatures for most of the cases; flame wall interac- tion for the heptane flame at pressures of 20 and 40 bar and 600 K wall temperature did not follow this trend. There was insufficient time to fully explore the phenomenon and determine an exact source. Refining the spa- tial grid has little effect on the peak heat flux (changes were limited to less than 10%); it is possible that there are problems with computational time stepping or the generated manifolds. However, the trend of heat flux in- creasing with pressure is similar to that seen in work of Boust et. al., albeit with lower values of the peak heat flux as discussed earlier. 28 3.2. High Pressure Premixed Flames 3.2.4 Quenching Distance To calculate the quenching distance based on the temperature profile, a sim- ilar approach to that of Westbrook et. al. was used, who calculated quench distance for methane and methanol flames under pressures and equivalence ratios in the range of those in this study[54]. It was found that either calcu- lating the location of maximum heat release or the location where the flame temperature was 1500 K gave similar results for methane and methanol; the second method was used in the present study for both methane and heptane due to its simple implementation. Comparison of values from this study with previously published results was complicated by the different definitions of δq and the varied methods by which it is measured. For example, in Westbrook’s study[54] δq was defined by the near-wall temperature profile and rate of heat release, in Boust’s study[10] it is obtained from photographic analysis and in Cleary and Farrell’s[13] it is defined by the concentration of intermediate species; this last measurement argues that δq for the stagnation case is larger than for sidewall quenching which contradicts other studies. The present values are provided as comparison of the trend and order of magnitude. The distances from this study follow the correct trend and are of the right order of mag- nitude. As no studies on single wall quenching exist for heptane flames, no comparison was possible. The results show reduction of quenching distance with increasing pressure and wall temperature. Boust Heat Flux Model Using the values of Sl and QW calculated previously (and α from the CFD simulation), δq for each of the cases was determined, and the results are shown in Figures 3.4, 3.5 and 3.2.4, together with the δq results using the temperature profile when available. The values obtained for quench distance at low pressure for methane compare well with Boust’s results. Cleary and Farrel’s results were not cited as they show the laminar stagnation quench distance to be larger than the sidewall distance (contradicting most other studies), though the order of magnitude is similar. This was surprising as the heat flux values were considered incorrect. It must be stressed that these results are very sensitive to the nondimensionalization of QW which varies from study to study; for example Popp and Baum report a peak wall heat flux of 0.51 MW/m2 and a flame power of 0.94 MW/m2 for methane at 1 bar, leading to a non-dimensional flux of 0.54. This is higher than the Φ = 0.2 from 29 3.2. High Pressure Premixed Flames  0  50  100  150  200  250  300  350  400  450  500 100 101 δ q  [ µ m ] P [bar] T Profile, 300 K T Profile, 600 K Thermal Model 300 K Thermal Model, 600 K Westbrook[54] Boust[10] Figure 3.4: Methane head on quench distance  0  50  100  150  200  250  300  350  400  450  500  550 100 101 δ q  [ µ m ] P [bar] T Profile, 300 K T Profile, 600 K Thermal Model 300 K Thermal Model, 600 K Boust[10] Figure 3.5: Methane sidewall quench distance results 30 3.2. High Pressure Premixed Flames  0  50  100  150  200  250  300 100 101 δ q  [ µ m ] P [bar] Head On T Profile, 300 K Sidewall T Profile, 300 K Head On Thermal Model 300 K Sidewall Thermal Model, 300 K Figure 3.6: Heptane quench distance results 31 3.3. Strained Diffusion Flames the present study. In an ideal implementation, both QW and QΣ would be correctly calculated and give more confidence to the model. The results so far are still very promising, as these issues can certainly be addressed through improvement of the heat flux calculation. 3.3 Strained Diffusion Flames With confidence that the thermal model described by Boust could be used with high pressure premixed flames, the next step in validation was deter- mining how the model behaved for diffusion flames under strain. It was initially thought to investigate each of these effects (diffusion and strain) individually, but the results from DeLataillade suggested that both effects could be looked at in one experiment[15]. The case of head on quenching with strain was investigated in a manner similar to DeLataillade’s study, with a strained diffusion flame being subjected to stagnation point flow. 3.3.1 Initial and Boundary Conditions The initial flame was created by creating a Gaussian temperature profile, centred 0.04 m away from the wall with a peak height of 2400 K and width of 2 mm, to ensure ignition of the mixture; the walls and near-wall gases are initially at a temperature of 300 K. The initial configuration for the velocity field was of an inviscid stagnation flow, with the velocity field u defined as u = a(x̂i− yĵ), (3.3) where a denotes the strain rate; x, y, î and ĵ are the coordinates and unit vectors of the domain respectively. The inlet boundaries were set to time constant, spatially varying values consistent with the initial velocity field; the outlet boundary was given typical outflow conditions (constant pressure and zero velocity gradient). The pressure was initially uniform at 10 bar and was kept constant at the stagnation point. The initial reactant profiles were set using an error function, to be con- sistent with the Gaussian temperature distribution. The mass fraction of fuel was set to transition from the stoichiometric value to zero, with nitro- gen used to make up the remainder of the gas in the pure fuel and oxidizer regions. A schematic of the initial temperature and species mass fraction profiles, as well as the velocity field for a strain rate of a = 1, are shown in Figure 3.7; the maximum velocity magnitude for the selected case is 7 m/s. 32 3.3. Strained Diffusion Flames  0  500  1000  1500  2000  2500 −0.05 −0.04 −0.03 −0.02 −0.01  0  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 T [ K ] Y i x [m] T YN2 YO2 YCH4 (a) Temperature and Major Species (b) Velocity Field Figure 3.7: Initial conditions for strained diffusion flame 33 3.3. Strained Diffusion Flames The previously described initial conditions do not completely define a flame at the initial time, but provide enough information so that the mix- ture is ignited; a consistent flame is established within the first dozen time steps (well before flame wall interaction). This method was tested by only specifying the temperature profile in a stoichiometric mixture and observing the flame wall interaction The mixture was ignited and the flame propagated in a manner similar to the high pressure study done earlier. This would have been a better method of initiating the flame for the high pressure study cases, but the required pre-processor had not been installed with the default con- figuration of OpenFOAM (the utility was only found later). Based on the species mass fraction profiles, the equivalence ratio at the peak temperature will be lean (φ < 1), but the exact value will be determined during the simulation by the quenching model. 3.3.2 Heat Flux Calculation As mentioned earlier in the literature review, the values of δq and QW are not significantly affected by flow turbulence. However, turbulence modelling is necessary to properly simulate the fluid flow. To determine when the flow transitions into turbulence, the Reynolds number based on domain size and velocity U = aX, where X is the domain size was used; it therefore is defined as Re = ρaX2 µ . (3.4) Turbulence modelling was activated for Re ≥ 105, using the k −  tur- bulence model with the standard constants. To calculate the wall heat flux, the typical formula is QW = k ∂T ∂n , (3.5) where k is the thermal conductivity of the fluid and the derivative in n denotes it taken normal to the wall. One problem in calculating this quantity in turbulent CFD simulations is that the near-wall velocity and temperature gradients are unresolved. One way to account for these are law of the wall models which model the flow as having a near-wall laminar region and a “outer” turbulent region. The literature on the development of these wall models is well established, and they have been applied in non-reacting flow simulation. 34 3.3. Strained Diffusion Flames Problems can arise with applying law of the wall models to reacting flows, as they typically assume isobaric and isothermal conditions. The isobaric assumption is not as restrictive in this case, as the process of quenching and heat transfer occurs at a roughly constant pressure. More troubling is the isothermal assumption, since the temperature ratio between the wall and the fully combusted mixture can be as large as 6. The included post-processor for the wall heat flux of reacting flows attempts to circumvent this problem by using an alternate formulation, QW = αt ∂h ∂n , (3.6) where n is the direction normal to the wall and αt is the turbulent thermal diffusivity. When performing this calculation, the built in post-processor evaluates the gradient normal to all cell faces using near-wall cell center distance and boundary values (this gives first order accurate heat flux con- sistent with the results of Section 3.2). When the mesh is non-orthogonal (as is the case during engine simulation), correction factors are applied that are specified by the user at run time to limit influence of non-orthogonality on the gradient. For all simulations in this work, the correction was such that the magnitude of the non-orthogonal portion is limited to be less than or equal to the magnitude of the orthogonal part. αt is calculated as αt = µt/Prt , (3.7) where µt and Prt are the turbulent dynamic viscosity and Prandtl number respectively. These turbulent quantities model enhanced momentum and energy transfer due to turbulent mixing. In standard wall function theory, µt is calculated by the formula µt = { µW y + 1 κ log(Ey +)− 1) : y+ > y+lam µW : y + ≤ y+lam (3.8) where κ and E are model constants; the following relations for y+ and y+lam: y+ = C0.25µ yδ √ kW ρWµW (3.9) y+lam = log(Ey + lam)/κ (3.10) where yδ is the near-wall cell distance. In the current version of OpenFOAM, wall functions can be specified at run-time and the default value of the constant overridden by the user. For 35 3.3. Strained Diffusion Flames Table 3.4: k −  wall function summary Coefficient Value Coefficient Value Cµ 0.09 σ 1.3 C1 1.44 Prt 1 C2 1.92 E 9.2 C3 -0.33 κ 0.41 σk 1 the present study the default wall functions of the k −  turbulence model are used, with standard numerical constants (listed in Table 3.4). From section 3.2, the peak heat flux for laminar premixed stoichiometric flames calculated was at least 30% lower than that reported in the literature. Nevertheless, the trend over pressure and wall temperature was consistent with those reported; as the trend in heat flux was most important, this method was used for the case of strained diffusion flames, in part to evaluate the method for turbulent flow. It should be noted that through discussion with other users on the OpenFOAM message boards that the reactingFoam solver (on which the current solver was based) exhibits temperature incon- sistencies in the laminar limit. 3.3.3 Flame Power Calculation For the current implementation of the model, the Boust definition of QΣ was used which required calculation of Yfuel, ρ, Sl , and ∆H. In a premixed flame Sl and ∆H can be determined a priori ; this is not the case in a strained diffusion flame. It will be necessary to calculate these quantities either directly from the simulation or through simplified models. Heat of Reaction Heat of reaction is in general weakly dependent on pressure and more sen- sitive to temperature. As the flame propagates towards the wall, the tem- perature of gases adjacent to the wall will rise, resulting in a larger value of ∆H being used in the thermal model (up to temperatures of ∼2200 K). In the general case where mixing of fuel and oxidizer will be inhomogeneous throughout the domain, flames in near stoichiometric mixtures will have a higher flame temperature, and the compression and expansion of gases will also change the local temperature. 36 3.3. Strained Diffusion Flames In the paper where the Boust model is proposed, the global heat of reaction is used in the definition of flame power (that is, the amount of energy per unit mass of fuel for the overall reaction)[10]. It was proposed that it would be faster to calculate the heat of reaction over a temperature range sufficient to cover all operating conditions and store it as a table. Values of ∆H could then be interpolated from this table for use in the thermal model, using the same TGLDM method as was used in the main CFD computations. Over a temperature range of [500-3500] K, the change in ∆H was less than 2%, so it may be possible to neglect this variation; for completeness it was retained in this version of the model. Flame Speed The method used in Section 3.2 during post-processing to calculate the speed of planar, stoichiometric premixed flames is more difficult to implement at run-time since the discrete time step is much smaller and the flame front may no longer be a planar surface; it would also likely be computationally intensive and thus increase simulation time. To simplify and expedite flame speed calculation, the semi-empirical Gülder model was used to calculate the laminar flame speed, which has the form Sl = ZWφ n exp[−ξ(φ− 1.075)2] ( T T0 )α( P P0 )β , (3.11) with P0 = 1.013 bar and T0 = 300 K being the reference conditions. φ is the equivalence ratio, a common way to express the amount of fuel in a given mixture relative to a stoichiometric mixture. It is defined as φ = 1 λ , (3.12) where λ is the relative air-fuel ratio. For this particular implementation of the transport and chemistry equations, a conserved scalar approach is used to track the relative amounts of fuel and oxidizer. The variable ζ1 is used to monitor methane fuel. The relation between the conserved scalar ζ1 and the equivalence ratio is φ = ζ1 ζ1 ,st , (3.13) where ζ1,st is the value of ζ1 when the mixture is stoichiometric. The sub- script 1 denotes variables that are related to methane combustion. Though methane is the only fuel used in these simulations, additional fuels will be 37 3.3. Strained Diffusion Flames introduced and will require other indices. The remaining parameters are experimentally determined constants. This model has been well validated against large sets of experimental data[19] and is also included in Open- FOAM. The constants of the Gülder model for methane reacting with air are summarized in Table 3.5. Table 3.5: Gülder laminar flame model Coefficient Value Coefficient Value Z 1 ξ 5.180 W 0.422 α 2.000 η 0.150 β -0.500 From the form of the Gülder model in Equation 3.11 and Table 3.5, it is evident that the laminar flame speed is highly dependent on the amount of fuel present in the mixture; because of the factor exp[−ξ(φ− 1.075)2] in the correlation, and particularly the coefficient ξ. With such a large multiplica- tion factor in the Gaussian term, the width of the peak at half maximum is 0.32; thus a large variation in flame speed as a result of inhomogeneous mixtures is expected. This is not a major concern for the current cases, but will be much more important during engine simulation. The Gülder model was already included with OpenFOAM in such a way that the parameters of the model can be changed at run-time; this simplifies diagnostics and adjustments that may have to be made. Since OpenFOAM has implemented flame speed models as dynamically linked libraries, it is also possible to implement new models in the future without having to re- write the solver. 3.3.4 Simulation Results By inspecting the form of the thermal quenching model, it is noted that at high strain rates (when the peak heat flux begins to decrease), Peb will begin to rise again, which can lead to an increased δq (depending on the current local value of δl ). In Delataillade’s study, it appears that δq levels off at a value of δl with increasing strain rate, though it is difficult know with certainty as there was very little data at higher strain rates. Thus the trend of the δq with strain rate from the thermal model would be an important issue, as it may limit its applicability. Figure 3.8 shows peak heat flux and quenching distance results in a manner similar to DeLataillade; here the values are fully dimensional, and 38 3.3. Strained Diffusion Flames are calculated for a uniform mesh spacing of ∆x = 1mm; this is a typical value for simulations and convergence tests found it to have a similar order of convergence as with the premixed quiescent flames. For the prescribed pressure the critical strain rate ac = t −1 f was calculated as 7586 s −1; the results were plotted against the normalized strain rate atf (3.8).  0  20  40  60  80  100  120  140  160  180  200  0  0.2  0.4  0.6  0.8  1  1.2 0.0e+00 1.0e+05 2.0e+05 3.0e+05 4.0e+05 5.0e+05 6.0e+05 7.0e+05 P e b Q W  [ W / m 2 ] a tf QW,max Peb Figure 3.8: Quenching model results for strained diffusion flames Peak heat flux reaches a maximum at a strain rate near atf = 1, similar to DeLataillade’s study. Strain rates much higher than atf = 1.1 caused instability in the solver, and were not considered due to the lack of available literature and time. It is interesting to note that even though the mixture near the flame is lean (φ ' 0.5 − 0.8), the maximum QW still occurs at atf for the stoichiometric mixture. Peb for the maximum heat flux was calculated at 6.6, and δq followed a similar trend with a minimum value of 0.9 mm. For a stoichiometric mixture this would be identified as more like sidewall quenching, but since the mixture is lean overall δq is expected to be larger. For the current study the form of the Boust model was retained, as even for the highest peak heat flux the Peclet number will still be non-zero, positive and give an indication the value of δq . A possible amendment to the Boust model to account for strained flames is discussed in Appendix 39 3.3. Strained Diffusion Flames D. Boust suggests that the flame stretch that occurs in sidewall quenching may be secondary to the thermal losses[10], but more experiments would be needed to determine the effects of strain. A flowchart for the basic form of the computational version of the thermal model is shown in Figure 3.9; a two-zone mass exchange model is included to account for the diffusion and oxidation of fuel from the quench layer. The development and implementa- tion of the model will be discussed in the following chapter. S l ,Q Σ CFD Simulation Obtain (P,T w ,α, etc) From boundary cells Calculate S l  and Q Σ Calculate Q W  and     Use φ to calculate Pe and δ q Apply two­zone mass exchange Adjust CH4 in boundary cells Φ Figure 3.9: Quench distance model 40 Chapter 4 Model in Engine Simulation 4.1 Overview of Algorithm To implement the Boust thermal model for computation in OpenFOAM, it is broken down into the algorithm shown in Program 4.1 Program 4.1 Algorithm for computing quench distance for cells ∈ wallCells do qW = hFlux() {Calculate the heat flux at wall cells} if heat flux has peaked then SigQ[cell] = flamePower(T[cell],Yfuel[cell],p[cell],rho[cell]) {Calculate the flame power based on local temperature, fuel amount, pressure and density} Pe[cell] = SigQ[cell]/qW[cell]-1 {Calculate local burning Peclet num- ber} delQ[cell] = Pe[cell]*alpha[cell]/Sl[cell] {Calculate current quenching distance} Ycor = fuelCorrect(delQ,Yfuel[cell]) {Find how much fuel is left in the quench layer, accounting for diffusion} Yfuel[cell] = Yfuel[cell] - Ycor[cell] {Apply the correction} end if end for 4.2 Validation under Engine Conditions Because the quenching model would be employed in simulation under a va- riety of engine conditions, it was decided to do an initial evaluation with a high pressure case where there was quite a bit of impingement early on in the cycle so that the spatial and temporal evolution of quench distance δq over the piston surface could be analyzed. The test cycle used was the Euro- pean Stationary Cycle (ESC) and is reproduced in Figure 4.1. It shows the operating conditions that Westport’s HPDI engines are tested against. The 41 4.2. Validation under Engine Conditions percentages next to each mode number indicate the percentage weighting factor for each mode. 0 25 50 75 100 10 50 75 En gi ne  L oa d (% ) Engine Speed (%) 1 15 % 2 8 % 3 10 % 4 10 % 5 5 % 6 5 % 7 5 % 8 9 % 9 10 % 10 8 % 11 5 % 12 5 % 13 5 % A B C Figure 4.1: Engine condition map, adapted from[24] A high load, mid-range speed condition (Mode 8) was chosen for the initial validation. The injection was performed with a common rail system where the injection angle of the main natural gas jet was 18°. Table 4.1 shows the relevant case statistics for this mode; note that the amount of injected fuel is for the simulation domain, which is one-third the amount for a full cylinder. The fuel injection profile was a simple square wave pulse whose length was such that the proper amount of fuel was injected; the start of injection is written in terms of degrees after top dead centre (ATDC). The engines studied in this work use a form of pilot ignition, with a small amount of diesel fuel being injected into the cylinder that begins to burn at the desired conditions. Once this flame has been established, the main injection of natural gas occurs which then burns and provides the majority of the power during the expansion stroke. In this case, the pilot fuel is heptane which auto-ignites at a lower temperature and pressure than methane does. Since the amount of heptane injected is so small and its combustion so localized, it is assumed that there is no contribution to wall 42 4.2. Validation under Engine Conditions Table 4.1: Validation mode statistics Parameter [Units] Value Bore x Stroke [mm x mm] 137.2 × 169 Engine Speed [RPM] 1486 EGR [n/a] 0 Injection Pressure [bar] 250 Injected Fuel (CH4) [mg/stroke] 80 Swirl [n/a] 0 Compression Ratio [n/a] 17 Injection Timing, diesel [°ATDC] -13.2 Injection Timing, natural gas [°ATDC] -0.5 UHC from its burning. The simulation process required periodic remeshing of the domain to avoid problems with high aspect ratio cells, and it was found that the majority of the impingement of the flame occurred during the interval of [0,60] CA°ATDC. For later comparison, various quantities related to the quenching distance were tracked at the locations specified in Figure 4.2. Figure 4.2: Location of Peclet number sampling Finally, the turbulence model for engine simulation was changed from standard k−  to the RNG k− , as it was found that it led to better solver 43 4.3. Natural Gas HPDI stability. The model constants are show in Table 4.2. Table 4.2: RNG k −  wall function summary Coefficient Value Coefficient Value Cµ 0.0845 α 1.39 C1 1.42 η0 4.38 C2 1.68 β 0.012 C3 -0.33 E 9.5 αh 1 κ 0.41 αk 1.39 Prt 0.85 4.2.1 Order of Accuracy In the previous sections the simulation domain was regular and structured, from which the order of accuracy and convergence of the numerical schemes could be directly determined by grid refinement. For this and future en- gine cases, the mesh is parametrically generated by a simple C++ program based on general geometric and operating parameters; as such a grid refine- ment study would be project unto itself. The nominal accuracy of Open- FOAM solvers is listed in the User Guide (available on the OpenFOAM webpage[35]), and the results are quoted here. The nominal accuracy of the solvers chosen was second order in space (with upwind schemes and lim- iters employed for specific variables as elaborated upon in Appendix B) and first-order, bounded implicit in time respectively. Because of the unsteady nature of combustion simulation and large vari- ations (in both space and time) of the solution variables, the time step selected played an important role in the stability of the solution. Through analysis of cases solved using the same solver for previous internal stud- ies at Westport, it was found that a Courant number of 0.25 provided an acceptable trade-off between stability and solver speed. 4.3 Natural Gas HPDI 4.3.1 Compression Ignition and HPDI Before discussing the results of the proposed engine conditions, a brief overview of HPDI combustion will be presented to provide context, with additional information regarding the HPDI system of Wesport’s engines; 44 4.3. Natural Gas HPDI much of the information is taken from Heywood’s Internal Combustion En- gine Fundamentals[22]. Compression ignition (CI) engines are known for their high efficiency due higher compression ratios. CI engines compress air in the combustion chamber to high pressure, then introduce fuel near the top dead center (TDC) position of the piston, which then begins to ignite from the heat of the compressed gases. The speed and efficiency at which these ignition and combustion processes occur depends on the type of fuel, the method by which it is introduced into the combustion chamber, the manner in which it is mixed with the air and the condition of the air at the TDC position (pressure, temperature, etc.). Many methods have been used to control the combustion that occurs in CI engines; at present only the method of fuel injection will be discussed. Direct injection has been used since the early days of CI engines, with liquid fuel injected through a nozzle with small orifices such that it is vapourized and begins to burn as it is introduced into the combustion chamber. With direct injection engines, mixing of the fuel vapour and compressed air is critical, as the inhomogeneity of the mixture may prevent complete com- bustion due to regions being too rich or lean. Indirect injection methods have an auxiliary chamber designed so that fuel-air mixing is enhanced and combustion “spills out” into the main chamber. As injector technology has improved, high pressure injection has been seen a means by which to en- hance combustion in CI engines, as it generally leads to better mixing and combustion overall. 4.3.2 Westport HPDI Westport’s HPDI engines use a system for that may best be described as “dual common rail”. The ability to use natural gas fuel (compose primarily of methane) with a compression ignition engine was seen as an opportu- nity to exploit the abundance and renewability of methane fuels, as well as the reduction in undesirable emissions such as carbon dioxide and nitrogen oxides. The difficulty in using methane in a compression ignition engine is the lower higher autoignition temperature of gaseous methane (approxi- mately 1400 K) compared to that of typical diesel fuel (1200 K). This would require substantially higher compression ratios (and much stronger engine components) which would be prohibitively expensive. To address this problem, Westport has developed a system where a small quantity of diesel fuel is injected which begins burn and forms high temper- ature flame kernel. Shortly thereafter natural gas is injected which is then ignited by the high temperature region and combustion proceeds through- 45 4.4. Flame Power Variation Figure 4.3: Westport fuel system diagram out the remainder of the chamber. This system thus achieves the benefits of methane combustion for high compression ratios. Additional informa- tion regarding the development of Westport’s engines has been provided through references in the bibliography ([21, 34] in particular give an excel- lent background in the Westport HPDI system and its advantages, while [37, 36] provide an overview of the effects of injection parameter change on emissions from this system). Figure 4.3 shows a schematic of the Westport components that are added to an OEM engine; both are available through the Westport webpage. 4.4 Flame Power Variation As previously defined, the flame-power non-dimensionalization wall heat flux and is dependent on local mixture properties. Density and fuel mass frac- tion are easily retrieved from the CFD simulation, but flame speed and heat of reaction were not readily available. To address flame speed, the semi- 46 4.4. Flame Power Variation empirical Gülder model (described in the previous chapter) was used, since it was considered valid over the range of temperatures, pressures and mix- ture fractions encountered during this particular implementation of HPDI combustion. Simulations results from strained diffusion flames using a preliminary nu- merical implementation of the thermal quenching model showed promising results, but it must be remembered that ρ, Sl , ∆H and φ will certainly vary greatly during the complex flow and combustion present in engine simula- tion. As such it was considered valuable to examine the quantities used to calculate QΣ throughout the simulation at the locations previously defined. Figure 4.4 shows density plots of ρ, Sl , φ (here denoted as ’eqr’) and QΣ over the surface of the piston, taken 10◦ after first contact of the flame with the piston. A white contour is drawn where φ = 1. It is clear the highest values of QΣ are within a thin region near the stoichiometric surface, and that this region is also where ρ transitions from low to high values. To have a better idea of how sensitive QΣ is to local conditions, the spatial variation the variation of ρ and Sl are examined both as a function of time. 4.4.1 Density Figure 4.5 shows the variation of density with time at the four locations being probed. Though there is a gentle decrease in density due to the expansion of the combustion chamber, it is clear that the density rapidly decreases as the flame front passes through the cell; this is easily explained by treating the gas as ideal and increasing the temperature while maintaining constant pressure (which, according to the literature is a reasonable assumption for cold wall quenching). Ideally, the density used for calculation of the flame power would be that just prior to the sharp increase in temperature. Unfortunately detecting that increase at run time is made complicated by having to define a criterion for what constitutes a sharp increase before the increase is even encountered; this criterion may also change from location to location and over time. For the moment, the value of ρ was simply taken as returned from the simulation; discussion on a possible criterion is located in Appendix D. 4.4.2 Flame Speed Looking at the form of the Gülder correlation introduced in Section 3.3.3, the flame speed should be most sensitive to changes in equivalence ratio, followed 47 4.4. Flame Power Variation by temperature and lastly by pressure (restricting ourselves to these three quantities for our analysis). Since the dependence on pressure is only to the power −0.5, it was decided to have a more detailed look at the variation of flame speed with φ and T . Figure 4.6 shows the results of the flame speed calculated from the Gülder model for a crank angles 10°ATDC, being mass averaged and plotted as contours according to the local temperature and equivalence ratio; the numbered contours in the background represent the amount of UHC produced from a reaction at the specified φ and T after 1 ms. From the figure, it appears the flame speed is concentrated around φ = 1 and two temperatures. The first grouping is near the autoignition temper- ature of ∼ 1400 K; the second is near the adiabatic flame temperature of ∼ 2200 K. Both groupings are consistent with the discussion in Section 3.3.3 regarding the functional form of the Gülder model and its dependence on temperature. As the temperature of wall surfaces for simulation was kept constant, the variation of temperature was considered less important than the mixing and was investigated less thoroughly in this study. While post-processing the injection phase, the stoichiometric surface was visualized to give an indication of the location of the jet within the domain. It was further noted that the peak heat flux would occur on a given cell as the stoichiometric surface approached. This led to speculation that the rele- vant flame speed for subsequent calculations was the stoichiometric laminar flame speed Sl ,st . As the flame speed from the Gülder model is used for nondimensionalization (both through QΣ and δl , it was decided to compare two different ways of calculating QΣ : • Calculate QΣ from local variables at cell centre • Calculate QΣ from stoichiometric flame speed and stoichiometric fuel mass fraction Local φ, T and p The values of QΣ and QW at the probed locations are plotted together in Figure 4.7. During quenching, QΣ peaks very near the moment of maximum QW . Investigating further it was determined that this is due to the dependence of the laminar flame speed on the equivalence ratio. This is consistent with the φ− T plot from Section 4.4.2, where the flame speed is clustered about φ = 1. Although ρ and ζ1 decrease as time moves forward, the effect isn’t as dramatic as the flame speed increase. For this reason it is important to 48 4.4. Flame Power Variation (a) ρ (b) Sl Figure 4.4: QΣ related contours 49 4.4. Flame Power Variation (c) φ (d) QΣ Figure 4.4: QΣ related contours (cont.) 50 4.4. Flame Power Variation  10  20  30  40  50  60  70  0  5  10  15  20  25  30 ρ [ k g / m 3 ] CA [° ATDC] L0 L1 L2 L3 Figure 4.5: Density variation at quenching 51 4.4. Flame Power Variation 0 1 2 3 4 5 6 1000 1500 2000 2500 3000 φ T[K] 1 1 1 1 10 10 10 50 50 100 100 200 200 500 500 900 900 Figure 4.6: φ− T analysis of Sl 52 4.4. Flame Power Variation  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08  1.4e+08  0  5  10  15  20  25  30 Q  [ W / m 2 ] CA [°] ATDC QW QΣ (a) Location 0  0  1e+07  2e+07  3e+07  4e+07  5e+07  6e+07  7e+07  8e+07  9e+07  1e+08  0  5  10  15  20  25  30 Q  [ W / m 2 ] CA [°] ATDC QW QΣ (b) Location 1 Figure 4.7: Local flame power and heat flux comparison 53 4.4. Flame Power Variation  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08  1.4e+08  0  5  10  15  20  25  30 Q  [ W / m 2 ] CA [°] ATDC QW QΣ (c) Location 2  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08  0  5  10  15  20  25  30 Q  [ W / m 2 ] CA [°] ATDC QW QΣ (d) Location 3 Figure 4.7: Local flame power and heat flux comparison (cont.) 54 4.4. Flame Power Variation track the maximum QΣ encountered and store it for use in the quenching calculations. Figure 4.8 shows Peb at the different locations as a function of crank angle. −2  0  2  4  6  8  10  12  14  0  5  10  15  20  25  30 P e b CA [° ATDC] L0 L1 L2 L3 Figure 4.8: Peb from local variables The initial value of Peb = −1 is because QΣ = 0 until the arrival of fuel and the flame front. One difficulty with implementing the thermal model for diffusion flames is that because flame power only peaks very near the moment of quenching, there is a limited window where Peb represents a physically meaningful quenching distance. Thus there exists the need for it to be monitored over time while running the simulation. Though not a problem in post-processing, if this tracking is implemented incorrectly at run-time it could result in erroneous results being carried forward for several iterations. φ = 1, Local T and p Since Sl peaks near stoichiometry, QΣ based on local temperature and pres- sure but for a stoichiometric mixture was tested. The results are shown in Figure 4.9. 55 4.4. Flame Power Variation  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08 −30 −20 −10  0  10  20  30 Q  [ W / m 2 ] CA [° ATDC] QW QΣ (a) Location 0  0  1e+07  2e+07  3e+07  4e+07  5e+07  6e+07  7e+07  8e+07  9e+07 −30 −20 −10  0  10  20  30 Q  [ W / m 2 ] CA [° ATDC] QW QΣ (b) Location 1 Figure 4.9: Stoichiometric Flame Power and Heat Flux Comparison 56 4.4. Flame Power Variation  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08 −30 −20 −10  0  10  20  30 Q  [ W / m 2 ] CA [° ATDC] QW QΣ (c) Location 2  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08 −30 −20 −10  0  10  20  30 Q  [ W / m 2 ] CA [° ATDC] QW QΣ (d) Location 3 Figure 4.9: Stoichiometric flame power and heat flux comparison (cont.) 57 4.4. Flame Power Variation Here the flame power still fluctuates with the temperature and pressure, but is still quite large through the main injection. Examining Location 0, the flame speed peaks with the temperature and the arrival of the flame; the variation with pressure is less dramatic. The flame power is defined at all times and as such the Peclet number can be tracked from even before the time of quenching; of interest is that fact that the maximum stoichiometric flame power is slightly less than the local flame power; this is consistent with the constant φ0 from the Gülder model that represents a peak flame speed for a slightly rich mixture. Because the heat flux monotonically increases during quenching, it can be used as a way to determine whether the quenching calculation is actually valid. This can be further explained by looking at the calculated Peclet number for these locations as a function of time (Figure 4.10).  0  5  10  15  20  25  30  35  40  45 −30 −20 −10  0  10  20  30 P e b CA [°] L0 L1 L2 L3 Figure 4.10: Peb based on stoichiometric flame speed At all locations, the Peclet number reaches a minimum at the moment of peak heat flux (and hence quenching). The magnitude is quite small (in some cases of head on quenching ∼ 1), but the relative magnitudes of the Peclet number between locations of different quenching configurations are consistent. According to this definition, quenching is represented by the 58 4.5. Peb and δl Calculation and Filtering minimum Peclet number. This interpretation is suitable for post-processing of already completed cases, since the discretely recorded time steps may not coincide with the exact moment of quenching. This definition would tend to over predict the flame power (as well as δq) if done at run-time at many locations, as a mixture might not actually reach stoichiometry there. It is important to note that for both methods of calculating QΣ as previ- ously described may continue to fluctuate after quenching due to the trans- port and consumption of fuel. If not accounted for, this change will cause the Peb (and hence δq) to fluctuate in a way that does not represent flame quenching. This was evident in later portions of the expansion stroke; the flame was already quenched but Peb continued to rise due to the change in equivalence ratio as a result of continuous diffusion. In the end, the most effective compromise was to use the cell centre values and filter the results according to a selection procedure described in the next section. 4.5 Peb and δl Calculation and Filtering Because quenching is a transient phenomena, care must be taken to appro- priately account for the changes in the quantities QW , QΣ , P eb , δl and δq . As seen in the previous sections, heat flux rapidly increases when the flame front arrives at the wall, then subsequently subsides as the large tempera- ture gradient disappears. As the timescale of the quenching (0.5 ms) is much larger than that of the computational simulation (1 µs), it was assumed that there was little advantage to any time based interpolation. It was decided to use wall heat flux to determine when quenching had occurred, though nearly all quantities made rapid transitions from one value to another during quenching. The advantage of using QW is that once the flame has quenched, the recorded maximum of QW is unlikely to be encountered again as that would require the flame to retreat and then quench a second time. As QW rises, all other relevant quantities (QΣ , P eb , δl , δq) are also recalculated and compared to the last recorded maximum (QΣ ) or minimum (δl , δq). One exception is Peb , which is always calculated via Equation 2.7; so long as QW and QΣ are at their maximum values, Peb will always be the minimum value. To illustrate this filtering process, Figure 4.11 shows QW , QΣ and Peb at L0 for crank angles after which the piston reaches top dead center (TDC), with the time varying trace and tracked maxima/minima being shown on the same plot. The filtering process for δl requires more explanation. Initially, the lam- 59 4.5. Peb and δl Calculation and Filtering  0  2e+07  4e+07  6e+07  8e+07  1e+08  1.2e+08  1.4e+08  0  5  10  15  20  25  30 Q  [ W / m 2 ] CA [° ATDC] QW QW,max QΣ QΣ,max (a) QW and QΣ −1 −0.5  0  0.5  1  1.5  2  0  5  10  15  20  25  30 P e b CA [° ATDC] L0, Raw L0, Filtered (b) Peb Figure 4.11: Peb filtering 60 4.6. Diffusion/Oxidation Model inar flame thickness is ' 0 on all the wall boundaries since there is very little fuel present for most initial conditions. Wall heat flux begins to rise as the flame approaches and δl based on local quantities is calculated, and then its value is determined as the minimum of • the current δl , • the smallest δl encountered so far and • the shortest cell centre distance to the wall (∆W ). This procedure ensures that the resulting δq represents the smallest value encountered and is always smaller than the local grid spacing. Additionally, these quantities are only updated when the heat flux is increasing, as δq is a quantity that is only defined for flame quenching, and should not be calculated before the flame arrives or after it has quenched. Relative Magnitude of Peb and δq The minimum Peb for premixed quiescent flames at low pressures should be about 2-3, but beyond that no assumptions can be made concerning what the minimum Peb should be for any given patch (as that would require knowing the geometry of flame/wall interaction, which is the point of im- plementing the model). In Figure 4.12 the calculated Peb is shown over the range of [10,30] CA◦ ATDC using the run-time formulation of QΣ and the proposed filtering, at the locations indicated in Figure 4.2; these are plotted all together for easier comparison of their relative magnitudes. The value of Peb helps describe the differences in quenching between the different locations. L0 and L2 are near the initial locations of impingement of flame jets on the piston surface, and correspond to the lowest values of Peb , whereas L3 is a location where the quenching configuration is closer to sidewall, and Peb is correspondingly higher. L1 is located deep in the piston bowl and would be expected to have a Peb much higher. The Peclet number itself is useful for comparison of quenching configurations, but is only an intermediate step to obtaining the quenching distance. 4.6 Diffusion/Oxidation Model Once the thickness of the quenching layer is calculated, the initial amount of fuel left in the quenched region can be estimated; the volume of the quench 61 4.6. Diffusion/Oxidation Model −2  0  2  4  6  8  10  12  14  10  15  20  25  30 P e b CA [° ATDC] L0 L1 L2 L3 Figure 4.12: Variation of Peclet number with location region on a given face is determined by VUHC = δq∆A, (4.1) where ∆A is the area of the wall face of the cell adjacent to the wall bound- ary. The amount of fuel in this region is then calculated as mUHC ,B = ρBYUHC ,BVUHC (4.2) where YUHC ,B is the mass fraction of fuel at the boundary face. From the literature review, it is known that an amount of unburned fuel is subsequently oxidized, resulting in a reduction of the total contribution to UHC emissions. To simulate this process, a simple two zone model was used with the initial amount of unburned fuel being calculated via the mass fraction. The literature shows that the timescale for the oxidization of fuel in the quenching layer is 1 ms or less at atmospheric pressure; with higher pressure this is likely to take longer due to changes in mass diffusion. Since the computational time step is never larger than 1E-5 seconds, it was decided that a simple implementation of Fick’s second law of diffusion would be 62 4.6. Diffusion/Oxidation Model sufficient to model the process; stated for this situation is ∂YUHC ,B ∂t = D ∂2YUHC ,B ∂x2 . (4.3) where YUHC ,B being the amount of fuel at the boundary and D is the mass diffusion coefficient. This equation can be discretized by making the follow- ing substitutions: • ∂t ' ∆t, the CFD time step • ∂x2 ' 2δ2q • ∂YUHC ,B ' YUHC ,cell−YUHC ,B , where YUHC ,cell is the cell centre value These substitutions assume that the largest gradient in concentration occurs at the interface between the quench layer and the bulk region. A schematic of a near-wall cell is shown in Figure 4.13. The second derivative is assumed to be calculated at the center of the quench region, assuming that the largest gradients occur over the length scale of the quenching layer. Strictly speaking this estimate does not rep- resent a rigorous discretization with known error and convergence charac- teristics; rather it represents an estimate of the order of magnitude of the diffusion rate. By making the aforementioned substitutions and rearranging Equation 4.3, the change in the fuel concentration at the wall is obtained: ∆YUHC ,B = −2D δ2q (YUHC ,cell − YUHC ,B )∆t. (4.4) This amount is added to YUHC ,B at every time step, after first checking that the change does not violate conservation of mass. The key parameter in this calculation is the mass diffusivity D. Since nitrogen is in abundance, D will be chosen to represent fuel diffusing into nitrogen. D is related to the viscosity via the Schmidt number Sc, similar to the way that Pr relates the thermal diffusivity and viscosity, and is defined here as: D = µ Scρ (4.5) In all simulations, Sc = 1 since the global CFD simulation assumes neg- ligible differential diffusion due to concentration gradients. Though this may seem contradictory as the model is meant to describe species diffusion, it was deemed necessary to obtain consistent transport properties. Initially, a mass exchange model was considered that would allow fuel in the quenching layer 63 4.6. Diffusion/Oxidation Model q Y CH4 ,ρc  x Y CH4 ,B ,B Physical Field Numerical SolutionY Figure 4.13: Schematic of post-quench diffusion model to diffuse out into the free stream to be oxidized and for the total amount of fuel in the cell to be reduced and subsequently affect the global reaction rate. To see the significant difference, the mass fraction of fuel trapped in 64 4.6. Diffusion/Oxidation Model the quench layer (assuming distribution of fuel is uniform) was calculated and compared with the mass fraction at the center of the computational cell. The results of the diffusion model are shown in Figure 4.14; YCH4 and YCH4,B are traced and compared for each location on the piston.  0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0  5  10  15  20  25  30 Y i CA [° ATDC] YCH4, L0 YCH4,B, L0 YCH4, L1 YCH4,B, L1 YCH4, L2 YCH4,B, L2 YCH4, L3 YCH4,B, L3 Figure 4.14: Comparison of cell center and quench layer values of YCH4 The amount of residual fuel left in the quenching region is reduced at a slower rate as compared to the burning region. Though this part of the quenching model was not evaluated rigorously, the timescale of the reduc- tion is consistent with published literature; with an engine speed of ∼1500 RPM, each crank angle degree is equivalent to approximately 0.1 ms. What is interesting to see is that the amount of reduction of the wall hydrocar- bons is much slower than that in the cell region. The irregular changes in YCH4 at Locations 0, 2 and 3 are due to local changes in equivalence ratio, progress variable and cell temperature. When the TGLDM method reports the species mass fractions, it uses local conditions to interpolate from the previously generated manifolds and determine if the mixture is reacting; when the required conditions are met the field is updated. The resulting mass fractions are not fed back into the CFD calculations and have no im- pact on the overall calculation. 65 4.6. Diffusion/Oxidation Model The reason the flame slows down before encountering a cold wall is be- cause colder temperatures near the wall inhibit intermediate chemical re- actions from proceeding. The unburned fuel near the boundary is a conse- quence of incomplete reactions, and as such there is no feedback mechanism for the unburned fuel to affect reaction in the bulk region. As a result it was deemed unnecessary to have the wall quenching diffusion model change the amount of fuel in the near-wall cells. Finally, the diffusion model was only activated if the cell centre temperature was greater than the combustion cut- off temperature (1400 K), as the unburned fuel in the quench layer would only be oxidized at sufficiently high temperatures. The final computer code for the wall quenching and diffusion models is reproduced in Appendix C. 66 Chapter 5 Engine Cases 5.1 Engine Conditions To test the thermal model’s predictions and estimate how piston quenching contributes to unburned hydrocarbons, several test cases were simulated using real engine operating conditions. A high and low load condition were chosen from the European Stationary Cycle Emission test[24] as shown in Figure 4.1, with the constant engine parameters listed in Table 5.1; Pinj is the injection pressure of the main natural gas jet. The engine timings tested are listed in Table 5.1, where τa and τr represent the advanced and retarded injection timings of both the heptane and natural gas jets; the amount of advance/retardation is 5°from the baseline for all such cases. These conditions represent high and low engine load conditions at a moderate speed. The injection profile for both diesel and natural gas jets was a square wave pulse whose length was adjusted (based on the injection pressure) to inject the correct amount of fuel. The amount of fuel injected is only for the sector that is being simulated and represents one-ninth of the amount for a complete cylinder. 5.2 Geometry and Mesh For all engine simulations, the same piston geometry as the initial validation in Chapter 4 was used; the effect of piston geometry would be studied by examining how flames interact with different portions of the piston surface. The simulation domain was an azimuthal sector of the combustion cham- ber, assuming it to have an axially symmetric arrangement of nine diesel and nine natural gas jets with zero interlace angle (the angle between the centre lines of the methane and heptane gas jets). This configuration reduces the simulation domain from a wedge of 120°(necessary for an injector with 6 diesel and 9 methane jets) to one of 40°; a sample mesh for a crank angle of 60°after top dead centre (ATDC) is shown in Figure 5.1(a). The cold wall boundaries for this mesh are shown in Figure 5.1(b);the temperature on all 67 5.2. Geometry and Mesh Table 5.1: Engine cases (a) Load conditions Parameter [measure] Mode 8 Mode 9 Speed [RPM] 1500 1500 Swirl 1.8 1.8 EGR [%] 19 24 Engine Load [%] 100 25 Fuel Mass Injected [mg] 26.5 5.5 Baseline Timing, diesel [°ATDC] -11.4 -11.4 Baseline Timing, natural gas [°ATDC] -5 -2 (b) Injection conditions Pinj [bar] High Load Low Load 150 N/A τ 250 τ, τa , τr τ 600 τ, τr N/A said boundaries was kept constant at 500 K (this is an estimate based on coolant temperatures and not representative of the real temperature distri- bution over the combustion chamber walls). The region of the cylinder that is compressed and expanded results in significant changes in the cell aspect ratio, requiring periodic remeshing of the domain. The cycle was divided into 5 phases summarized in Table 5.2. Remeshing led to the potential side effect of having quenching quantities (QW , QΣ , etc.) fluctuate when mapping data from one phase to the next. After examining various operating conditions it was found that this was most likely to occur if QW was rising at the wall patch when the data was remapped. The worst case fluctuations were 20-30%, but this was rare and fluctuations were typically less than 10%. Research is being done to replace the remeshing procedure with the addition-subtraction of cell layers throughout the compression-expansion so that the aspect ratio can be better controlled, but it had not yet been implemented as of the publication of this work. 68 5.3. Engine Simulation Results (a) Mesh density (b) Wall boundaries Figure 5.1: Computational domain Table 5.2: Phase breakdown Phase Start CA [°] End CA [°] Cells Notes 1 270 290 41824 Compression 2 290 334 30656 Compression 3 334 390 32768 Injection & Ignition 4 390 420 47728 Expansion 5 420 500 57376 Expansion 5.3 Engine Simulation Results 5.3.1 Comparison with Experiments To show the validity of the numerical simulations, the results of numerical simulation are shown against experimental data collected previously for a similar engine, for the standard baseline injection pressure and timing cases. All experimental data was obtained from earlier testing of a Westport GX engine performed in the test cell facilities at Wesport Innovations, where bulk averages of the in cylinder pressure and heat release rate, as well as emissions and fuel consumption data were recorded. The heat release rate for the numerical simulations was calculated using 69 5.3. Engine Simulation Results the standard formula HRR(θ) = γ γ − 1p dV dθ + 1 γ − 1V dp dθ . (5.1) V is the engine volume and is a function of the crank angle θ (as is the pressure p); γ is the ratio of specific heats cp/cv and is taken to be 1.35. The heat release rate is a measure of the progress of combustion in the bulk region of the cylinder, and does not give much information about combustion at the boundaries. Figures 5.2 through 5.3 show similar comparisons for fuel consumption, power and emissions for both modes. The emissions comparison is not exact as the measurements are taken at different points in the cycle; the experi- mental results are for exhaust gases outside the cylinder, whereas the numer- ical results represent gases just before they leave the combustion chamber (at the opening of the exhaust valve). The experimental results for power represent brake mean effective pressure, while the numerical calculation re- turns the gross indicated mean effective pressure. Since the gross indicated mean effective pressure does not account for frictional losses, it is expected to be higher than brake mean effective pressure for both cases. The trend in emissions for both species are consistent with experiment, though in both cases is over- and under-predicted for the high and low load conditions respectively. This is partly explained by the cut off temperature for reactions in the TGLDM model; in this version the cut-off temperature is quite high, but the most current version recognizes that there will be con- tinued mixing as the exhaust valve is opened and gases leave the combustion chamber. This difference in TGLDM does not affect flame propagation or other properties to an appreciable degree, and was not considered during this study. 5.3.2 Numerical Quenching Results Presenting the results of the engine simulation and determining the extent and severity of quenching was a non-trivial task, as it must account for the spatial variation and time evolution of the flame and its interaction with the piston surface. To better understand the data, the overall combustion characteristics are presented for each mode (bulk-averaged in-cylinder pres- sure and heat release rate) as a function of crank angle, along with a brief discussion of what their evolution in time may mean for quenching. Four different visualizations will then be presented for each case: 70 5.3. Engine Simulation Results 0 5 10 15 20 M8 Base M9 Base M EP [ba r] Numerical Experimental (a) Mean Effective Pressure 0 50 100 150 200 250 M8 Base M9 Base IS FC [g/ kW .h] Numerical Experimental (b) Indicated Specific Fuel Consumption Figure 5.2: Fuel consumption and power comparison 71 5.3. Engine Simulation Results 0 0.2 0.4 0.6 0.8 1 1.2 1.4 M8 Base M9 Base N O x[g /kW .h] Numerical Experimental (a) NOx 0 0.5 1 1.5 2 M8 Base M9 Base UC H 4 [g/ kW .h] Numerical Experimental (b) Unburned Methane Figure 5.3: Emissions comparison 72 5.3. Engine Simulation Results • A time series of images showing density plots of ζ1 (the conserved scalar associated with methane transport). All images have a ζ1 range of [0:YCH4,st ], where YCH4,st is the mass fraction of methane present under stoichiometric conditions; as such regions where the flame has already interacted with the wall remain lighter coloured. • A time series of scatter plots of the Peclet number plotted against the progress variable Ypg1 (which tracks the progress of the reaction of methane) and φ. Each face of the piston surface where quenching occurs is shown with symbols whose colour and size depend on the relative amount of UHC in the quenching layer and Peb respectively. This view shows how quenching is related to the progress of combustion and the local fuel mixture. • A time plot of the amount of UHC in the quenching layer (MCH4,B) from the piston and from all wall surfaces, reported as a percentage of the total unburned fuel in the cylinder at the time. This shows how much UHC wall quenching contributes to engine out emissions. • A time plot of the quench distance, area-averaged over both all sur- faces combined and over the piston surface alone. The averages can be compared between the cases to identify or confirm trends in quench- ing with overall engine parameters. Only locations where the Peclet number provides a physically meaningful quench distance are included in calculation of the average. 5.3.3 Mode 8 Results Mode 8 Combustion Summary The bulk average pressure and heat release rate for the baseline Mode 8 cases are presented in Figure 5.4. The increased injection pressure has a significant impact on the heat release rate, which peaks earlier and at a larger value for the 600 bar case as compared to the 250 bar case; the extent of combustion is also changes, dropping sharply around 10°ATDC. From this information, it can be inferred that post-quench diffusion will be much more rapid and occur sooner in the cycle for higher injection pressures, and that the overall amount of UHC will be smaller. 73 5.3. Engine Simulation Results  0  50  100  150  200  250 -60 -40 -20  0  20  40  60 p  [ b a r ] CA [° ATDC] 250 Bar 600 Bar -100  0  100  200  300  400  500  600 -60 -40 -20  0  20  40  60 H R R  [ J / ° ] CA [° ATDC] 250 Bar 600 Bar Figure 5.4: Mode 8 combustion results 74 5.3. Engine Simulation Results Pinj = 250 bar The first case examined has engine parameters very similar to the validation in Chapter 4. The time series of images in Figure 5.5 show the injection of natural gas into the combustion chamber. There is noticeable impingement on the piston bowl, from which the greatest contribution to quenching from faces with lower Peb is expected. As the chamber expands the extent of quenching on the flat portion of the piston increases and an area of fuel trapped deep within the piston bowl is formed. From the initial validation section is known that Peb will be larger in that region, and a richer mixture there may result in a larger contribution of UHC from those faces. The scatter plots of Peb are shown in Figure 5.6. The results are con- sistent with the previous discussion while examining the ζ1 density plots, but also reveals new features. As combustion proceeds the cloud of points migrates to the right; as fuel is consumed and the air and fuel mix, the cloud of points fades from a rich mixture (φ > 1) to stoichiometry or leaner (φ ≤ 1). It is also clear quenching is occurring over several different ge- ometries (as Peb shows no clearly discernible trend); the seemingly uniform distribution over Ypg1 also reveals that combustion is still quite incomplete near the boundary for this operating condition. Because of the filtering of Peb the vertical distribution eventually reaches equilibrium as the flame has interacted with majority of the piston. The size of the points begins to di- minish after 70°ATDC, though whether this is due to an actual decrease in UHC in the quenching layer or an increase in UHC in the entire cylinder is not discernible from this visualization; information regarding wall quenching over the entire boundary is required, which shall be explored next. Figure 5.7 shows the amount of UHC in the quenching layer as a per- centage of the total UHC in the cylinder; because of the rapid decrease in UHC during the latter portion of the expansion, a logarithmic scale will be used. As expected, no wall quenching is seen at all until after TDC, after which a rapid increase in wall hydrocarbons is observed as the flame interacts with the piston wall. At its peak wall quenching on the piston contributes as much as 10% of the UHC in the piston, which is eventually reduced to around 2%. The total UHC from quenching follows a similar trajectory, though the final percentage stays around 10%, meaning that the other surfaces contribute a significant portion to wall quenching. Comparing to the results of Amano and Okamoto (the most relevant results available), the maximum peak contribution (not accounting for diffusion and oxidation of fuel after quenching) is of the same order of magnitude (near 2% of fuel input, compared to between 1-3%[4]); strictly speaking engines of the same 75 5.3. Engine Simulation Results (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.5: Mode 8 Pinj = 250 bar flame propagation 76 5.3. Engine Simulation Results  0  5  10  15  20  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  5  10  15  20  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.6: Mode 8 Pinj = 250 bar Peb distribution 77 5.3. Engine Simulation Results  0  5  10  15  20  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  5  10  15  20  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.6: Mode 8 Pinj = 250 bar Peb distribution (cont.) 78 5.3. Engine Simulation Results  0.0001  0.001  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.7: Mode 8 Pinj = 250 bar global quenching 79 5.3. Engine Simulation Results type and under the same conditions should be compared. The advantage of these results is that no assumptions about the spatial distribution or time evolution of δq were made, and the results for Peb and δq can be used on their own to evaluate a particular piston geometry or operating condition. The area weighted averages of the quench distance are shown in Figure 5.8. At under 400 µm (both the average on the piston and on the en- tire boundary), these values are comparable with the stoichiometric quench distance for atmospheric premixed flames. Since stoichiometric flames at the peak pressure should be roughly an order of magnitude lower[31], it is clear that the inhomogeneous mixing and strain rates have strong effects on quenching.  50  100  150  200  250  300  350  400  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.8: Mode 8 Pinj = 250 bar global quench distance Pinj = 600 bar An injection pressure of 600 bar is not typically used in production engines; though the specified injector is capable of handling such pressures, the in- jection system is rated to a maximum of 250 bar. Higher pressure injection was investigated as experimental evidence suggests it can enhance mixing, speed up the burning rate and reduce the formation of soot for diesel en- gines (tests have been performed with injection pressures of 500 bar[30] up 80 5.3. Engine Simulation Results to 2000 bar [29, 55]). The density plots of ζ1 for Pinj = 600 bar are shown in Figure 5.9. The degree of initial impingement is much greater than the 250 bar case, but it is unclear how quenching will ultimately be affected. There will undoubtedly be higher strain rates and improved mixing that could lower Peb and δq , but global quenching may still increase due to more flame-wall interaction. Additional information is obtained by visualizing Peb , shown in Figure 5.10. The scatter plots shows a migration of points towards Ypg1 = 1 and an overall lower average Peb . Though the overall Peb number appears to be lower, the contribution of wall quenching at 40°ATDC to the total amount is much higher as evidenced by the overall larger size of the points. However, post-quench diffusion and oxidization are also more pronounced for this case. The final amount of UHC from wall quenching may still be comparable to the 200 bar case; it is necessary to look at the total amount of wall quenching UHCs as it evolves over time (see Figure 5.11). Quenching from the piston contributes about two thirds of the UHC from wall quenching, which is unsurprising due to the amount of impingement that occurs on the piston surface. Though impingement has increased for this case, the total amount at the end of the stroke will be less; here the final contribution is 1%, compared to the 10% seen in the 250 bar case. Thus it appears a higher injection pressure increases the extent of quenching, but at the same time reduces the amount of fuel trapped in the quenching layer due to a thinner layer and faster reaction occurring in the bulk region (which promotes post-quench diffusion and oxidization). Examining the average quench distance in Figure 5.12, both averages are lower than the 250 bar injection pressure case. Figure 5.12 also illustrates how quickly quenching occurs under this op- erating condition. The average quench distance over the piston is roughly constant after 28°ATDC, more than 20 CA°before the 250 bar case; we also see the averages lie very close together, once again confirming that the quenching in this case is predominantly on the piston. 5.3.4 Mode 9 Cases Mode 9 Combustion Summary The baseline combustion results in Figure 5.13 show that even though much less fuel is being introduced into the combustion chamber, the resulting heat release rates are both similar in shape to the high injection pressure case for 81 5.3. Engine Simulation Results (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.9: Mode 8 Pinj = 600 bar flame propagation 82 5.3. Engine Simulation Results  0  2  4  6  8  10  12  14  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  2  4  6  8  10  12  14  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.10: Mode 8 Pinj = 600 bar Peb distribution 83 5.3. Engine Simulation Results  0  2  4  6  8  10  12  14  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  2  4  6  8  10  12  14  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.10: Mode 8 Pinj = 600 bar Peb distribution (cont.) 84 5.3. Engine Simulation Results  0.001  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.11: Mode 8 Pinj = 600 bar global quenching  80  100  120  140  160  180  200  220  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.12: Mode 8 Pinj = 600 bar global quench distance 85 5.3. Engine Simulation Results  0  10  20  30  40  50  60  70  80 -60 -40 -20  0  20  40  60 p  [ b a r ] CA [° ATDC] 150 Bar 250 Bar -50  0  50  100  150  200  250 -60 -40 -20  0  20  40  60 H R R  [ J / ° ] CA [° ATDC] 150 Bar 250 Bar Figure 5.13: Mode 9 combustion results 86 5.3. Engine Simulation Results Mode 8, though with much lower peak magnitudes. As such combustion is both rapid for these operating conditions it can be expected that the final amount of UHC will be significantly lower than for the high load cases, both as a percentage and as a raw amount. Pinj = 150 bar Both of the Mode 8 cases use relatively rich equivalence ratios (due to their high load). In Mode 9, less fuel is injected during combustion, leading to the possibility of larger δq (due to the lean flames with lower wall heat flux); at the same time lean mixtures may result in less fuel in the quenching layer. Figure 5.14, shows a high degree of impingement that occurs for this mode and injection pressure. At 10°ATDC, stagnation quenching occurs at the same location as it has for the Mode 8 cases; more interesting is the sidewall quenching on the flat surface of the piston; it initially seems very similar to the Mode 8, 600 bar injection pressure case, with the flame penetrating in the narrow region between the piston and the cylinder head. From the Peb scatter plots in Figure 5.15 it is clear sidewall quenching is more prominent under this operating condition (particularly at 40°ATDC). The filtering process keeps the final quenching Peb constant after QW has peaked, so the proportion of large to small Peb remains constant throughout the remainder of the stroke. Note the dramatic decrease in the amount of fuel in the quenching layer that occurs afterwards; in particular, a rapid decrease in the contribution from the piston surface to the total amount of UHC is seen as as the reaction proceeds to completion. The trend for wall UHC is markedly different for this mode compared to both Mode 8 cases; whereas the wall UHC from Mode 8 appear to be distributed over the entire chamber surface, here the majority of wall UHC comes from the piston surface. Also of importance is the subsequent diffusion of the wall UHC on the piston; from Figure 5.16 it is clear the rate at which diffusion occurs is much faster. This is likely due to two related reasons. First, the overall lower pressure means that the mass diffusion coefficient will be larger. Second, the reaction in the cell center proceeds rapidly, meaning there will be a larger gradient at the quench layer interface, increasing the rate at which fuel trapped in the quench layer diffuses out. The final amount of wall hydrocarbons from the the piston is 0.14%, over ten times less than the Mode 8 250 bar injection pressure case. Even on a logarithmic scale it is difficult to see any difference between the contribution from all surfaces and that from only the piston; thus it appears that for a leaner combustion case at lower load the piston contributes the majority of wall UHC. The actual 87 5.3. Engine Simulation Results (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.14: Mode 9 Pinj = 150 bar flame propagation 88 5.3. Engine Simulation Results  0  2  4  6  8  10  12  14  16  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  2  4  6  8  10  12  14  16  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.15: Mode 9 Pinj = 150 bar Peb distribution 89 5.3. Engine Simulation Results  0  2  4  6  8  10  12  14  16  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  2  4  6  8  10  12  14  16  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.15: Mode 9 Pinj = 150 bar Peb distribution (cont.) 90 5.3. Engine Simulation Results  0.0001  0.001  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.16: Mode 9 Pinj = 150 bar global quenching  80  100  120  140  160  180  200  220  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.17: Mode 9 Pinj = 150 bar global quench distance 91 5.3. Engine Simulation Results contribution to total UHC is still smaller than at high load. Figure 5.17 shows the quench distance averages, which are lower than both of the Mode 8 Cases analyzed so far. Piston wall quenching is most important here; the two traces follow the same trend and are very close together. Pinj = 250 bar The final case examined for normal timing was Mode 9 with a higher injec- tion pressure. Based on the results from all preceding cases, it is expected that the Peb distribution will be clustered around lower values, while the contributions to wall and global UHC should be lower if the trend observed for increased injection pressure continues. From the density plots in Figure 5.18 the flame continues to enter the region between the flat portion of the piston and the cylinder head, though the interaction with the cylinder head surface appears stronger. The scatter plots in Figure 5.19 show that there are more faces on which the flame wall interaction is primarily sidewall, though the magnitude is smaller than that seen when Pinj = 150 bar. This case shows an even greater amount of post-quench diffusion that the previous case. At 40°ATDC the contribution to global UHC is already greatly diminished. After determining that there was little change to the global pressure distribution, the only remaining explanation is the grouping of the point cloud. The reaction near the piston surface proceeds quickly (similar to the Mode 9 150 bar case), thus encouraging the wall UHC to diffuse back into the bulk region. For both Mode 9 cases, the peak contribution of wall surfaces to the total UHC is larger than the Mode 8 cases (∼ 50% vs. ∼ 1%); however, the final contribution is an order of magnitude lower; the UHC contribution here is 0.1% at the end of simulation. This case also continues the trend of the piston surface contributing the most to the wall UHC. Figure 5.21 shows the average quench distances, and the trend of higher contribution from the piston and lower quenching distances in higher pres- sure injection is continued. One difference is that even though the trend for both averages is similar, there is a larger difference than in the low injection pressure case. This could be explained by the earlier analysis of ζ1 where we see the flame enter the flat region between the piston and cylinder head surfaces. 92 5.4. Engine Timing Changes (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.18: Mode 9 Pinj = 250 bar flame propagation 93 5.4. Engine Timing Changes  0  1  2  3  4  5  6  7  8  9  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  1  2  3  4  5  6  7  8  9  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.19: Mode 9 Pinj = 250 bar Peb distribution 94 5.4. Engine Timing Changes  0  1  2  3  4  5  6  7  8  9  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  1  2  3  4  5  6  7  8  9  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.19: Mode 9 Pinj = 250 bar Peb distribution 95 5.4. Engine Timing Changes  0.001  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.20: Mode 9 Pinj = 250 bar global quenching  90  100  110  120  130  140  150  160  170  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.21: Mode 9 Pinj = 250 bar global quench distance 96 5.4. Engine Timing Changes 5.4 Engine Timing Changes As a final assessment of UHC emissions originating from cold wall quench- ing on the piston surface, the timing of injection was adjusted with both late and early injection. This would allow the geometry of interaction be- tween the flame and the piston surface to change while keeping the physical geometry constant. Changing the timing also affects other variables (peak pressure, temperature, etc.), so the effects of flame wall geometry cannot be completely isolated. 5.4.1 Normal Injection Pressure Normal Injection Pressure Timing Sweep Combustion Summary A comparison of the pressure and heat release rate for all the different tim- ings at normal pressure is shown in Figure 5.22. The change in timing does not affect the peak magnitude of the heat release rate so much as shift the location and how long it is sustained. Most interesting is that retarding the timing actually increases the heat release from later combustion. This alone cannot be taken to indicate better combustion, as the pressure data clearly shows a much lower cylinder pressure (and likely less efficient combustion). Mode 8 Pinj = 250 bar Advanced Timing Figure 5.23 shows a clear difference from normal injection; the initial im- pingement of the jet/flame on the piston results in an even larger rich region of fuel being concentrated in the lower bowl. This region typically experi- ences incomplete combustion, thus a large δq in this region is expected; there is also more interaction at the piston liner. However, the peak cylin- der pressure here is about 10% higher than the normal injection case, which can result in smaller δq . The distribution of Peb shown in Figure 5.24 has a larger maximum value, but the amount of fuel in the quenching layer is lower throughout the entire interaction. The distribution of points appears similar to the normal injection timing, though there are two distinct groupings (a vertical grouping along Ypg1 = 0.5, the other a triangular grouping whose topmost edge is along Peb = 18) appear at 70°ATDC and move together for the remainder of the simulation. The vertical grouping is associated with the piston bowl, while the triangular one represents points on the flat region of the piston, pointing to the increased coupling of geometry and quenching in this case. 97 5.4. Engine Timing Changes  0  20  40  60  80  100  120  140  160  180  200 -60 -40 -20  0  20  40  60 p  [ b a r ] CA [° ATDC] 250 Bar 250 Bar τa250 Bar τr -50  0  50  100  150  200  250  300  350 -60 -40 -20  0  20  40  60 H R R  [ J / ° ] CA [° ATDC] 250 Bar 250 Bar τa250 Bar τr Figure 5.22: Mode 8 altered timing combustion results 98 5.4. Engine Timing Changes (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.23: Mode 8 Pinj = 250 bar τa flame propagation 99 5.4. Engine Timing Changes  0  5  10  15  20  25  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  5  10  15  20  25  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.24: Mode 8 Pinj = 250 bar τa Peb distribution 100 5.4. Engine Timing Changes  0  5  10  15  20  25  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  5  10  15  20  25  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.24: Mode 8 Pinj = 250 bar τa Peb distribution 101 5.4. Engine Timing Changes  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.25: Mode 8 Pinj = 250 bar τa global quenching 102 5.4. Engine Timing Changes The amount of UHC throughout combustion in Figure 5.25 confirms the reduced contribution of the piston; the peak and final percentages are 26 and 2.4% respectively, less than half that of normal injection. It is clear that diffusion proceeds far more rapidly for the fuel in the quench layer than normal injection, with the contribution rapidly falling after 20°ATDC. The peak cylinder pressure for early injection case is nearly 10% higher than the normal injection case. The higher pressure cases discussed thus far have a slower reduction of boundary fuel. In light of the discussion of Figure 5.24 and examining the scatter plots on a finer timescale it was found that locations on the piston where fuel was initially trapped and Peb is large are both adjacent to cells where the reaction proceeds rapidly. Thus the increased diffusion is due to the large concentration gradient in fuel. This suggests that using the piston to redirect the flame does not have to incur a UHC penalty as long as combustion in the bulk region is complete. The average quenching distances in Figure 5.26 show the piston average to be similar to the normal injection timing, but that the contribution from all other surfaces becomes much more important; the average doesn’t reach a stable value during the simulation time, meaning that wall quenching is still occurring in other parts of the cylinder.  100  150  200  250  300  350  400  450  500  550  600  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.26: Mode 8 Pinj = 250 bar τa global quench distance 103 5.4. Engine Timing Changes Mode 8 Pinj = 250 bar Retarded Timing When retarding the timing it is expected that less flame wall interaction will be observed along with a correspondingly smaller contribution to UHC from the piston. Figure 5.27 shows that with later injection the flame/wall interaction has shifted from the region just below the bowl to the flat region of the piston. This is expected to increase the amount of sidewall quenching. From the density plots in Figure 5.27, the rich region in the expanding part of the cylinder has moved higher up and is closer to the cylinder head boundary. With this shift the other regions where cold wall quenching occurs will contribute more to the total amount of UHC. The scatter plots in Figure 5.28 show the effect of late injection on Peb ; the image for 10°ATDC shows that no quenching has occurred yet. On average the distribution is still very similar to the normal injection case, though there exists the tendency for the points to converge into two groups (a vertical grouping at Ypg1 = 0.5 and a horizontal triangular region whose lower edge is along Peb = 1 − 2). The migration of points to Ypg1 = 1 is closer to the normal injection case, indicating that the incomplete combustion may lead to more fuel trapped in the quenching layer at the end of the stroke. Observation of the scatter plots for all recorded time steps indicate that the main interaction ends at 40-50°ATDC, with diffusion proceeding for the remainder of the stroke. The amount of fuel trapped in the quenching layer appears to be larger at the end of quenching, but the final amount can be most clearly seen from Figure 5.29. The amount of UHC from wall quenching has a larger contribution from surfaces other than the piston, though the later decay appears to be dominated by the piston surface, as the two traces are virtually parallel to each other. This may be explained by the fact that in the retarded injection case the peak cylinder pressure is 20% less than during normal injection, which results in a higher diffusion coefficient overall. The larger distribution of points near Ypg1 also appears to confirm that the late cycle combustion does little to enhance burning at the wall. Comparing the average quenching distances from Figure 5.30 with the previous timings, the piston average is only slightly lower than both the advanced and normal timings. The average over all surfaces exhibits a new behaviour where there is a maximum near 50°ATDC. Though not entirely clear, the author believes this to be a result of the later interaction of the flame with piston liner. 104 5.4. Engine Timing Changes (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.27: Mode 8 Pinj = 250 bar τr flame propagation 105 5.4. Engine Timing Changes  0  5  10  15  20  25  30  35  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  5  10  15  20  25  30  35  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.28: Mode 8 Pinj = 250 bar τr Peb distribution 106 5.4. Engine Timing Changes  0  5  10  15  20  25  30  35  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  5  10  15  20  25  30  35  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.28: Mode 8 Pinj = 250 bar τr Peb distribution 107 5.4. Engine Timing Changes  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.29: Mode 8 Pinj = 250 bar τr global quenching  0  50  100  150  200  250  300  350  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.30: Mode 8 Pinj = 250 bar τr global quench distance 108 5.4. Engine Timing Changes  0  50  100  150  200  250 -60 -40 -20  0  20  40  60 p  [ b a r ] CA [° ATDC] 600 Bar 600 Bar τr -100  0  100  200  300  400  500  600  700 -60 -40 -20  0  20  40  60 H R R  [ J / ° ] CA [° ATDC] 600 Bar 600 Bar τr Figure 5.31: Mode 8 high injection pressure, altered timing combustion results 109 5.4. Engine Timing Changes 5.4.2 High Injection Pressure High Injection Pressure Timing Sweep Combustion Summary The comparison of the high load, high injection pressure cases in Figure 5.31 show less sensitivity to timing in the heat release rate (though both the magnitude and duration actually increase with injection delay), but a sub- stantial difference in the peak cylinder pressure. As seen earlier, the increase in heat release rate may not improve combustion near the wall, though with a higher injection pressure the possibility exists that the improvement in combustion extends to the near-wall region. The preceding analysis has shown that changing the geometry of flame wall interaction (in this case achieved by altering the injection timing) can have a significant impact on the amount of UHC cause by cold wall quench- ing. Although many of the trends still apply (e.g. more diffusion for lower pressure, higher injection pressure leading to lower Peb and quenching dis- tance), there may be unanticipated results caused by the coupled effects of the flame wall interaction, fuel/air mixing and combustion in the bulk region. Mode 8 Pinj = 600 bar Retarded Timing Figure 5.32 shows that with the change in injection timing, the flame wall interaction is now mostly located near the bowl of the piston, similar to the 250 bar case. As the piston expands the two rich regions that are typical of this geometry continue to form (one in the piston bowl and one in the expanding region above the flat portion of the piston). Based on the scale over which the density plot is produced, the mixture is lean in the majority of the combustion chamber, which coupled with a 20% lower peak pressure may lead to larger Peb . It is quite surprising to see that the Peb distribution is far lower than any of the cases seen thus far. A greater proportion of points are on the higher end of the scale early on in the expansion stroke (evident in the 40°ATDC image), though the amount clearly diminishes quickly. Comparing the Peb distributions with the normal timing case, the difference in the bowl of the piston is largely unaffected; rather the late injection means that the interaction between the flame and the piston over the flat portion has changed. From examination of the time series of scatter plots in Figure 5.33, it was observed that the distribution of Peb was more uniform in the late injection case, so that the higher maximum Peb in the normal injection is confined 110 5.4. Engine Timing Changes (a) 10°ATDC (b) 40°ATDC (c) 70°ATDC (d) 100°ATDC Figure 5.32: Mode 8 Pinj = 600 bar τr flame propagation 111 5.4. Engine Timing Changes  0  1  2  3  4  5  6  7  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (a) 10°ATDC  0  1  2  3  4  5  6  7  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (b) 40°ATDC Figure 5.33: Mode 8 Pinj = 600 bar τr Peb distribution 112 5.4. Engine Timing Changes  0  1  2  3  4  5  6  7  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (c) 70°ATDC  0  1  2  3  4  5  6  7  0.4  0.5  0.6  0.7  0.8  0.9  1 P e b Ypg1  0  0.5  1  1.5  2 φ (d) 100°ATDC Figure 5.33: Mode 8 Pinj = 600 bar τr Peb distribution 113 5.4. Engine Timing Changes  0.001  0.01  0.1  1  10  100  0  20  40  60  80  100  120  140 M C H 4 , B [ %  T o t a l  U H C ] CA [° ATDC] UHC - All Surfaces UHC - Piston Figure 5.34: Mode 8 Pinj = 600 bar τr global quenching 114 5.5. Summary of Results to a small region at the corner where the flat portion of the piston meets the liner; this is not surprising as it is known that two surface quenching typically exhibit higher Peb than single surface configurations[17, 33]. Although the total contribution to quenching to UHC seen in Figure 5.34 is comparable to the normal timing case, the difference between the final contribution of the piston is strikingly different; at 0.08% is is nearly four times less than the normal injection case. The overall cylinder pressure is lower in the late injection case, which typically leads to a larger mass diffusion coefficient. Examining the difference between YCH4,B and YCH4 over the piston surface revealed a much higher difference than during normal injection which leads to more diffusion of fuel out of the quenching layer. The average quench distance on the piston and remainder of the combustion chamber surfaces for this case is shown in Figure 5.35. The same trend of late interaction with piston liner that was seen with the 250 bar retarded injection pressure case occurs here as well.  0  50  100  150  200  250  300  350  20  40  60  80  100  120  140 δ q  [ µ m ] CA [° ATDC] Average - All Surfaces Average - Piston Figure 5.35: Mode 8 Pinj = 600 τr bar global quench distance 5.5 Summary of Results For each individual case the final amount of UHC was reported as a per- centage of the total. To enable comparison between cases the amount of 115 5.5. Summary of Results UHC from wall quenching was calculated and is collected in Figure 5.36. The four Mode 8 cases exhibit both the most UHC contribution and the largest changes with operating conditions. As such it is recommended any optimization involving changes to the operating conditions and combustion chamber geometry focus on the results of Mode 8.  1e-05  0.0001  0.001  0.01  0.1  1 Advanced Normal Retarded U H C [ m g / k W  h ] Timing Mode 8 Pinj 250 barMode 8 Pinj 600 barMode 9Pinj 150 barMode 9 Pinj 250 bar Figure 5.36: Unburned hydrocarbon comparison In the end, it appears that for HPDI engines the contribution of wall hydrocarbons near cycle end is 10% or less of total UHC emissions, with the 2% observed for the Mode 8 normal injection timing being the highest contribution. As a percentage of the total fuel injected, this is roughly 0.1%, with the values for Mode 9 being even lower (0.0005% of injected fuel or less). The usefulness of this model is not only in providing a quantitative estimate of the amount of UHC from wall quenching, but also in the spatial- and time evolution of the UHC. From the normal injection cases it was observed that the piston geometry can become important at higher loads; from the altered timing cases it is clear that flame/wall interaction can not only impact the quenching distance but the subsequent post-quench diffusion through the mixing of fuel and oxidizer in the bulk region. There still remain more detailed questions about exactly how the shape 116 5.5. Summary of Results of the piston can affect mixing and combustion by the breakup or redi- rection of the fuel jet[2, 6, 28, 49, 57]. The results of this work can be used to determine the usefulness of such measures, as the reduction in UHC from changing the combustion properties may be larger or smaller than the increase in wall UHC. By having a clear estimate of the unburned hydrocar- bons from wall quenching, engine designers can determine if the trade off in UHC is justifiable. The previous quenching model used for simulations of this type assumed a uniform quenching layer over all surfaces (regardless of flame wall inter- action), and a similar diffusion/oxidation model that has no lower cutoff temperature. The thickness of the uniform quenching layer was initially selected based on a uniform wall temperature of 500 K using Cleary and Farrell’s correlation for premixed flames during sidewall quenching and sub- sequently adjusted. The final amount of UHC predicted by this model from cases under the same conditions as this study are compared in Figure 5.37.  1e-05  0.0001  0.001  0.01  0.1  1 Uniform Thermal Model U H C [ m g / k W  h ] Models Mode 8 Pinj 250 barMode 8 Pinj 600 barMode 9Pinj 150 barMode 9 Pinj 250 bar Figure 5.37: Quench model comparison From the figure it is clear that the order of magnitude of UHC prediction is consistent for both models, but the improved components of the thermal model have different effects on different operating modes. The prediction of 117 5.5. Summary of Results the Mode 8 cases has increased while that of the Mode 9 cases has decreased. This is consistent with the results from previous sections that establish the rapid reduction of boundary fuel in the lower load cases and the relative magnitudes of the average quench distance. What is surprising is the simi- larity of the results given that the average quench distance from the thermal model is less than half of the selected uniform quench distance. 118 Chapter 6 Conclusions The effects of wall quenching and piston geometry on unburned hydrocar- bon emissions from a high pressure, direct injection natural gas engine were investigated using the a thermal model based on the work of Boust et. al. which was adapted for use in CFD simulations. This model was chosen as a survey of the available literature on single wall quenching found that although quenching is a complex phenomenon involving the inhibition of in- termediate chemical reactions, thermo- and fluid dynamics interaction and differential diffusion, it can be approximately modelled using energy con- servation and thermal arguments. To properly estimate the amount of fuel trapped in the quenching layer, it is necessary to include a model for the diffusion and oxidation of wall UHC that occurs subsequent to flame wall interaction. The thermal model was first evaluated at engine relevant pressures, in a range from 1 to 40 bar using numerical simulation techniques. The quench- ing distances for stoichiometric premixed flames of methane and heptane were each determined by both observing the temperature profile of the near- wall region during quenching, and by post-processing the simulation results with the thermal model. Though the peak heat flux was under predicted, it did follow a non-linearly increasing trend similar to that seen in the litera- ture. The results from the model were closer to experiment than those from post-processing the temperature profile. Strained diffusion flames were also simulated to see if numerical model could accurately predict the dependence of the peak heat flux and strain rate. The Peclet number number Peb was calculated using the Boust model during run-time and was found to give results with an expected order of magnitude, though direct comparison is difficult as experimental observations of the quench distance for strained diffusion flames are not currently available. The numerical implementation of the thermal model was further refined by simulating an engine case and observing the effects on quantities relating to the calculation of wall heat flux, flame power, Peclet number and quench distance (QW , QΣ , P eb and δq respectively). From these results a filtering algorithm was developed to calculate Peb and δq at the moment of quench- 119 6.1. Future Work and Applications ing. The spatial distribution and time evolution of Peb and δq over the piston surface was examined, with higher Peb and δq found at locations of sidewall quenching when compared to head-on quenching. A two-zone mass exchange model was implemented to account for the UHC trapped in the quench layer and was found to reduce the concentration of UHC near walls consistent with the transport properties of the mixture. The thermal model was used to determine the amount of UHC from the piston surface for several different cases. Both high- and low load conditions at the same engine speed and different injection pressures were evaluated; the effect of injection timing was also investigated with one of the high load cases. In general, it was found that the final contribution of UHC from quenching on the piston was around 10% for the high load case and less than 1% for low load cases (at normal injection timing), and the higher injection pressures make the contribution from the piston surface more im- portant while reducing total overall wall UHC. Changing the timing resulted in altered combustion characteristics, fuel-air mixing and flame-wall inter- action. Through the thermal model it was found that these coupled effects can influence the final contribution of UHC from wall quenching in unex- pected ways, and that results from the model can be helpful in engine design applications. 6.1 Future Work and Applications The results obtained by modelling show that attempts at reducing UHC emissions from natural gas HPDI engines by controlling wall quenching can only result in marginal improvements, certainly not enough to warrant full scale optimization of the piston geometry to reduce wall quenching. With the new quenching model in place, it is now possible to evaluate approaches that seek to improve fuel-air mixing by changing the flame-wall interaction, since it is possible to quantify the amount of increase or decrease in UHC from wall quenching and compare it to the changes due to changes in oper- ating conditions and combustion chamber geometry. Thus the new thermal quenching model is a useful tool for preliminary engine design using compu- tational methods. Because this quenching model is not tied to any particular combustion model, it is also possible to use this model in situations where the mixture is premixed or partially premixed, not just in HPDI combustion. One of the weaknesses of the model as it is currently implemented is that the wall heat flux QW is mesh dependent, which makes it difficult to validate it with new and different combustion models. A way to improve 120 6.1. Future Work and Applications the performance and widen the applicability of the model is to find ways to decrease the dependence of QW on mesh size, possibly by finding or developing near-wall models to account for the large temperature gradients that occur at the wall and being more sensitive to the wall temperature. For this study, the laminar flame speed was used to calculate the flame power and thermal flame thickness. It would be beneficial to compare this formulation to one using the turbulent flame speed; though the literature supports the use of the laminar flame speed under typical combustion condi- tions (including turbulent combustion), this comparison can lead to further insight into the behaviour of turbulent flames quenching on cold surfaces. 121 Bibliography [1] Exhaust gas hydrocarbons genesis and exodus. SAE Paper 620122, 6(620122):192, 1964. [2] A compound technology for hcci combustion in a di diesel engine based on the multi-pulse injection and the bump combustion chamber. SAE Paper 2003-01-0741, (2003-01-0741), 2003. SAE 2003 World Congress. [3] Andrew A. Adamczyk and George A. Lavoie. Laminar head-on flame quenching-a theoretical study, 1978. [4] Toshiji Amano and Kazuhisa Okamoto. Uuburned hydrocarbons emis- sion source from engines. SAE Paper 2001-01-3528, SP-1644(2001-01- 3528), 1 2001. SAE International Fall Fuels and Lubricants Meeting and Exhibition. [5] J. Andrae, P. Björnbom, and L. Edsberg. Numerical studies of wall effects with laminar methane flames. Combustion and Flame, 128(1- 2):165 – 180, 2002. [6] Rolf Egnell Andreas Vressner and Bengt Johansson. Combustion cham- ber geometry effects on the performance of an ethanol fueled hcci en- gine. SAE Paper 2008-01-1656, (2008-01-1656), 2008. SAE Interna- tional Powertrains, Fuels and Lubricants Congress. [7] C. Angelberger, T. Poinsot, and B. Delhaye. Improving near-wall com- bustion and wall heat transfer modeling in si engine computations. SAE Paper 972881, (972881), 1997. SAE International Fall Fuels and Lubri- cants Meeting and Exhibition. [8] M. Bellenoue, T. Kageyama, S. A. Labuda, and J. Sotton. Direct mea- surement of laminar flame quenching distance in a closed vessel. Ex- perimental Thermal and Fluid Science, 27(3):323 – 331, 2003. [9] R.J. Blint and J.H. Betchel. Hyrdrocarbon combustion near a cooled wall. SAE Paper 820063, (820063), 1982. SAE International Congress and Exposition. 122 Bibliography [10] B. Boust, J. Sotton, S.A. Labuda, and M. Bellenoue. A thermal formu- lation for single-wall quenching of transient laminar flames. Combustion and Flame, 149(3):286 – 294, 2007. [11] Bastien Boust, Julien Sotton, and Marc Bellenoue. Unsteady heat transfer during the turbulent combustion of a lean premixed methane- air flame: Effect of pressure and gas dynamics. Proceedings of the Combustion Institute, 31(1):1411 – 1418, 2007. [12] G. Bruneaux, T. Poinsot, and J.H. Ferziger. Premixed flame/wall inter- action in a turbulent channel flow: budget for the flame surface density evolution equation and modelling. Journal of Fluid Mechanics, 349(- 1):191–219, 1997. [13] David J Cleary and Patrick V Farrell. Single-surface flame quench- ing distance dependence on wall temperature, quenching geometry and turbulence, 1995. SAE International Congress and Exposition. [14] W.A. Daniel. Flame quenching at the walls of an internal combustion engine. Symposium (International) on Combustion, 6(1):886 – 894, 1957. Sixth Symposium (International) on Combustion. [15] A. De Lataillade, F. Dabireau, B. Cuenot, and T. Poinsot. Flame/wall interaction and maximum wall heat fluxes in diffusion burners. Pro- ceedings of the Combustion Institute, 29(1):775–779, 2002. [16] O. Ezekoye, R. Greif, and R.F. Sawyer. Increased surface temperature effects on wall heat transfer during unsteady flame quenching. Sympo- sium (International) on Combustion, 24(1):1465 – 1472, 1992. [17] P.W. Fairchild, R.D. Fleeter, and F.E. Fendell. Raman spectroscopy measurements of flame quenching in a duct-type crevice. Symposium (International) on Combustion, 20(1):85 – 90, 1985. Twentieth Sym- posium (International) on Combustion. [18] X. J. Gu, M. Z. Haq, M. Lawes, and R. Woolley. Laminar burning velocity and markstein lengths of methane-air mixtures. Combustion and Flame, 121(1-2):41 – 58, 2000. [19] Ömer Gülder. Correlations of laminar combustion data for alternative s.i. engine fuels. SAE Papers 841000, (841000), 8 1984. West Coast International Meeting and Exposition. 123 Bibliography [20] M. Gummalla, D. G. Vlachos, and M. A. Delichatsios. Bifurcations and structure of surface interacting methane-air diffusion flames. Combus- tion and Flame, 120(3):333 – 345, 2000. [21] James Harrington, Sandeep Munshi, Costi Nedelcu, Patric Ouellette, Jeff Thompson, and Stewart Whitfield. Direct injection of natural gas in a heavy-duty diesel engine, May 2002. [22] J.B. Heywood. Internal combustion engine fundamentals. McGraw-Hill series in mechanical engineering. McGraw-Hill, 1988. [23] J. Huang and W.K. Bushe. Simulation of transient turbulent methane jet ignition and combustion under engine-relevant conditions using con- ditional source-term estimation with detailed chemistry. Combustion Theory and Modelling, 11:977 – 1008, 2007. [24] EcoPoint Incorporated. Emission test cycles:european stationary cycle. http://www.dieselnet.com/standards/cycles/esc.html, 2010. [25] Nobuhiko Ishikawa and Melvyn C. Branch. A simple model of transient thermal flame quenching, 1977. [26] E. W. Kaiser J. A. LoRusso, G. A. Lavoie. An electrohydraulic gas sampling valve with application to hydrocarbon emissions studies. SAE Paper 800045, 80(800045), 1980. Automotive Engineering Congress and Exposition. [27] S.A. Labuda M. Bellenoue J. Sotton, B.Boust. Head-on quenching of transient laminar flame: Heat flux and quenching distance measure- ments. Combustion Science and Technology, 177, 2005. [28] Gilles Bruneaux Julian T. Kashdan, Sylvain Mendez. On the origin of unburned hydrocarbon emissions in a wall-guided, low nox diesel combustion system. SAE Paper 2007-01-1836, (2007-01-1836), 1 2007. 2007 JSAE/SAE International Fuels and Lubricants Meeting. [29] Prashanth Karra and Song-Charng Kong. Diesel emission characteris- tics using high injection pressure with converging nozzles in a medium- duty engine, 04 2008. [30] Shinji Kobayashi, Takayuki Sakai, Toshio Nakahira, Masanori Komori, and Kinji Tsujimura. Measurement of flame temperature distribution in d.i. diesel engine with high pressure fuel injection, February 1992. 124 Bibliography [31] Sergei Labuda, Maxime Karrer, Julien Sotton, and Marc Bellenoue. Experimental Study of Single-Wall Flame Quenching at High Pressures. Combustion Science and Technology, 183(5):409–426, 2011. [32] Jordan K. Lampert, M. Shahjahan Kazi, and Robert J. Farrauto. Pal- ladium catalyst performance for methane emissions abatement from lean burn natural gas vehicles. Applied Catalysis B: Environmental, 14(3-4):211 – 223, 1997. [33] George A. Lavoie. Correlations of combustion data for s. i. engine calculations - laminar flame speed, quench distance and global reaction rates, 1978. [34] Guowei Li, Patric Ouellette, Silviu Dumitrescu, and Philip G. Hill. Optimization study of pilot-ignited natural gas direct-injection in diesel engines, October 1999. [35] OpenCFD Limited. Openfoam®- the open source computational fluid dynamics (cfd) toolbox. http://www.openfoam.com/, 2009. [36] G. P. McTaggart-Cowan, W. K. Bushe, S. N. Rogak, P. G. Hill, and S. R. Munshi. Injection parameter effects on a direct injected, pilot ignited, heavy duty natural gas engine with egr, October 2003. [37] G. P. McTaggart-Cowan, H. L. Jones, S. N. Rogak, W. K. Bushe, P. G. Hill, and S. R. Munshi. The effects of high-pressure injection on a compression–ignition, direct injection of natural gas engine. Journal of Engineering for Gas Turbines and Power, 129(2):579–588, 2007. [38] Sylvain Mendez, Julian T. Kashdan, Gilles Bruneaux, Benoist Thi- rouard, and Franck Vangraefschepe. Formation of unburned hydrocar- bons in low temperature diesel combustion. SAE Paper 2009-01-2729, 2(2):205–225, 2009. SAE Powertrains, Fuels and Lubricants Meeting. [39] Intergovernmental Panel on Climate Change. Climate Change 2007: Impacts, Adaptation and Vulnerability. Intergovernmental Panel on Cli- mate Change, 2008. [40] Thierry Poinsot and Denis Veynante. Theoretical and Numerical Com- bustion, Second Edition. R.T. Edwards, Inc., 2005. [41] T.J. Poinsot, D.C. Haworth, and G. Bruneaux. Direct simulation and modeling of flame-wall interaction for premixed turbulent combustion. Combustion and Flame, 95(1-2):118 – 132, 1993. 125 Bibliography [42] S. B. Pope and U. Maas. Simplifying chemical kinetics: trajectory- generated low-dimensional manifolds. FDA 93-11, 1993. [43] P Popp and M Baum. Analysis of wall heat fluxes, reaction mechanisms, and unburnt hydrocarbons during the head-on quenching of a laminar methane flame. COMBUSTION AND FLAME, 108(3):327–348, FEB 1997. [44] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes 3rd Edition: The Art of Sci- entific Computing. Cambridge University Press, New York, NY, USA, 3 edition, 2007. [45] C.D. Rakopoulos, G.M. Kosmadakis, and E.G. Pariotis. Critical eval- uation of current heat transfer models used in cfd in-cylinder engine simulations and establishment of a comprehensive wall-function formu- lation. Applied Energy, 87(5):1612 – 1630, 2010. [46] G. Rozenchan, D.L. Zhu, C.K. Law, and S.D. Tse. Outward propaga- tion, burning velocities, and chemical effects of methane flames up to 60 atm. Proceedings of the Combustion Institute, 29(2):1461 – 1470, 2002. [47] A.J. Smallbone, W. Liu, C.K. Law, X.Q. You, and H. Wang. Experi- mental and modeling study of laminar flame speed and non-premixed counterflow ignition of n-heptane. Proceedings of the Combustion In- stitute, 32(1):1245 – 1252, 2009. [48] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eite- neer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. S. Song, W. C. Gardiner, V. V. Lissianski, and Z. Qin. GRI-Mech 3.0. 1996. [49] Wanhua Su, Rongwen Lin, Hui Xie, and Shao xi Shi. SAE Paper 971616, (971616), 1997. SAE International Spring Fuels and Lubri- cants Meeting and Exposition. [50] S.R. Vosen, R. Greif, and C.K. Westbrook. Unsteady heat transfer during laminar flame quenching. Symposium (International) on Com- bustion, 20(1):75 – 83, 1985. Twentieth Symposium (International) on Combustion. [51] M. Wang, J. Huang, and W.K. Bushe. Simulation of a turbulent non- premixed flame using conditional source-term estimation with trajec- 126 Bibliography tory generated low-dimensional manifold. Proceedings of the Combus- tion Institute, 31(2):1701 – 1709, 2007. [52] J. Warnatz, U. Maas, and R.W. Dibble. Combustion: Physical and Chemical Fundamentals, Modeling and Simulation, Experiments, Pol- lutant Formation. Springer, 2001. [53] J. T. Wentworth. Effect of combustion chamber surface temperature on exhaust hydrocarbon concentration, 2 1971. [54] Charles K. Westbrook, Andrew A. Adamczyk, and George A. Lavoie. A numerical study of laminar flame wall quenching. Combustion and Flame, 40:81 – 99, 1981. [55] Haruyuki Yokota, Takeyuki Kamimoto, Hidenori Kosaka, and Kinji Tsujimura. Fast burning and reduced soot formation via ultra-high pressure diesel fuel injection, February 1991. [56] Paul Zelenka, Wolfgang Cartellieri, and Peter Herzog. Worldwide diesel emission standards, current experiences and future needs. Applied Catalysis B: Environmental, 10(1-3):3 – 28, 1996. [57] Y. Zhu and H. Zhao. Computational study of the effects of the throat diameter of the piston bowl for the performance and emissions of a high-speed direct-injection diesel engine. Proceedings of the I MECH E Part D Journal of Automobile Engineering, 220(1):111–124, 2006. 127 Appendix A Thermal Quench Model Presented below is a summary of the derivation of the thermal model by Boust et. al. that was used as the basis for thermal quenching model. A full discussion about the model and its applicability can be found in the paper “A thermal formulation for single-wall quenching of transient laminar flames”[10]. At the moment the flame is seen as a one dimensional structure; at the moment of quenching, the near the wall is partitioned into three distinct regions: • The reaction zone, where combustion is actively occurring, which is taken to be as thick as the flame. • The preheat zone, which is the region directly ahead of the flame propagation. The gases in this region are being heated by the flame, and the size of this region is on the order of the laminar flame thickness δl . • The quench layer, where no reaction occurs and the unburned fuel is located. This region is of thickness δq , which is the desired quantity. A schematic of these regions is reproduced from Boust’s paper, which also shows the temperature profile as the flame quenches on the cold wall in Figure A.1. The following subscripts are used to denote quantities at the following regions or conditions: • u: Unburned gases • ad: Adiabatic (referring to the adiabatic flame temperature) • f : Flame temperature tq represents the time at which flame quenching occurs. The first law of thermodynamics for a control volume at the flame can be written as Qu + ρuSlcp(Tf − Tu) = ρuSlYfuel∆H. (A.1) 128 Appendix A. Thermal Quench Model δq δl Tad T f TW Tu x QW Qu Qu ρSlc p(T f−Tu) ρSlΔH Y fuel t=tq : quenching t<tq : propagation PreheatZone ReactionZone Figure A.1: Simplified model of head-on quenching[10]. Reprinted and adapted from Combustion and Flame, Volume 149, Issue 3, B. Boust, J. Sotton, S.A. Labuda, M. Bellenoue, A thermal formulation for single-wall quenching of transient laminar flames, p.333, ©2007 Combustion and Flame with permission from Elsevier At the moment of quenching, this energy balance can be applied with the following approximations: • Quenching occurs when heat losses from the unburned region to the wall exceed those from the reaction zone to the unburned region: Qu = QW • The wall heat flux is estimated to first order by conduction through a thin slab of gas: QW = k Tf−TW δq • The unburned gas temperature approximately equal to the wall tem- perature: Tu ≈ TW 129 Appendix A. Thermal Quench Model Incorporating these approximations into (A.1) yields the equation QW + ρuSlcp QW δq k = ρuSlYfuel∆H. (A.2) Equation (A.2) can be easily manipulated to obtain the final form of the thermal relation: 1 + δq δl = QΣ QW (A.3) Peb = 1 Φ − 1 (A.4) A similar argument is made for sidewall quenching, with the assumption (and subsequent validation through Boust’s results) that the flame stretching near the wall only has second order effects on the quench distance. 130 Appendix B Numerical Model Details This section lays out the details of the OpenFOAM based numerical solver used for the majority of the results. It draws heavily from the OpenFoam user guide which is part of the standard OpenFOAM installation. The main details and structure of the solver will be discussed here in general terms, as the user guide and included source code documentation provide much more detail. In addition there are numerous tutorials available on the Internet that describe how to build solvers in OpenFOAM. Introduction to OpenFOAM At its core, OpenFOAM is a set of open source libraries and tools written in C++ than can be used to create numerical solvers and perform post- processing for many ordinary and partial differential equations. It bridges the gap between commercial solvers (which provide ease of use but lack opportunities for customization for non-standard problems) and custom- coded solvers that require extensive adaptation and can require a significant amount of development and validation before being usable. In grossly oversimplified terms, a numerical solver of any kind is required to perform the following tasks: • Discretize the solution domain • Initialize the solution variables • Perform a numerical calculation • Update solution variables • Enforce boundary conditions The last three steps are repeated until criteria set forth by the user are met. The previous steps are the absolute minimum a solver should perform; there may be additional steps involving post-processing, input-output operations and changes to the computational domain. OpenFOAM allows the rapid development of solvers by providing standard libraries with many well known 131 Appendix B. Numerical Model Details numerical schemes, thermophysical models and other various other useful model, as well as the ability to implement custom models as additional libraries. As an illustrative example, the solution of the Navier-Stokes equations is shown. This set of equations can be written as ∂ρU ∂t +∇ · φU−∇ · µ∇U = −∇p. (B.1) After the creation of and discretization of the domain and initialization of solution variables (both accomplished at run time by reading input files specified by the user) and loading the thermophysical and chemistry models. The computational time step is calculated based on the mesh and previously loaded models, as well as additional inputs specified by the user (the reader is pointed to the OpenFOAM User Guide for details). The actual code for solving the equations is rather simple, Program B.1 Code for solving Navier-Stokes equations 1 s o l v e 2 ( 3 fvm : : ddt ( rho , U) 4 + fvm : : div ( phi , U) 5 − fvm : : l a p l a c i a n (mu, U) 6 == 7 − f v c : : grad (p) 8 ) ; The code bears a striking resemblance to the mathematical form of the original equation, which is one of the key features of the OpenFOAM plat- form. Instead of spending their time implementing standard methods from scratch, the user is free work from a high-level of abstraction and focus on the physics. Should problems arise with the numerical models, the open- source nature of the OpenFOAM code allows inspection and modification of the relevant models and methods. Numerical Schemes and Models The solver built for this research was based on the reactingFoam program which is included in the standard OpenFOAM installation, with the built in chemistry solver only being used for slowly reacting species such as NOx. 132 Appendix B. Numerical Model Details The reaction of hydrocarbons is handled by an external Fortran implemen- tation of the TGLDM method, detailed in Huang and Bushe’s work[23]. To best summarize the numerical schemes and solution methods em- ployed for the solver, both the fvSchemes and fvSolution files are included in this section. The standard OpenFOAM header has been removed for bet- ter reproduction in print; Reference can be made to the OpenFOAM User Guide for interpretation of the token for a particular scheme or method. For the PISO method for solving for the pressure, the number of non- orthogonal correctors was varied between 1 and 3 dependent on the initial mesh quality. Several good references exist for details of each numerical schemes; one particularly salient work is Numerical Recipes: The Art of Scientific Computing [44]. 133 Appendix B. Numerical Model Details Program B.2 fvSchemes file 1 ddtSchemes 2 { 3 default Euler ; 4 } 5 6 gradSchemes 7 { 8 default Gauss l i n e a r ; 9 grad (p) Gauss l i n e a r ; 10 } 11 12 divSchemes 13 { 14 default Gauss l i m i t e d L i n e a r 1 ; 15 div ( phi , rho ) Gauss l i m i t e d L i n e a r 1 ; 16 div ( phi ,U) Gauss upwind ; 17 div ( phid , p) Gauss upwind ; 18 div ( phiU , p) Gauss l i n e a r ; 19 div ( phi , k ) Gauss upwind ; 20 div ( phi , e p s i l o n ) Gauss upwind ; 21 div ( phi , Yi ) Gauss l i m i t e d L i n e a r 1 ; 22 div ( phi , Ypg1 Gauss l i m i t e d L i n e a r 1 ; 23 div ( phi , Ypg2) Gauss l i m i t e d L i n e a r 1 ; 24 div ( phi , Yno) Gauss l i m i t e d L i n e a r 1 ; 25 div ( phi , h ) Gauss l i m i t e d L i n e a r 1 ; 26 div ( phi , zeta1 ) Gauss l i m i t e d L i n e a r 1 ; 27 div ( phi , pzeta ) Gauss l i m i t e d L i n e a r 1 ; 28 div ( ( muEff∗dev2 ( grad (U) .T( ) ) ) ) Gauss l i n e a r ; 29 div (U) Gauss l i n e a r ; 30 } 134 Appendix B. Numerical Model Details 1 lap lac ianSchemes 2 { 3 default Gauss l i n e a r l i m i t e d 0 . 5 ; 4 l a p l a c i a n ( muEff ,U) Gauss l i n e a r l i m i t e d 0 . 5 ; 5 l a p l a c i a n ( DkEff , k ) Gauss l i n e a r l i m i t e d 0 . 5 ; 6 l a p l a c i a n ( Deps i lonEf f , e p s i l o n ) Gauss l i n e a r l i m i t e d 0 . 5 ; 7 l a p l a c i a n ( muEff , zeta1 ) Gauss l i n e a r l i m i t e d 0 . 5 ; 8 l a p l a c i a n ( muEff , pzeta ) Gauss l i n e a r l i m i t e d 0 . 5 ; 9 l a p l a c i a n ( alphaEff , h ) Gauss l i n e a r l i m i t e d 0 . 5 ; 10 l a p l a c i a n ( muEff , Yi ) Gauss l i n e a r l i m i t e d 0 . 5 ; 11 l a p l a c i a n ( muEff , Ypg1) Gauss l i n e a r l i m i t e d 0 . 5 ; 12 l a p l a c i a n ( muEff , Ypg2) Gauss l i n e a r l i m i t e d 0 . 5 ; 13 l a p l a c i a n ( muEff , Yno) Gauss l i n e a r l i m i t e d 0 . 5 ; 14 l a p l a c i a n ( ( rho |A(U) ) ,p) Gauss l i n e a r l i m i t e d 0 . 5 ; 15 l a p l a c i a n ( rhoD , k ) Gauss l i n e a r l i m i t e d 0 . 5 ; 16 l a p l a c i a n ( rhoD , e p s i l o n ) Gauss l i n e a r l i m i t e d 0 . 5 ; 17 } 18 19 inte rpo la t i onSchemes 20 { 21 default l i n e a r ; 22 } 23 24 snGradSchemes 25 { 26 default l i m i t e d 0 . 5 ; 27 } 28 29 f luxRequ i red 30 { 31 p ; 32 } 135 Appendix B. Numerical Model Details Program B.3 fvSolution file 1 s o l v e r s 2 { 3 rho 4 { 5 s o l v e r PCG; 6 p r e c o n d i t i o n e r DIC ; 7 t o l e r a n c e 1e−05; 8 r e l T o l 0 ; 9 } ; 10 U 11 { 12 s o l v e r PBiCG; 13 p r e c o n d i t i o n e r DILU ; 14 t o l e r a n c e 1e−05; 15 r e l T o l 0 ; 16 } ; 17 p 18 { 19 s o l v e r PCG; 20 p r e c o n d i t i o n e r DIC ; 21 t o l e r a n c e 1e−06; 22 r e l T o l 0 ; 23 } ; 24 Yi 25 { 26 s o l v e r PBiCG; 27 p r e c o n d i t i o n e r DILU ; 28 t o l e r a n c e 1e−05; 29 r e l T o l 0 ; 30 } ; 31 Ypg1 32 { 33 s o l v e r PBiCG; 34 p r e c o n d i t i o n e r DILU ; 35 t o l e r a n c e 1e−05; 36 r e l T o l 0 ; 37 } ; 136 Appendix B. Numerical Model Details 1 Ypg2 2 { 3 s o l v e r PBiCG; 4 p r e c o n d i t i o n e r DILU ; 5 t o l e r a n c e 1e−05; 6 r e l T o l 0 ; 7 } ; 8 Yno 9 { 10 s o l v e r PBiCG; 11 p r e c o n d i t i o n e r DILU ; 12 t o l e r a n c e 1e−05; 13 r e l T o l 0 ; 14 } ; 15 Ych4 16 { 17 s o l v e r PBiCG; 18 p r e c o n d i t i o n e r DILU ; 19 t o l e r a n c e 1e−05; 20 r e l T o l 0 ; 21 } ; 22 Yco 23 { 24 s o l v e r PBiCG; 25 p r e c o n d i t i o n e r DILU ; 26 t o l e r a n c e 1e−05; 27 r e l T o l 0 ; 28 } ; 29 Yc2h2 30 { 31 s o l v e r PBiCG; 32 p r e c o n d i t i o n e r DILU ; 33 t o l e r a n c e 1e−05; 34 r e l T o l 0 ; 35 } ; 137 Appendix B. Numerical Model Details 1 h 2 { 3 s o l v e r PBiCG; 4 p r e c o n d i t i o n e r DILU ; 5 t o l e r a n c e 1e−05; 6 r e l T o l 0 ; 7 } ; 8 k 9 { 10 s o l v e r PBiCG; 11 p r e c o n d i t i o n e r DILU ; 12 t o l e r a n c e 1e−06; 13 r e l T o l 0 ; 14 } ; 15 e p s i l o n 16 { 17 s o l v e r PBiCG; 18 p r e c o n d i t i o n e r DILU ; 19 t o l e r a n c e 1e−06; 20 r e l T o l 0 ; 21 } ; 22 zeta1 23 { 24 s o l v e r PBiCG; 25 p r e c o n d i t i o n e r DILU ; 26 t o l e r a n c e 1e−05; 27 r e l T o l 0 ; 28 } ; 29 zeta2 30 { 31 s o l v e r PBiCG; 32 p r e c o n d i t i o n e r DILU ; 33 t o l e r a n c e 1e−05; 34 r e l T o l 0 ; 35 } ; 138 Appendix B. Numerical Model Details 1 pzeta1 2 { 3 s o l v e r PBiCG; 4 p r e c o n d i t i o n e r DILU ; 5 t o l e r a n c e 1e−05; 6 r e l T o l 0 ; 7 } ; 8 pzeta2 9 { 10 s o l v e r PBiCG; 11 p r e c o n d i t i o n e r DILU ; 12 t o l e r a n c e 1e−05; 13 r e l T o l 0 ; 14 } ; 15 } 16 PISO 17 { 18 nCorrectors 2 ; 19 nNonOrthogonalCorrectors 3 ; 20 momentumPredictor yes ; 21 fluxGradp no ; 22 } 139 Appendix C Quench Model Code The following are portions of code developed for this work, which execute are part of the larger solver for engine simulation. The variables have been named in a manner that should make their purpose self-evident; comments have been added to clarify the purpose of code that is specific to Open- FOAM. Quench Distance Code This code calculates Peb and δq for all faces of all boundaries that are designated as solid walls. Boundary Diffusion Code The following code the amount of fuel at the boundaries, assuming that wall quenching has occurred. 140 Appendix C. Quench Model Code Program C.1 qDist.H 1 // C a l c u l a t e the e q u i v a l e n c e r a t i o and 2 // l o c a l laminar f lame speed 3 equiRat = zeta1 / s t o i c h Y f u e l ; S l = uSl ( ) ( eqr ) ; 4 5 // C a l c u l a t e f lame power in the i n t e r i o r domain 6 f o r A l l ( flamePower . mesh ( ) . c e l l s ( ) , c e l l i ) 7 { 8 flamePower [ c e l l i ]=delH (T[ c e l l i ] ) ∗ rho [ c e l l i ]∗ Sl [ c e l l i ]∗ zeta1 [ c e l l i ] ; 9 } 10 11 // C a l c u l a t e the w a l l heat f l u x and copy over 12 // to a compat i b l e f i e l d 13 hFlux = −f v c : : i n t e r p o l a t e ( turbulence−>a lphaEf f ( ) ) ∗ f v c : : snGrad (h) ; 14 15 f o r A l l (qW. boundaryField ( ) , patch i ) 16 { 17 qW. boundaryField ( ) [ patch i ] = hFlux . boundaryField ( ) [ patch i ] ; 18 } 19 20 //For a l l boundar ies t h a t are w a l l s 21 f o r A l l ( wallPatchNames , wordi ) 22 { 23 // Find the current patch 24 l a b e l curPatch = mesh . boundaryMesh ( ) . f indPatchID ( wallPatchNames [ wordi ] ) ; 25 26 // Find the face−ID of the f i r s t f a c e in the patch 27 const fvPatch& p = patches [ curPatch ] ; 28 const polyPatch& pp = p . patch ( ) ; 141 Appendix C. Quench Model Code 29 //For a l l w a l l pa tches 30 f o r A l l ( hFlux . boundaryField ( ) [ curPatch ] , patch i ) 31 { 32 s c a l a r deltaLTemp = 0 . 0 ; 33 // I f the heat f l u x i s non−zero 34 i f ( hFlux . boundaryField ( ) [ curPatch ] [ patch i ] > QwMin) 35 { 36 // C a l c u l a t e f lame power 37 flamePower . boundaryField ( ) [ curPatch ] [ patch i ] = delH (T[ own [ pp . s t a r t ( )+patch i ] ] ) ∗ rho [ own [ pp . s t a r t ( )+patch i ] ] ∗ Sl [ own [ pp . s t a r t ( )+patch i ] ] ∗ zeta1 [ own [ pp . s t a r t ( )+patch i ] ] ; 38 39 // I f the curren t heat f l u x i s l a r g e r than the l a s t 40 i f ( hFlux . boundaryField ( ) [ curPatch ] [ patch i ] > qWMax. boundaryField ( ) [ curPatch ] [ patch i ] ) 41 { 42 // Update maximum w a l l heat f l u x 43 qWMax. boundaryField ( ) [ curPatch ] [ patch i ] = hFlux . boundaryField ( ) [ curPatch ] [ patch i ] ; 44 45 // Update maximum flame power 46 flamePowerMax . boundaryField ( ) [ curPatch ] [ patch i ] = max( flamePower . boundaryField ( ) [ curPatch ] [ patch i ] , flamePowerMax . boundaryField ( ) [ curPatch ] [ patch i ] ) ; 47 48 // C a l c u l a t e l o c a l thermal f lame t h i c k n e s s 49 i f ( S l [ own [ pp . s t a r t ( )+patch i ] ] > SMALL) 50 { 51 deltaLTemp = turbulence−>alpha ( ) [ own [ pp . s t a r t ( )+patch i ] ] / S l [ own [ pp . s t a r t ( )+patch i ] ] ; 52 } 142 Appendix C. Quench Model Code 53 // Update thermal f lame t h i c k n e s s 54 i f ( deltaLTemp != 0 .0 && deltaL . boundaryField ( ) [ curPatch ] [ patch i ] > SMALL) 55 { 56 de ltaL . boundaryField ( ) [ curPatch ] [ patch i ] = min ( de l taL . boundaryField ( ) [ curPatch ] [ patch i ] , min ( de l taL [ own [ pp . s t a r t ( )+patch i ] ] , deltaLTemp ) ) ; 57 } 58 else 59 { 60 de ltaL . boundaryField ( ) [ curPatch ] [ patch i ] = min ( de l taL [ own [ pp . s t a r t ( )+patch i ] ] , deltaLTemp ) ; 61 } 62 63 } 64 // C a l c u l a t e P e c l e t number accord ing to Boust 65 Peb . boundaryField ( ) [ curPatch ] [ patch i ] = flamePowerMax . boundaryField ( ) [ curPatch ] [ patch i ] /qWMax. boundaryField ( ) [ curPatch ] [ patch i ] − 1 ; 66 } 67 else 68 { 69 Peb . boundaryField ( ) [ curPatch ] [ patch i ] = −1; 70 } 71 // C a l c u l a t e quenching t h i c k n e s s 72 i f (Peb . boundaryField ( ) [ curPatch ] [ patch i ] > 0) 73 { 74 qDist . boundaryField ( ) [ curPatch ] [ patch i ] = Peb . boundaryField ( ) [ curPatch ] [ patch i ]∗ deltaL . boundaryField ( ) [ curPatch ] [ patch i ] ; 75 } 76 } 77 } 143 Appendix C. Quench Model Code Program C.1 updateBoundary.H 1 s c a l a r de l tax =0.0005; 2 s c a l a r D = 1.51 e−5; 3 s c a l a r Sc = 1 . 0 ; 4 totalMch4B = 0 ; 5 v o l S c a l a r F i e l d nearWall = wa l lD i s t ( mesh ) . y ( ) ; 6 7 f o r A l l ( wallPatchNames , wordi ) 8 { 9 l a b e l curPatch = mesh . boundaryMesh ( ) . f indPatchID ( wallPatchNames [ wordi ] ) ; 10 const fvPatch& p = patches [ curPatch ] ; 11 const polyPatch& pp = p . patch ( ) ; 12 13 f o r A l l (Ych4Bond . boundaryField ( ) [ curPatch ] , patch i ) 14 { 15 de l tax = qDist . boundaryField ( ) [ curPatch ] [ patch i ] ; 16 i f ( de l tax > 0 .0 && de l tax < nearWall [ own [ pp . s t a r t ( )+patch i ] ] ) 17 { 18 D = (1 . 0 / Sc/ rho . boundaryField ( ) [ curPatch ] [ patch i ] ) ∗ turbulence−>mu( ) . boundaryField ( ) [ curPatch ] [ patch i ] ; 19 s c a l a r vBond = mesh . boundary ( ) [ curPatch ] . magSf ( ) [ patch i ]∗ de l tax ; 20 s c a l a r dYch4Bond = −2∗D/( de l tax ∗ de l tax ) ∗( Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ]− Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] ) ∗dt ; 21 s c a l a r Ych4BondTemp = dYch4Bond+Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] ; 22 23 i f (Ych4BondTemp>1) 24 { 25 i f (Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] > Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] ) 26 Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] = min (Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] , Ych4BondTemp) ; 144 Appendix C. Quench Model Code 27 else 28 Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] = 1 ; 29 } 30 else i f (Ych4BondTemp<Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] ) 31 { 32 Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] = Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] ; 33 } 34 else 35 { 36 Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ]= Ych4BondTemp ; 37 } 38 39 Ych4Bond . boundaryField ( ) . boundaryInte rna lF i e ld ( ) [ curPatch ] [ patch i ] = Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] ; 40 totalMch4B = totalMch4B+Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ]∗ vBond∗ rho . boundaryField ( ) [ curPatch ] [ patch i ] ; 41 } 42 else 43 { 44 Ych4Bond . boundaryField ( ) [ curPatch ] [ patch i ] = Ych4 . boundaryField ( ) [ curPatch ] [ patch i ] ; 45 } 46 } 47 } 145 Appendix D Refinements to Model The following section represent two possible avenues for improvement of the quench model under general conditions. Effect of Strain Rate on Flame Power As cited earlier from DeLataillade’s work, the peak heat flux for a strained diffusion flame can be almost ' QΣ as usually defined[15]. From the Fig- ure 2.2(b), there is almost a linear correlation between the product of strain rate and flame time (atf ) and the peak heat flux and burning Peclet number (QW , P eb), at least in the range 0 ≤ atf ≤ 1. A simple linear correction to the QW seems possible, so that the new version of Boust model in dimen- sional variables may be written as δq δl = QΣ QW (atf + 1)− 1 (D.1) This form was chosen so that the original form of the Boust model is recovered when atf = 0. Also important to note that when atf = 1, then QW ' QΣ, which according to this relation should give δqδl ' 1, which is also desired. If the product atf > 1, then as QW < QΣ and the model will predict that the quenching distance begins to increase again. It is not exactly clear what truly happens after this point, since there is only one experiment where atf > 1, and as such predictions should not be taken as being indicative of the performance of the model until further experimental evidence is obtained. There also remains fact that the sidewall diffusion flame configuration has not been evaluated, and it is uncertain what effect (if any) strain may have on heat transfer. As the flow pattern is parallel to the quenching thickness, it seems at first that strain should not have a significant effect on the heat transfer. As such the above equation could also be amended by accounting for the direction of the flow field using a dot product. δq δl = QΣ QW (atfu · n̂ + 1)− 1 (D.2) 146 Appendix D. Refinements to Model where u is the velocity at the cell center and n̂ is the face normal of the bounding wall surface, defined as the outward normal to the domain. This way sidewall quenching is not unduly affected, while head on quenching has a larger effect recorded. Initial results from this amended model do give Peb consistent with the results of head on premixed stoichiometric flames, however more testing and analysis are required to see how this change affects the sidewall Peb. It would also first be necessary to determine the effect of strain on sidewall quenching. Gas Properties at Quenching Another way the model as it is currently implemented could be improved is in the tracking of gas properties during quenching. At the moment all quenching quantities are tracked and updated as long as the heat flux con- tinues to rise. This is assuming that a flame will quench once and only once on a particular wall. As the temperature, density and pressure all rapidly change, it is difficult to determine when gas properties are changing due to the incoming flame front and when they are changing due to chemical reac- tions at the near-wall cell. All of these have an as yet undetermined effect on QΣ At the moment, the opposite trends for density and temperature during quenching have been exploited, though no rigorous testing has been done to determine how much these effects alter the value of QΣ. One possible method of determining this is through the progress variable Ypg1. This can be seen while tracking the Ypg1 and QW at individual locations from the validation case, as shown in Figure D.1. Locations 0,1 and 2 all show similar trends in Ypg1; the values rises rapidly as the flame interacts with the wall, then subsequently falls as the spent fuel moves away from the wall. Location 3 has a more complex, with QW lagging Ypg1 slightly. It may be that quenching is occurring in a slightly different way, since location 3 is at the lip of the piston bowl where the flow is being redirected. Thus, is may be possible to also use Ypg1 as a measure of the approach of the flame front. Further analysis of the coupling between Ypg1 and QW should reveal if this is possible and the best approach for implementation. Due to time constraints this analysis was not performed. 147 Appendix D. Refinements to Model  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  0  5  10  15  20  25  30  0  5e+06  1e+07  1.5e+07  2e+07  2.5e+07  3e+07  3.5e+07  4e+07  4.5e+07  5e+07  5.5e+07 Y p g 1 Q W  [ W / m 2 ] CA [° ATDC] Ypg1 QW (a) Location 0  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  0  5  10  15  20  25  30  0  1e+06  2e+06  3e+06  4e+06  5e+06  6e+06  7e+06 Y p g 1 Q W  [ W / m 2 ] CA [° ATDC] Ypg1 QW (b) Location 1 Figure D.1: Ypg1 and QW at quenching locations 148 Appendix D. Refinements to Model  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  0  5  10  15  20  25  30  0  5e+06  1e+07  1.5e+07  2e+07  2.5e+07  3e+07  3.5e+07 Y p g 1 Q W  [ W / m 2 ] CA [° ATDC] Ypg1 QW (c) Location 2  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  0  5  10  15  20  25  30  0  2e+06  4e+06  6e+06  8e+06  1e+07  1.2e+07  1.4e+07  1.6e+07  1.8e+07  2e+07  2.2e+07 Y p g 1 Q W  [ W / m 2 ] CA [° ATDC] Ypg1 QW (d) Location 3 Figure D.1: Ypg1 and QW at quenching locations (cont.) 149

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072348/manifest

Comment

Related Items