Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

CFD modeling of injection strategies in a high-pressure direct-injection (HPDI) natural gas engine Kheirkhah, Pooyan 2015

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

Item Metadata

Download

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

Full Text

CFD modeling of injection strategies in a  High-Pressure Direct-Injection (HPDI) natural gas engine  by  Pooyan Kheirkhah  BSc., University of Tehran, 2012   A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE  in  THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering)   The University of British Columbia (Vancouver)  April 2015  © Pooyan Kheirkhah, 2015  ii  Abstract Direct injection of natural gas in compression-ignition engines offers benefits such as emissions reduction, fuel diversity, and energy security. However, in order to meet the upcoming stringent emissions regulations, further improvements in the performance of the High-Pressure Direct-Injection of Natural Gas (HPDI-NG) is needed. For this reason, different natural gas injection strategies and nozzle designs are numerically studied. The in-cylinder phenomena during the closed portion of the cycle is simulated using a Large-Eddy Simulation (LES) turbulence model and the Trajectory-Generated Low Dimensional Manifold (TGLDM) for chemistry coupling. Soot is modeled with a two-equation Hiroyasu model. To partially investigate the effect of LES variability, simulations with successive mesh refinements and infinitesimally varied inputs are carried out. Anticipated emission trends were observed for parametric sweeps with substantial variation of soot about the trend-line. This motivated the analysis of in-cylinder mixture, jet penetration and other more robust metrics. A novel paired-nozzle geometry was designed to increase the fuel-air mixing at the base of the jet, thus reducing soot. In reality, the paired jets increased exhaust PM. The CFD analysis revealed that the gas jet penetration was reduced compared to the baseline single-hole jet, while more air was entrained into the core of the jet. However, the effect of mixing due to impaired penetration dominates and results in more rich mixture and therefore more soot. CFD predicted the PM reduction benefits of “Late Post-Injection” (LPI) due to two major reasons: 1- the reduction in formed PM from the 1st pulse due to shortened pulse width, 2- negligible PM formation from the 2nd pulse for enough pulse separation. A second injection strategy, “Slightly Premixed Combustion” (SPC) also reduced PM in experiments. The CFD package had not been developed for such combustion regime, wherein the diesel-gas kinetic interactions should be resolved; hence perfect matching between the experimental and numerical combustion for SPC was not attained. Nevertheless, by iii  optimizing the injection timing to resemble the phasing of experimental Heat Release Rate (HRR) curve, to the “best extent possible”, more premixing, higher rate of penetration, and less rich-mixture mass was observed.   iv  Preface The author of this thesis conducted the simulations and subsequent data processing. The CFD tool used for this study, GOLD, has been developed at Westport Innovations Inc. by Dr. Jim Huang. At the author’s request, some adjustments to the source code was made by Jason Zheng at Westport Innovations Inc. to reproduce injection behaviour similar to the experimental data for different injection strategies. The paired-hole nozzles for HPDI operation was originally designed by Ehsan Faghani and Dr. Steven Rogak. The engine tests were conducted by the graduate students, Ehsan Faghani and Chris Mabson at UBC. The author wish to thank Ehsan and Chris for supporting this study by providing necessary experimental data and invaluable insight.   v  Contents Abstract ............................................................................................................................................ ii Preface ............................................................................................................................................. iv Contents ........................................................................................................................................... v List of Tables .................................................................................................................................. viii List of Figures ................................................................................................................................... ix Nomenclature .................................................................................................................................. xi Acknowledgements ........................................................................................................................ xv 1 Introduction ............................................................................................................................. 1 1.1 Engine emissions ............................................................................................................. 1 1.1.1 Environmental and health impacts ......................................................................... 1 1.1.2 Emission regulations and test cycles ....................................................................... 2 1.2 Compression ignition of natural gas ................................................................................ 3 1.2.1 Natural gas as a fuel ................................................................................................ 3 1.2.2 Direct injection of natural gas – HPDI technology ................................................... 4 1.3 Formation and control of emissions ................................................................................ 5 1.3.1 Pollutant formation in direct-injection CI engines .................................................. 5 1.3.2 Exhaust Gas Recirculation (EGR) ............................................................................. 6 1.3.3 Injection pressure .................................................................................................... 6 1.3.4 Group-hole nozzles .................................................................................................. 6 1.3.5 Multiple injections ................................................................................................... 7 1.3.6 Low-temperature combustion ................................................................................ 7 1.4 Application of CFD for studying IC engines ..................................................................... 8 1.5 Motivation, objectives, and scope .................................................................................. 9 1.5.1 Motivation ............................................................................................................... 9 1.5.2 Objectives and scope ............................................................................................. 10 2 Research methods and post-processing techniques ............................................................. 12 2.1 Engine simulation with OpenFOAM .............................................................................. 12 2.1.1 CFD grid and numerical scheme ............................................................................ 12 2.1.2 Turbulence modeling ............................................................................................. 14 2.1.3 Wall treatment ...................................................................................................... 18 vi  2.1.4 Diesel injection ...................................................................................................... 18 2.1.5 Natural gas injection .............................................................................................. 19 2.2 Modeling combustion chemistry in GOLD ..................................................................... 21 2.2.1 Natural gas combustion mechanism ..................................................................... 21 2.2.2 Diesel combustion mechanism .............................................................................. 22 2.2.3 Simplification of reaction mechanisms ................................................................. 22 2.2.4 Soot chemistry and mathematical modeling ........................................................ 26 2.3 The Single-Cylinder Research Engine (SCRE) ................................................................. 29 2.4 Calculation of CFD inputs .............................................................................................. 32 2.5 CFD post-processing methods ....................................................................................... 33 2.5.1 Heat release rates .................................................................................................. 33 2.5.2 Equivalence ratio – temperature maps ................................................................. 34 2.5.3 Equivalence ratio and mixture fraction distribution ............................................. 35 2.5.4 Rich and lean mixture evolution, air entrainment ................................................ 36 2.6 Numerical verification and sensitivity study ................................................................. 37 3 Results ................................................................................................................................... 39 3.1 GOLD validation ............................................................................................................. 39 3.1.1 Pressure, heat release rate, and emissions ........................................................... 39 3.1.2 GOLD sensitivity to mesh ....................................................................................... 45 3.1.3 GOLD sensitivity to input parameters and discussion of cyclic variability ............ 51 3.2 Analysis of the paired-nozzle injectors .......................................................................... 55 3.2.1 Heat release rates and final emissions at mode B-75 ........................................... 55 3.2.2 Penetration, air entrainment, and mixture formation .......................................... 59 3.2.3 Equivalence ratio - temperature maps .................................................................. 62 3.2.4 Mixture probability density function ..................................................................... 64 3.2.5 Jet interaction with cylinder and fire-deck ............................................................ 66 3.3 Analysis of the LPI strategy ............................................................................................ 68 3.3.1 Heat release rates and final emissions .................................................................. 68 3.3.2 PM formation in LPI strategy ................................................................................. 71 3.4 Analysis of SPC injection strategy .................................................................................. 73 3.4.1 Heat release rates and final emissions .................................................................. 73 3.4.2 Mixture formation ................................................................................................. 75 vii  4 Conclusions ............................................................................................................................ 80 5 Recommendations ................................................................................................................. 82 Bibliography ................................................................................................................................... 83 Appendices .................................................................................................................................... 90 Appendix A: Initial engine mixture calculation.......................................................................... 90 Appendix B: Φ-T map calculations ............................................................................................ 92    viii  List of Tables  Table 1- European emission regulations ......................................................................................... 2 Table 2- CFD grid information for J-36 and 7-7-0 injectors ........................................................... 13 Table 3- Charge composition for the tabulated manifolds ........................................................... 25 Table 4- B-75 Point over time coefficient of variance ................................................................... 30 Table 5- Design specifications of the simulated nozzles ............................................................... 32 Table 6- Initial charge condition and injection specifications for GOLD ....................................... 32 Table 7- Operating parameters at ESC mode B75 (M4) ................................................................ 39 Table 8- CFD inputs for baseline simulations ................................................................................ 40 Table 9- SCRE commanded inputs for EGR sweep test ................................................................. 42 Table 10- CFD inputs for EGR sweep simulations .......................................................................... 42 Table 11. Mean and standard error for the penetration of the jets ............................................. 50 Table 12- SCRE commanded parameters for operation under mode B-75. ................................. 55 Table 13- Inputs for CFD simulation of paired-nozzle injectors .................................................... 56 Table 14- The gap between the gas nozzle centreline and the cylinder head .............................. 66 Table 15- Performance of the optimized LPI strategy ................................................................... 68 Table 16- SCRE and GOLD injection specifications for LPI tests. ................................................... 69 Table 17- SCRE operating parameters in SPC strategy .................................................................. 74 Table 18- Sample initial charge composition for different working conditions ............................ 91   ix  List of Figures  Figure 1- U.S petroleum production and consumption and fossile fuel CO2 emission Statistics .... 1 Figure 2- Operating modes in European Stationary Cycle tests. ..................................................... 3 Figure 3- Westport HPDI-2 injector ................................................................................................. 5 Figure 4- CFD grid structure from top and side views ................................................................... 14 Figure 5- Pressure fluctuation in the nozzle near-field region ...................................................... 20 Figure 6- Reaction pathways for methane oxidation at 40 bar  .................................................... 22 Figure 7- physical and chemical time scales in a turbulent flame ................................................. 23 Figure 8- Low dimensional manifolds for a CO-H2-air system ...................................................... 24 Figure 9- schematic diagram of the different blocks in GOLD software ....................................... 26 Figure 10- Phenomenological picture of soot formation .............................................................. 27 Figure 11- Arrangement of holes in the baseline 7-7-0  and paired-hole injectors. ..................... 31 Figure 12- Equivalence ratio-temperature maps generated using CANTERA.. ............................. 35 Figure 13- PDF of equivalence ratio for the baseline and negative PSEP injection. ..................... 36 Figure 14- Pressure and rate of heat release predictions at the baseline mode .......................... 40 Figure 15- Effect of injection timing on HRR curves of 7-7-0 and J-36 injectors ........................... 41 Figure 16- Emission trends during an EGR sweep test .................................................................. 43 Figure 17- PM-NOx and PM-CO maps for single-nozzle injectors ................................................. 45 Figure 18- Sensitivity of CFD emissions to the mesh circumferential resolution .......................... 47 Figure 19- Gas injection in HPDI engine ........................................................................................ 48 Figure 20- Penetration sensitivity to the mesh resolution ............................................................ 49 Figure 21- Comparison between 7-7-0 jets with coarse and fine meshes .................................... 50 Figure 22- Effect of mesh resolution on mixture production ........................................................ 51 Figure 23- Effect of Injection pressure on emissions .................................................................... 53 Figure 24- Effect of combustion phasing and injection timing on emissions ................................ 54 Figure 25- The rate of heat release curves for HPDI paired-nozzle injectors ................................ 56 Figure 26- Misalignment between the gas and diesel holes in L-H-S-A injector ........................... 57 Figure 27- HRR graphs of repeated simulations with adjusted interlace angle. ........................... 58 Figure 28- Measured and simulated emissions for 7-7-0 and paired-nozzle injectors ................. 58 Figure 29- Paired-nozzle emissions added to the baseline emission maps .................................. 59 Figure 30- Penetration and entrainment plots for baseline and paire-nozzle injectors. .............. 60 Figure 31- Evolution of soot and rich and lean mixtures in the cylinder ...................................... 61 Figure 32- Status of the mixture on the φ-T planes for the 7-7-0 injector ................................... 63 Figure 33- Plots of fuel mass difference between baseline and paired-nozzle injectors .............. 64 Figure 34- The distribution function of the equivalence ratio ...................................................... 65 Figure 35- Fuel jet interaction with the cylinder wall ................................................................... 67 Figure 36- Penetration and rich-mixture evolution for injectors of small and large protrusion. . 67 Figure 37- Heat release rates for LPI simulations .......................................................................... 69 Figure 38- Emissions for single and late-post injections. .............................................................. 70 x  Figure 39- PM-NOx and PM-CO emission maps............................................................................. 71 Figure 40- Evolution of in-cylinder soot for single and post injections ......................................... 72 Figure 41- Regions of high PM concentration in a close-coupled versus split-flame regimes. ..... 73 Figure 42- Simulated and SCRE HRR curves for the standard J36 actuation delay ....................... 75 Figure 43-Simulation of SCRE at SPC regime with varied PSOI ..................................................... 75 Figure 44- Evolution of the rich and lean mixtures for the simulated SPC cases .......................... 76 Figure 45- fuel penetration under baseline and SPC J36 injection ............................................... 77 Figure 46- Mixture difference between positive and negative PSEPs on φ-T maps. .................... 78 Figure 47- PM evolution for J-36 at PSEP=0.3 ms, PSEP=-0.9 ms, and PSEP=-2.1 ms ................... 79 Figure 48- Flowchart for cylinder (EGR-diluted) charge composition ........................................... 90 Figure 49- Flowchart for binning the mixture on a φ-T plane. ...................................................... 92   xi  Nomenclature Abbreviations 7-7-0  HPDI injector with 7 gas and 7 aligned diesel holes aSOI  after Start of Injection aTDC  after Top Dead Centre CFD  Computational Fluid Dynamics CI  Compression Ignition CNG  Compressed Natural Gas COV  Coefficient of Variance CSE  Conditional Source-term Estimation DNS  Direct Numerical Simulation EGR  Exhaust Gas Recirculation EQR  Equivalence Ratio ESC  European Stationary Cycle EVO  Exhaust Valve Opening GIMEP  Gross Indicated Mean Effective Pressure GOLD  Westport in-house CFD package for in-cylinder engine simulation GRI  Gas Research Institute GRP  Gas Rail Pressure GSEP  Gas pulse separation GSOI  Gas Start of Injection HACA  Hydrogen Abstraction C2H2 Addition HC  Hydrocarbon emissions HCCI  Homogenous-Charge Compression-Ignition HPDI  High-Pressure Direct-Injection HRR  Heat Release Rate ICE  Internal Combustion Engines IHR  Integrated Heat Release ILDM  Intrinsically Low Dimensional Manifolds IQ  Index of Quality J36  HPDI injector with 9 gas holes and 7 diesel holes xii  KH-RT  Kelvin-Helmholtz Rayleigh-Taylor LES  Large-Eddy Simulation L-H-S-A  Large-Hole Small-Angle LNG  Liquefied Natural Gas LPI  Late Post-Injection LPM  Littre per Minute NOx  Oxides of Nitrogen NTC  Negative Temperature Coefficient PAH  Poly-Cyclic Aromatic Hydrocarbons PDF  Probability Density Function PM  Particulate Matter PPCI  Partially-Premixed Compression Ignition PPW  Pilot Pulse Width PRF  Primary Reference Fuel PSEP  Pulse Separation PSOI  Pilot Start of Injection RANS  Reynolds-Averaged Navier-Stokes SCRE  Single-Cylinder Research Engine SGS  Sub-Grid Scale S-H-S-A  Small-Hole Small-Angle SI  Spark Ignition SMPS  Scanning Mobility Particle Spectrometer SOI  Start of Injection SPC  Slightly Premixed Combustion TEM  Transmission Electron Microscopy TEOM  Tapered Element Oscillating Microbalance TGLDM  Trajectory-Generated Low Dimensional Manifold TKE  Turbulent Kinetic Energy uHC  unburned hydrocarbons ULSD  Ultra Low-Sulfate Diesel   xiii  Symbols ( ̃ )   Favre-averaged/filtered quantity ( ̅ )   Ensemble averaged - or spatially filtered quantity ( )𝑠𝑓   Related to soot formation ( )𝑠𝑜   Related to soot oxidation ( )𝑠   Related to net soot ( )𝛼   Related to species α ( )′′   Reynolds or sub-grid scale turbulent stresses ( )𝑡   Turbulent quantity ( )𝑓  Forward ( )𝑟  Reverse Λ𝐾𝐻   Kelvin-Helmholtz instability wave number 𝐵0   Kelvin-Helmholtz instability constant 𝐵1   Kelvin-Helmholtz break-up constant 𝑐  Molar concentration 𝐶Δ   Coefficient of eddy viscosity in LES 𝐶𝜇   Eddy viscosity coefficient in k-ɛ model 𝐶𝜏   Rayleigh-Taylor break-up time constant 𝐸𝑎   Reaction Activation Energy 𝑛𝑠   Number of species in a reacting system ℎ   Enthalpy Δ   Filter width in LES simulation 𝑀   Molar mass or injected fuel mass 𝑁𝑢   Nusselt number 𝑃𝑟   Prandtl number 𝑅  Gas constant 𝑅𝑒   Reynolds number 𝑆  System change due to chemical reactions – or strain rate tensor 𝑇   Temperature 𝑈   Velocity 𝑍   Mixture fraction 𝑗   Diffusive flux 𝑘   Turbulent Kinetic Energy 𝑝   Pressure 𝑦   Mass fraction of species 𝜁   Mixture Fraction 𝜃   Crank Angle 𝛾   Specific heat ratio 𝜈   Kinematic viscosity xiv  𝜌   Density 𝜏   Ignition/Injection delay, Stress tensor, or break-up time 𝜖   Turbulent kinetic energy dissipation rate 𝜙   Equivalence Ratio   xv  Acknowledgements I am delighted to take this opportunity to thank those who helped me to reach this far in my education. I was given the chance to study and research at UBC thanks to my supervisor, Dr. Steven Rogak; this was a true privilege and I am very thankful to him for this opportunity. I truly benefitted from his teachings and research expertise. I would also like to extend my gratitude to Dr. P. Kirchen who provided invaluable feedback to this research; I owe him many skills in research and teaching. I am especially thankful to Dr. J. Huang for his technical support and contribution to this work; without access to the resources he provided, this work was impossible. Through this two and a half years, I have had the benefit of working with outstanding researchers and engineers. This research was enriched through many critical evaluations from Dr. Gordon McTaggart-Cowan, Dr. Ning Wu, and Mr. Bronson Patychuck. I learned discipline, focus, and dedication in addition to invaluable technical knowledge from them. Many thanks goes to my fellow researchers: Ehsan Faghani and Chris Mabson, who provided experimental inputs for my work. Their expertise at engine experiments and modeling helped me to focus my efforts in the right track. Ehsan and Chris thank you for many useful discussions over coffees and lunches and practical lessons in the engine cell and skiing sessions! I had the pleasure of working and befriending with Amin Engarnevis, James Montgomery, and Ramin Dastanpour during my studies at this research group. Amin, the laughter, bus rides, beers, and discussions we shared was a mental fuel giving me energy and motivation to continue. I am truly delighted to have you as a friend. I am thoroughly thankful to my parents; they were always supportive of my educational pursuits. The hospitality and kindness of my aunt, uncle, and cousin have reminded me there is more to life than research. My brother, Peyman, thanks for your great sense of humor; nothing cheers me up more than talking to you. The presence of all of these xvi  people made me a better engineer, researcher, and a most important of all a better human.1  1 Introduction  1.1 Engine emissions 1.1.1 Environmental and health impacts Internal combustion engines (ICE) have applications in industry, power generation, and transportation. The petroleum consumption in transportation sector alone equals the net import of petroleum in the U.S. (Figure 1), and in Canada 30% of fossil fuel consumption is attributed to transportation [1–3]. ICE-emitted CO2 emissions in industry and transportation had fastest growth since 1970 (Figure 1) [4].     By-products of diesel combustion affect air quality and accelerate global warming [5]. The carbonaceous particulate emissions, soot, are found to have pathological effects such as cough, lung function deterioration and nausea [6] as well as carcinogenic and mutagenic effects [7,8]. As a result of such mechanistic evidence diesel exhaust is categorized as group A (carcinogenic to humans) in 2012 [9]. Moreover, according to atmospheric studies NOx, another major diesel emission, have devastating environmental and health effects, including ozone-generating photochemical smog, acidic rain, and respiratory system irritation [10]. Figure 1- U.S petroleum production and consumption statistics [1] (left), and global/U.S fossil fuel CO2 emissions [4] (right) 2  1.1.2 Emission regulations and test cycles Heavy-duty engine emissions standards are classified by the size and weight of the vehicles, load and power ratings, and their application. In this study, emissions of a heavy-duty engine meant for on-road transportation are investigated and therefore relevant standards are briefly explained here. In Europe, vehicles with the maximum laden mass of more than 3500 kg are classified as heavy-duty and are subject to emission limits outlined in Table 1. Notably, the PM and NOx mass have been dramatically reduced in the most recent regulations and the particle numbers are added to the Euro IV standard [11]. Table 1- European emission regulations [11]  Date CO HC NOx PM PN g/kWh 1/kWh Euro I 1992 ≤ 85 kW 4.5 1.1 8.0 0.612  1992 ≥85 kW 4.5 1.1 8.0 0.36  Euro II 1996.10 4.0 1.1 7.0 0.25  1998.10 4.0 1.1 7.0 0.15  Euro III 1999.10 EEV only 1.5 0.25 2.0 0.02  2000.10 2.1 0.66 5.0 0.10  Euro IV 2005.10 1.5 0.46 3.5 0.02  Euro VI 2013.01 1.5 0.13 0.40 0.01 8.1x1011  To certify with European standards, engines are run on the European Stationary Cycle (ESC) or European Transient Cycle (ETC). In ESC tests, the engine emissions and fuel consumption are measured in 13 operating modes including idling and factorials of engine speed and load. The net result is brake-specific emission and fuel consumption averaged over all the cycle tests. Figure 2 shows the map of ESC testing conditions with the averaging factors at operating condition [12]. The engine tests modeled in this thesis are similar to some of the ESC modes, Mode-4 or B-75 is especially important here because of its relatively high weighting factor and high PM emissions. 3      1.2 Compression ignition of natural gas 1.2.1 Natural gas as a fuel Natural gas is cheaper than conventional diesel and it has less tendency to form soot and NOx per unit output work due to less carbon-carbon bonds and lower flame temperature. These features have made natural gas an attractive fuel in transportation sector, especially given the introduction of the recent emission regulations worldwide [13]. Even though its lower energy density makes onboard natural gas storage a challenge, many heavy-duty trucks for the transportation of goods and medium to heavy-duty buses have been equipped with LNG and CNG tanks. Despite cleaner operation, natural gas spark ignition engines suffer from impaired brake power output. Natural gas port-injection deteriorates the ingestion volumetric efficiency due to the replacement of fresh air with fuel. Consequently, direct injection of natural gas with the aid of an ignition source such as a spark plug, hot surface glow plug, or pilot fuels have been developed [14]. Figure 2- Operating modes in European Stationary Cycle tests [12]. 4  1.2.2 Direct injection of natural gas – HPDI technology Compression ignition engines are inherently more efficient due mostly to higher compression ratios and throttle-less operation. Indeed, diesel engines are the engines of choice for heavy duty trucks due to higher efficiency, desirable load characteristic, durability, and reliability. Also, owing to the recent advances in the fuel injection systems, the diesel engines have found applications in light-duty vehicles as well. In spite of better fuel economy, diesel engines are prone to high levels of PM and NOx emissions due to the non-premixed nature of combustion in these engines [14,15]. This has made the introduction of natural gas in compression ignition engines an attractive theme, with many technologies developed around the idea of  mixed-fuel blends [16,17]. To achieve the full benefits of compression ignition, the natural gas should be directly injected in the cylinder near the top-dead-center. However, given its very low cetane number (equivalently, its long ignition delay), a separate source of ignition is required. A very reliable source of ignition is a high-cetane pilot fuel injected just before the injection of natural gas. Westport Innovations has developed this idea with the development of dual fuel injectors for diesel engines. Both gas and diesel are injected using the same injector, hence minimum or no modification to the engine head is required. These HPDI injectors consist of two concentric needles, a large gas needle encompassing the thinner diesel needle inside [18]. The diesel is injected from small holes at the tip of the injector followed by injection of gas from separate nozzles located higher in the injector body. Figure 3 shows the schematic diagram of Westport HPDI injector.  5     1.3 Formation and control of emissions The main regulated emissions from compression ignition engines include NOx, CO, methane and non-methane HC, and PM.  Their formation and control is briefly reviewed below.  After-treatment systems are used for nearly all heavy-duty engines today, but their cost and performance penalties [19,20] motivate the reduction of pollutants from the engine as much as possible.  1.3.1 Pollutant formation in direct-injection CI engines Different emissions are formed in different regions of the flame. The main route to nitrogen oxide formation includes high activation energy reactions; therefore, its peak concentration is observed around the stoichiometric, high-temperature envelope [21]. Carbon monoxide is mainly formed in oxygen-deficient regions; the majority of formed CO is oxidized later as more air is entrained into the jet; although, minor amounts survive this oxidation process [22]. Unburned hydrocarbons (UHC) are emitted due to mixture packets escaping flame surface due to turbulence intermittency, shear flame extinction, wall quenching, and bulk quenching at expansion stroke [21]. Particulate emissions consist mainly of fractal-like agglomerated carbonaceous particles coated with semi-volatile organic material. Traces of metallic particles, ash, and sulfates can be internally Figure 3- Westport HPDI-2 injector [18] 6  or externally mixed with the particle emissions. Typically, the soot formation has highest intensities in rich parts of the flame where abundant amounts of carbon molecules exist and the temperature is high enough to favour surface reactions leading to soot mass growth. A large portion (>90%) of the formed soot is later oxidized when mixing with more oxygen downstream of the jet [23]. 1.3.2 Exhaust Gas Recirculation (EGR) EGR systems return 10-20% of the exhaust gases to the intake manifold. Use of recirculated gases lowers the peak combustion temperature, thus greatly reducing NOx, yet increasing soot emissions due to reduced O2 concentration. Such behaviour in diesel engines is known as the “PM-NOx” trade-off. Because of this trade-off, it is especially important to use PM reduction strategies in engines using EGR. Different PM-NOx trade-offs are also created by sweeping parameters such as injection timing or boost pressure. 1.3.3 Injection pressure The average injection pressure in DI diesel engines has increased from 10-20 MPa to over 200 MPa over the past twenty years. Better atomization and fuel vaporization and higher rate of penetration are introduced as the effective soot reduction mechanism in high-pressure diesel injection [24]. Previous research on HPDI engines also showed improvement in the PM emissions with higher injection pressures, yet NO2  emissions were slightly increased as a result of more premixing and slightly advanced combustion [13]. 1.3.4 Group-hole nozzles  Grouping the injection holes into pairs instead of multi single-hole nozzles has been practiced for diesel injection. The effects of paired-nozzle included angle and spacing have been investigated. It has been consistently observed that jet penetration decreases for paired-nozzles; more reduction was reported for larger included angle and nozzle spacing. Nevertheless, as a result of better atomization, better fuel evaporation and mixture preparation was observed compared to a single-hole nozzle [25,26]. Engine optical studies have shown less flame luminosity when grouped nozzles are used, concluding improved atomization and air entrainment improved the mixture preparation and 7  reduced soot [27]. Other numerical studies found that the paired-nozzles do entrain more air if the injection has not interrupted; however, upon the termination of injection, the rate of air entrainment is much reduced compared to single-hole baseline nozzle [28] Overall, the jet penetration and atomization are found to have a trade-off behaviour; since the two are important mechanisms for mixture preparation the sensitive dependence of pollutant formation on hole configuration is emphasized [27,28]. 1.3.5 Multiple injections Multiple injection strategies are practiced in diesel engines for several reasons; pilot injections reduce the combustion noise and assist the overall NOx reduction [29]. The post injection are utilized for increasing exhaust temperature for the regeneration of diesel catalysts and filters, oxidizing the in-cylinder UHC and CO, and reducing particulate emissions. Some studies proposed increased mixing and higher temperatures, caused by the high-speed turbulent injection and combustion of the second pulse, as the effective soot reduction mechanism. Nevertheless, separating the effect of mixing and temperature on the performance of post injection is a formidable research task and no conclusive answer has been offered in the literature of post injection [30]. Alternatively, the effect of injection duration and fuel quantity in each pulse is proposed as the main soot limiting quantity for split injections. Reducing the soot replenishment in the head of the jet by interrupting the injection is deemed to be the main mechanism in this technique [31]. The literature does not provide clear guidance on whether or not post injection would benefit HPDI engines. Recent measurements at UBC demonstrated major PM reductions using late post-injection, without negatively affecting other emissions and only 1-2 percent fuel consumption penalty [32]. However, no explanation of the underlying soot reduction mechanism has been offered.  1.3.6 Low-temperature combustion Ideas have been developed to avoid rich regions encountered in non-premixed diesel combustion. In Homogenous Charge Compression Ignition (HCCI) a lean mixture is produced by very early injection of the fuel, which unlike SI or CI engines burns rapidly in a volumetric fashion 8  around TDC. Soot and NOx emissions are greatly reduced but high CO and UHC emissions are encountered [33]; cold starts and transient operation are difficult [34]. By allowing some degree of inhomogeneity by later fuel injection, a Partially Premixed combustion (PPCI) is realized which improves some of the mentioned deficiencies of HCCI operation [35,36]. HPDI offers independent control over the ignition of the main fuel through the use of diesel pilot; this feature has been exploited to overcome combustion control barriers in PPCI-like regimes for HPDI. Partially premixed natural gas combustion by injecting the gaseous fuel during the intake stroke to allow full premixing, then igniting it by small diesel injection was tested for low-moderate loads. To extend its applicability to full loads, a second quantity of natural gas was injected close to the TDC to provide the extra load [37]. An improved version of this strategy is developed by slightly advancing the gas injection relative to the pilot to avoid highly premixed knock issues. McTaggart-Cowan et al. found that implementing this negative PSEP strategy results in more rapid, intense combustion with significant reductions of PM at the expense of increased NOx and hydrocarbon emissions. Conventional PM-NOx trade-off could be beaten in this Slightly Pre-mixed Combustion (SPC) mode: by increasing EGR levels, NOx emissions are reduced to below the baseline levels while PM emission are unaffected [38]. 1.4 Application of CFD for studying IC engines The application of computational tools in engineering has increased due to increased computer power and better software [39,40]. In particular, CFD analysis of IC engines is becoming useful for research and development. GM reported the use of in-house CFD for the optimization of a direct-injection 3.6 L DOHC engine through adjustment of injection parameters and better understanding of mixture-chamber interaction [41]. The Engine Research Centre (ERC) in the University of Wisconsin Madison reported the use of a 1-dimensional gas dynamics code to model the charge induction process linked with a 3D full-cycle simulation using KIVA-3V for the optimization of valve lift profiles [42]. 9  At a more fundamental level, CFD has been used to understand the charge motion and different injection and combustion strategies in IC engines. Proper data reduction has been an essential part of this; especially in analysis of new combustion regimes, reduced data using   probability-like distribution of mixture and equivalence ratio-temperature (φ-T) maps have facilitated better understanding [25,43,44]. Turbulent Large Eddy Simulation (LES) studies have gained much popularity in recent years owing to the improved computer power, and the development of LES sub-models [45]. LES is particularly helpful in studying the inherently transient engine phenomena that were not reliably available with ensemble or Reynolds-averaged models. The chaotic nature of cylinder flow is resolved to the mesh-size level; therefore the issue of cyclic variability could be studied in more depth. Numerical studies of multiple consecutive motored cycles showed considerable charge motion fluctuation among cycles, with as much as 25 consecutive cycles needed for velocity averaging [46]. 1.5 Motivation, objectives, and scope 1.5.1 Motivation To reduce the in-cylinder emissions of a heavy-duty Cummins ISX 6-cylinder HPDI engine, several injection strategies were tested on a Single-Cylinder Research Engine (SCRE) [47–49]. The overall test matrices were focused on the following injector nozzle designs and injection strategies: 1- Paired-hole injector nozzle design: The concept of paired-hole injection of gas was developed by Faghani and Rogak [50], in which the surface area of a single nozzle is divided into two closely-spaced nozzles. The main concept was formed based on the theory of gaseous jets; upon the reduction of gas nozzle diameter, the rate of entrainment into the jet is increased in the cost of slower penetration. To counter this, the two pairs are closely spaced and the surface area of each is increased in one sample of the paired-nozzle combination to maintain the same penetration. In reality, the particulate emissions were higher in all modes with all the paired-hole nozzles.  2- Late Post Injection (LPI): a small pulse of natural gas (15% of overall cycle energy) is injected late in the cycle, to meet the cycle load and reduce the emission of PM. In the 10  consequent tests promising results were obtained; PM emissions were reduced by 75% while other emissions were reduced or remained close the baseline values. 3- Slightly Premixed Combustion (SPC): This strategy has previously been tested both on a single-cylinder and on a multi-cylinder production engine with promising results. Yet, the high rates of NOx and HC emissions were reported; thus in more recent tests by Ehsan Faghani these issues were tackled. It was found that the PM emissions were far less sensitive to EGR and EQR level as opposed to the conventional operating modes. The phenomena described above cannot be understood without knowing something about in-cylinder flows and concentrations – information that can be obtained (at least approximately) from CFD modelling. 1.5.2 Objectives and scope The objective of this work is to numerically simulate the mentioned injection strategies and nozzle designs, and to develop necessary post-processing techniques for answering the following questions: 1- Do paired-hole nozzles generate more PM because of undesired mixing emanated from paired-jet dynamics or unforeseen interactions with cylinder enclosure cause this behaviour. 2- Are the reductions in LPI strategy due to late stage PM oxidation from the second pulse, less PM formation in the first pulse, or the enhanced entrainment effect? 3- How is the effect of SPC mixing different than other modes of mixing in conventional HPDI combustion? Why are the results unaffected by EGR and EQR levels? These questions are tackled numerically, using a 3D CFD software developed on OpenFOAM solver code, which is modified and coupled with an in-house chemistry model for natural gas developed. The solver itself was developed and validated against experimental data at Westport Innovations Inc.; some case-specific adjustments are made for each injection strategy and several post-processing methods are employed for better understanding of the mixing process. In chapter 2, the background information related to the CFD tool is presented; then, the details of the numerical tool and simulation conditions are described in that chapter. In chapter 3, the 11  results of the simulated injection strategies are analysed using different methods such as φ-T maps, distribution plots and mixture evolution. Finally, major findings of this study are summarized in the last chapter with recommendations for future research.   12  2 Research methods and post-processing techniques This section describes the engine CFD simulation methods, single-cylinder engine research facility, the methods for post-processing and comparing the measurements and CFD predictions. Physical features including turbulent charge motion, chemical reactions, spray atomization, and pollutant formation need to be modeled in CFD simulation of ICEs. OpenFOAM is an open-source package developed based on C++ object-oriented language. It offers great flexibility in the choice of CFD mesh even for complicated 3D applications, robust numerical schemes, the ability to run parallel computations, and a rich variety of physical models applicable to IC engines [51]. OpenFOAM is chosen as the base CFD platform for this study given its flexibility for moving unstructured mesh, and its convenient engine-related libraries. In explaining the methods and tools, a brief review of background information is provided as necessary. 2.1 Engine simulation with OpenFOAM This thesis used GOLD, developed at Westport Innovations Inc. by Jim Huang [52]. GOLD is built on OpenFOAM (see section 1.4), with additional modules to simulate the combustion of diesel and natural gas in a non-premixed turbulent regime, described below. 2.1.1 CFD grid and numerical scheme  To avoid excessive computational costs, GOLD simulates a one-sector region of the combustion chamber including one1 gas nozzle and one diesel nozzle; hence periodic (cyclic) boundary conditions are chosen for the sector partitions (left and right boundaries). Moreover, instead of simulating an entire 4-stroke engine cycle, i.e. 720 degrees of crank shaft rotation, only the closed portion of the cycle starting from -90° aTDC continuing until EVO at 140°. Therefore, the extra burden of simulating intake and exhaust flows in the ports and around the valves are avoided, also, the simulation duration is reduced to only 230° compared to a full 720° cycle, i.e. a 70% reduction. For further efficiency, the domain is re-meshed near TDC such that extreme cell aspect ratios are avoided and highly-resolved CFD grid during gas injection is obtained. Throughout the cycle, the ratio of grid size between different regions of the combustion chamber is kept constant.                                                           1 It could be one-pair, when simulating paired-hole injectors. 13  In particular, the number of circumferential grids on the gas nozzle is a key parameter defining the mesh resolution and is refined during injection. The number and types of the CFD cells for the simulation of the two baseline injectors for different crank angle intervals are summarized in Table 2. Table 2- Number of CFD grids and circumferential resolution for J-36 and7-7-0 injectors  7-7-01 J362 # of cells dθnozzle [deg] dθchamber [deg] # of cells dθnozzle [deg] dθchamber [deg] -90° ~ -70° aTDC 127,374 1.55 2.1 124,894 1.35 1.92 -70° ~ -20° aTDC 97,214 95,414 -20° ~ -10° aTDC 51,922 51,172 -10° ~ 10° aTDC 87,400 1.03 1.4 93,262 0.9 1.24 10° ~ 20° aTDC 47,398 1.55 2.1 46,750 1.35 1.92 20° ~ 30° aTDC 65,494 64,438 30° ~ 40° aTDC 65,654 93,918 40° ~ 140° aTDC 128,882 126,368  Cylinder mesh refinement occurs between -10° to 10° aTDC as shown in Figure 4; between re-meshing intervals, the grid is compressed or expanded to match the cylinder volume, but no circumferential or radial refinement occurs. OpenFOAM’s mesh motion solver utilities are used to calculate the mesh motion in re-meshing intervals. A first-order-accurate scheme for time integration and a second-order-accurate upwind scheme for spatial integration are employed. A preconditioned conjugate gradient solver for symmetric matrices (PCG) is employed to solve the linear system of equations for pressure and density; for the rest of the parameters a preconditioned bi-conjugate gradient solver for asymmetric matrices is selected.                                                            1 Baseline injector for paired-hole nozzle study. This injector has 7 gas holes and 7 aligned diesel holes. 2 HPDI production injector with 9 gas and 7 diesel holes. In this thesis a 9-9 configuration is simulated to be able to use a sector mesh. 14        2.1.2 Turbulence modeling To fully resolve the turbulent flow, numerical girds of the size of the smallest eddies (Kolmogorov length) is needed [53]. To obtain a Direct Numerical Simulation (DNS) solution, the number of CFD cells scales with the 3rd power of Re number; thus requiring immense computer power. For practical flow problems, this is prohibitive, too costly, and often undesired due to the huge amount of data not necessarily needed for an engineering analysis. To avoid this, some modeling is needed in turbulent flow problems. In the Large-Eddy Simulation (LES) approach, the flow equations are spatially filtered; resulting in quantities resolved on scales larger than the filter size, while the effect of smaller eddies is modeled. When modeling turbulent reactive flows, the density variation is significant and cannot be ignored; also, direct averaging or filtering results in higher moments of density that complicates the analysis. To minimize these issues, density-weighted averages (Favre averaged) are used: Figure 4- CFD grid structure from top and side views. The mesh is refined during injection, i.e. in -10° to 10° aTDC (right). 15  ?̃? =𝜌𝑞̅̅ ̅̅?̅?  1 – a 𝑞′′ = 𝑞 − ?̃?  1 – b After Favre-averaging, the Navier-Stokes along with energy and species transport equations will read as follows: 𝜕?̅?𝜕𝑡+𝜕?̅?𝑢𝑖𝜕𝑥𝑖= 𝜌?̇?𝑙𝑖𝑞  Mass 2 – a 𝜕(?̅?𝑢𝑖)𝜕𝑡+𝜕(?̅?𝑢𝑖𝑢𝑗)𝜕𝑥𝑗= −𝜕?̃?𝜕𝑥𝑖+ 𝜕?̃?𝑖𝑗𝜕𝑥𝑗+ 𝜌?̇?𝑗 𝑢,𝑙𝑖𝑞 −𝜕(?̅?𝑢𝑖′′𝑢𝑗′′̃ )𝜕𝑥𝑗  Momentum 2 – b 𝜕(?̅??̃?𝛼)𝜕𝑡+𝜕(?̅?𝑢𝑗?̃?𝛼)𝜕𝑥𝑗= −𝜕𝐽𝑗𝛼𝜕𝑥𝑗+𝑊𝛼?̅̇?𝛼 + 𝜌?̇?𝛼 𝑙𝑖𝑞 −𝜕(?̅?𝑢𝑗′′𝑦𝛼′′̃ )𝜕𝑥𝑗  Chemical species 2 – c 𝜕(?̅?ℎ̃)𝜕𝑡+𝜕(?̅?𝑢𝑗ℎ̃)𝜕𝑥𝑗= −𝜕𝐽𝑗ℎ𝜕𝑥𝑗+𝐷?̅?𝐷𝑡+ (𝜏𝑖𝑗𝜕𝑢𝑗𝜕𝑥𝑖)̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅+ 𝜌?̇?ℎ 𝑙𝑖𝑞 −𝜕(?̅?𝑢𝑗′′ℎ′′ ̃ )𝜕𝑥𝑗  Energy 2 – d  In Eq.2, the source terms 𝜌?̇?𝑙𝑖𝑞, 𝜌?̇?𝑗 𝑢,𝑙𝑖𝑞, 𝜌?̇?𝛼 𝑙𝑖𝑞, and 𝜌?̇?ℎ 𝑙𝑖𝑞 are due to mass, momentum, and energy exchange from droplets, possibly from diesel injection. 𝐽ℎ and 𝐽𝛼 are the enthalpy and species diffusive fluxes obeying the Fourier and Fick’s diffusion law. The term 𝜏𝑖𝑗 is the viscous stress tensor, 𝑊𝛼 and ?̅̇?𝛼 are the molar weight and the filtered molar rate of production (or consumption) of species 𝛼 due to chemical reactions. For a given reaction in Eq3, the rate of production of the species 𝛼 is determined by Eq.4, wherein the reaction rates 𝑘𝑙  are described by an Arrhenius correlation as in Eq.5.  Σ𝛼=1𝛼=𝑁𝑠𝜈𝑙𝛼′ 𝑀𝛼 ↔ Σ𝛼=1𝛼=𝑁𝑠𝜈𝑙𝛼′′𝑀𝛼  3 ?̇?𝛼 = Σ𝑙=1𝑙=𝐿 {(𝜈1𝛼′′ − 𝜈1𝛼′ ) [𝑘𝑙,𝑓(𝑇)Π𝛽=1𝑁𝑠 𝑐𝛽𝑣𝑙𝛽′− 𝑘𝑙,𝑟(𝑇)Π𝛽=1𝑁𝑠 𝑐𝛽𝑣𝑙𝛽′′]}  4 𝑘𝑙,𝑓(𝑇) = 𝐴𝑙,𝑓𝑇𝑏𝑙,𝑓 exp {−𝐸𝑙,𝑓𝑅𝑈𝑇}  5 The last terms on the right-hand-side of Eq.2 – a(b-d) are the sub-grid scale stresses which pose a closure problem on the system of equations and need modeling. A wide range of turbulence models use the concept of turbulent viscosity, Eq6-a, wherein the anisotropic turbulent stresses are assumed to be proportional to the filtered strain-rate tensor (?̃?𝑖𝑗), Eq.6-b.  16  𝑢𝑖′′𝑢𝑗′′̃ −23𝑘𝛿𝑖𝑗 = −2𝜈𝑡?̃?𝑖𝑗  6 – a ?̃?𝑖𝑗 =12(𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖)  6 – b  An extra transport equation for the turbulent kinetic energy is solved by using 1-equation (for TKE) LES model in OpenFOAM. This is recommended for engine simulations as the Sub-Grid Scale (SGS) kinetic energy affects the reaction rates [45]. In this LES approach, eddy viscosity and the rate of dissipation of TKE are scaled with the SGS kinetic energy itself, Eq.7 [54]. The dependency of the turbulent viscosity on Δ, the filter width, means that the effect of nonlinear advection terms on the LHS of Eq.2 grow by mesh refinement; thus more nonlinear chaotic motion is unfolded. In other words, LES approaches are grid-dependant; particularly, the solution approaches the DNS results as the filter and mesh size become smaller. Throughout this study, the spatial filter width (Δ) is (cell volume)1/3, and for the constant in Eq.7 CΔ =0.185 is used. 𝜈𝑡 = 𝐶Δ√𝑘𝑆𝐺𝑆Δ  7 – a 𝜖𝑆𝐺𝑆 =𝐶𝜖𝑘𝑆𝐺𝑆32Δ  7 – b  Given the highly non-linear dependence of the chemical reaction rates on the temperature and species compositions, Eq.4-5, calculating their filtered values poses a major modeling challenge on combustion simulations. A successful approach is the Conditional Moment Closure (CMC), in which the filtered rates of reactions are calculated based on the conditionally averaged temperature, density, and species given a conditioning variable. Since temperature and chemical compositions have strong dependence on the mixture fraction, it is used as the conditioning variable [55]. Mixture fraction, 𝜁, is a conserved scalar representing the mass fraction of the fuel stream (material exiting the injection nozzle), or equivalently the mass fraction of the carbon and hydrogen atoms from the fuel (and other atoms if the fuel consists of non-hydrocarbon species) in the cylinder. Notably, mixture fraction is an inherently conserved scalar without any sink or source terms resulting in its simplified transport equations. Invoking the CMC assumption, the conditionally filtered reaction rates are expressed in terms of the conditionally filtered 17  temperature, species mass fractions, and density, Eq.8. The instantaneous filtered rate of reaction at any location is then calculated given the conditionally filtered reaction rates, calculated in Eq.8, and the mixture fraction probability density function (pdf), Eq.9. ?̇?(𝑦𝛼, 𝑇, 𝜌)|𝜁̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ≈ ?̇?𝛼( 𝑦𝛼|𝜁̅̅ ̅̅ ̅̅ , 𝑇|𝜁̅̅ ̅̅ , 𝜌|𝜁̅̅ ̅̅ )  8 ?̅̇?𝛼(𝑥𝑘, 𝑡) = ∫ 𝑃𝑍(𝑥𝑘, 𝑡; 𝜁)(?̇?|𝜁)̅̅ ̅̅ ̅̅ ̅𝑑𝜁10  9  The probability distribution of the mixture fraction is obtained by solving transport equations for the filtered (mean) and variance of the mixture fraction [56,57]. Given its mean and variance and assuming an apriori function, the mixture fraction distribution is obtained. It is shown that assuming a β-distribution function could reproduce the DNS mixture fraction data for a given mean value and variance [56]. In the CMC approach, the conditionally averaged quantities are obtained by solving separate transport equations for each. However, this adds to the number of transport equations and the dimensionality of the system, compromising the computational time. To mitigate this, Bushe and Steiner suggested calculating the conditionally averaged quantities using an ensemble of points in a jet radial plane, using an apriori knowledge from DNS data sets that the conditionally averaged quantities show the least gradient in such planes. Knowing the averaged flow quantities from the transport equations 2a-d, the associated conditionally averaged quantities are calculated by inverting the following integral equation, Eq.10, wherein the subscript A and superscript (n) denote the ensemble of points and the n-th point in the ensemble, respectively [56,57]. Such conditionally averaged quantities are then used to calculate the conditionally averaged reaction rates in Eq.8. ?̃?(𝑥(𝑛), 𝑡) =  ∫ 𝑃(𝑥(𝑛), 𝑡; 𝜂) (𝑇|𝜂)𝐴,𝑡̅̅ ̅̅ ̅̅ ̅̅ ̅̅ 𝑑𝜂10  10 Despite the simplification offered by the CSE assumption, calculation of the reaction rates for a complex chemical mechanism is still time-consuming. To simplify this aspect of the combustion modeling, mechanism reduction methods are used. Details of the chemical mechanisms are explained in section 2.2; the mechanism simplification through TGLDM method and the link between CSE calculation and TGLDM tables are explained in section 2.2.3. 18  2.1.3 Wall treatment Fully resolving the near-wall flow imposes high costs on turbulent simulations; consequently turbulent wall models are adapted to aid maintaining a reasonable near-wall mesh resolution. Spalding’s law of the wall, the standard wall-treatment in OpenFOAM turbulence calculations, is used throughout this work with standard model constants [58]. Using this model, the first cell adjacent to the wall could be located within the buffer layer (y+ < 30), yet the effect of wall shear stress on the momentum could be captured [59]. For the turbulent kinetic energy at the solid boundaries, a zero gradient boundary condition is assumed. 2.1.4 Diesel injection  The diesel injection and subsequent spray breakup and atomization in GOLD are handled by “dieselFoam” solver. Several injector types are supported for diesel injection calculations, the type “commonRailInjector” is selected in the current work. By specifying the injection pressure (rail pressure) and nozzle diameter; the injection velocity is calculated according to Eq.11. The finite needle opening and closing durations are also accounted for through adjustable mass flow rate profiles. 𝑈 = √2(𝑃𝑖𝑛𝑗−𝑃𝑎𝑚𝑏)𝜌  11  OpenFOAM provides libraries linked with dieselFOAM to account for spray atomization and break-up. In GOLD simulations, primary atomization is not explicitly modeled; instead, the diesel is exiting the nozzle as “blobs” with 80% of the nozzle diameter. These blobs undergo Kelvin-Helmholtz instability causing the separation of primary drops. The radii of the formed droplets and the rate of change of blobs’ radii are determined from Eq.12-13 wherein ΛKH and τKH are the wavenumber and breakup time. The numerical values of coefficients B0 and B1 are 0.61 and 40 respectively [60]. 𝑟𝑐 = 𝐵0ΛKH   12 𝑑𝑟𝑑𝑡=𝑟−𝑟𝑐𝜏𝐾𝐻  13 19  𝜏𝐾𝐻 =3.726𝐵1𝑟Ω𝑅𝑇Λ𝐾𝐻  14  Secondary breakup, governed by the Rayleigh-Taylor (RT) instability, occurs if the instability wavelengths are smaller than the drop radii, Eq.15. Further breakup happens if the time after the start of RT instability is greater than the breakup time, Eq.16. The numerical value of model coefficients CRT and Cτ are 1 and 0.1 according to [60].  𝜆𝑅𝑇 =2𝜋𝐶𝑅𝑇𝐾𝑅𝑇  15 𝜏𝑅𝑇 =𝐶𝜏Ω𝑅𝑇  16  OpenFOAM’s standard Ranz-Marshal model for droplet heat transfer, and aerodynamic drag for droplet motion are used, Eqs.17,18. If the droplet Re is greater than a limiter value, Relimiter, a constant drag coefficient, Cd,limiter, is used; these constants have numerical values 1000 and 0.44 [61]. 𝑁𝑢 = 2.0 + 0.6𝑅𝑒0.5Pr0.33 17 𝐶𝑑 =24𝑅𝑒(1 + 1.67𝑅𝑒0.67)  18 2.1.5 Natural gas injection The gas injector uses zero-dimensional thermodynamics model to calculate the sound speed and Mach number at the nozzle outlet based on stagnation pressure at the injector gas plenum. The gas flow is subject to block heating in the plenum, thus its temperature at the injector nozzle is higher than the fuel rail. A default constant value of 370 K for the gas flow at the nozzle is assumed throughout this study. During the gas injection, the cylinder pressure varies from 100 bar up to 150 bar; as a result the injector-cylinder pressure ratio varies between 2.5 and 1.5. Given the critical pressure ratio for natural gas is 1.85, the nozzle is choked during most of injection. In reality this means a sonic (Mn=1) flow at the nozzle, along with expansion fans emanating from it corners to match jet pressure with the cylinder. The reflected expansion waves form compression waves, which for highly under-expanded nozzles leads to a normal shock in 20  the nozzle downstream, called the Mach Disk [62,63]. For moderate pressure ratios (below 3) Li et al. showed that the boundary reflected compression waves do not include shock and can repeat until dissipated by viscous effect [62].  In GOLD, The velocity at the nozzle is calculated based on the plenum-cylinder pressure ratio, Eq.19 (a-c). The natural gas specific heat ratio and gas constant are 1.3 and 506.5 J/kg.k, and once substituted in Eq.19-a, gives the nozzle Mach number followed by the calculation of nozzle velocity. The actual nozzle velocity is calculated by taking into account the nozzle blockage effect by discharge coefficient, Cd, Eq.19-c. Since only a slightly supersonic flow at the nozzle is observed, no CFD divergence is encountered; however, this does not hold for higher injection pressure ratios, e.g. very early at the compression stroke.  𝑃0𝑃= [1 +12(𝛾 − 1)𝑀𝑛2]𝛾𝛾−1  19 – a 𝑈𝑛∗ = 𝑀𝑛 × √𝛾𝑅𝑇𝑛  19 – b 𝑈𝑛𝑎 = 𝐶𝑑𝑈𝑛∗   19 – c      Figure 5- Pressure fluctuation in the nozzle near-field region 21  A rather coarse mesh is used in GOLD which cannot resolve the pressure waves; still, modest pressure waves are observed downstream of the nozzle, Figure 5. The gas injector file takes the total injection mass, start of injection, and needle ramp-up and ramp-down times; the needle closing occurs such to deliver the total mass by the end of injection. 2.2 Modeling combustion chemistry in GOLD 2.2.1 Natural gas combustion mechanism Natural gas is mostly a mixture of methane, ethane, and propane (CH4, C2H6, C3H8). The methane mole fraction is greater 90% [64]; hence, the combustion mechanism of methane is critical in the modeling of natural gas combustion [65–67]. It has been shown that the kinetic behaviour of the CH4 substantially varies in different temperature and pressure limits [68]. To study the behaviour of methane combustion under engine-relevant conditions, Huang et al. measured the ignition delay of lean, stoichiometric, and slightly rich mixtures in a shock tube facility. The shock-tube data were used to improve an earlier CH4 oxidation mechanism. Different oxidation paths at low (1050 K) and high (1250 K) temperatures were found as shown in Figure 6 [69]. The effect of heavier alkane compounds, namely ethane and propane, was investigated in a later work by Huang et al. With the same experimental method, the results were used for the modification of the GRI 2.1 mechanism to reproduce ignition delay data [70]. The C2-C4 reactions were taken from other studies, and the final mechanism contained 55 species and 278 elementary reactions. After fine-tuning the rate coefficients, the mechanism was able to predict the ignition delay of different mixture compositions. Ethane and propane were found to have accelerating effects on ignition through methylperoxy and methylhydroperoxy chemistry. This improved natural gas mechanism is used for combustion calculations in GOLD. 22    2.2.2 Diesel combustion mechanism To simplify the modelling, heptane is often used to represent the behaviour of the conventional diesel fuel. The heptane mechanism from Lawrence Livermore National Lab covers the pressures ranging from 1 to 40 atm, temperatures from 550 to 1700 K and a wide range of dilution ratios. This mechanism was improved to model the phenomena such as auto-ignition, cool flame and negative temperature coefficient (NTC) particular to diesel fuels [71]. LLNL diesel mechanism with 170 species and 1500 reactions is used for diesel combustion calculations in GOLD. 2.2.3 Simplification of reaction mechanisms Hydrocarbon combustion involves tens to hundreds of intermediate species and hundreds or thousands of reactions. Coupling such a complicated reaction mechanism with flow simulation in turbulent regimes calls for demanding computational resources and considerable human time for implementation. Therefore, many simplification in the reaction system is needed [72]. Figure 7 shows different reactions on a time spectrum, also compared with the time scales encountered in turbulent flows. Fast reactions with time scales less than 10-5 seconds soon reach partial equilibrium; also, intermediates with much higher rates of consumption than production obtain a steady-state value. On the other hand, for reactions with time scales comparable to eddy turn-Figure 6- Reaction pathways for methane oxidation at 40 bar [69] 23  over, the equilibrium or steady-state assumption causes major computational error; hence the physical and chemical phenomena should be resolved concurrently.    Given the drawbacks of steady-state and partial equilibrium assumptions, the idea of using tabulated species mass fractions prior to the flow simulation is used; replacing reaction rate calculations to table look-up. However, such tables have many dimensions in a reacting system, as large as the number of all reacting species. Thus, tabulating the species mass fractions requires excessive use of memory, an impractical strategy even for most powerful computers. Mass and Pope [73] showed the existence of converging paths in the state space of a reacting system. In other words, in spite of different initial conditions, the trajectories soon “bunch” around a common “intrinsic” manifold long before reaching equilibrium; thus greatly reducing the degree of freedom of the system. GOLD chemistry is simplified using a manifold approach. Figure 8 shows a manifold for the CO-H2-air system with reaction trajectories corresponding to very different initial conditions. The conventionally reduced mechanism only agrees with the manifold path in high, near-equilibrium temperatures. Manifold methods offer several benefits compared to conventional mechanism reduction, namely: Figure 7- physical and chemical time scales in a turbulent flame [21] 24  1- There are no assumptions such as partial equilibrium; therefore, the method is valid for a wide range of conditions. 2- Given the fact that the number of parameters is greatly reduced, the manifold could be tabulated prior to the simulation of the flow. The necessary calculations are done once prior to the simulations; therefore, significantly reducing the flow computations. 3- The method offers flexibility for the number of progress variables; depending on the desired accuracy, one or more variables could be used to parameterise the manifold.      In the generation of the “attracting manifolds”, it is guaranteed that progressions within the manifold are in the direction of slow time scales and equilibrate the reactions with the fast time scales. Two implementations are developed based on the manifold methods, Intrinsically Low Dimensional Manifolds (ILDM) [73], and Trajectory Generated Lower Dimensional Manifolds (TGLDM) [74]; the latter is implemented in the current solver, thus is explained in more detail here. For a 2-dimensional TGLDM, it is assumed that there are 2 independent variables that could describe the manifold. The rates of reactions for other species are calculated based on the value of these two variables evolving with time scales equal to the mixing process. Using TGLDM requires the following steps: Figure 8- Low dimensional manifolds for a CO-H2-air system. Reaction trajectories (arrows), conventionally reduced trajectory (thin line), and the manifold are shown in H2O-CO2 (left) and H-CO2 (right) planes [73]. 25  1- Choosing the appropriate variables to parameterize the manifold. In many free reacting jet studies, the mass fraction of the two species, CO2 and H2O are chosen as the progress variables [73,74]. In the current solver, the enthalpy difference and temperature are used as the progress variables. This choice is made to ensure smooth heat-release calculations. 2- Pre-processing or generating the manifold prior to any simulation. This is only done once and the data is stored in a table-like format. Tables are generated for different constant values of mixture fraction and the entire process is done for two pressure levels and 4 different initial oxygen concentration values, Table 3. The cylinder initial composition is calculated based on the engine oxygen-based overall equivalence ratio and the EGR flow rate. The calculation method and further details are provided in Appendix I, Figure 48, and sample compositions for three EGR levels at oxygen-based equivalence ratio of 0.62 summarized in Table 18. As shown in the appendix, the calculation method slightly over-predicts the oxygen mass fraction in the charge, but the error is less than 0.6% in all cases. Generating tabulated manifolds is a time-consuming task, and therefore, is carried out only for a finite number of oxygen mass fractions in the charge. In this work TGLDM manifolds are generated only for 7 different charge compositions, labeled I-VII in Table 3. For each simulation, the manifold with the closest level of oxygen to the calculated cylinder charge oxygen mass fraction is chosen for the combustion calculations. Table 3- Charge composition for the tabulated manifolds      I II III IV V VI VII O2 0.2330 0.2200 0.2149 0.2102 0.2050 0.2020 0.1983 N2 0.7670 0.7637 0.7625 0.7613 0.7600 0.7592 0.7583 CO2 0.0000 0.0090 0.0125 0.0158 0.0194 0.0215 0.0240 H2O 0.0000 0.0073 0.0101 0.0128 0.0156 0.0173 0.0194  As stated in section 2.1.2 the filtered flow quantities are calculated using the CSE assumption. When TGLDM method is used, the CSE equations should be solved only for the progress variables; knowing the conditionally averaged progress variables, the mass fraction of intermediate species are retrieved from the TGLDM tables. The mass fraction of intermediate species, along with other 26  conditionally averaged quantities are used for the calculation of the source and sink terms in the transport equations. Figure 9 describes the logical relation between different modules of the combustion computational model. Notably, the link between the TGLDM module and the CSE calculation for the chemical reaction rate calculations are shown in the figure. The CSE module calculates the conditional average of the progress variables, given the solution of transport equations from Open-FOAM and the mixture fraction PDF. By passing the progress variables to the TGLDM module, the intermediate’s mass fractions and other related quantities are looked up and passed to the main routine. The conditional reaction rates are calculated given the conditionally averaged quantities, which is then integrated over the mixture fraction space to yield the net reaction source terms in the transport equations.  Figure 9- schematic diagram of the different blocks in GOLD software, modified from Huang [52].  2.2.4 Soot chemistry and mathematical modeling Soot consists of carbon and small amounts of hydrogen. The evolution of soot is a complex process starting from gas phase chemistry which after several progressions results in the emergence of solid particles; different steps leading to the formation of soot particles are 27  depicted in Figure 10. The gaseous phase of soot formation starts with fuel pyrolysis; during which the ring-structured aromatic hydrocarbons are formed [23].    After initial formation, the PAHs form more rings and grow in size mainly through the Hydrogen Abstraction C2H2 Addition mechanism, HACA for short, as indicated in Eq.20 in a phenomenological form  [75,76]. 𝐶2𝐻2 + 𝑛𝐶𝑠 → (𝑛 + 2)𝐶𝑠 + 𝐻2  20  The aromatic compounds grow in size and mass by further chemical reactions and also through collision with each other. After reaching the size of 1-2 nm, such large compounds form soot nuclei. Soot nuclei grow in size via surface growth and coagulation, eventually forming spherical primary particles ranging from 20-70 nm in diameter. The primary particles also collide and agglomerate forming fractal-like soot particles with a diverse range in size from small hundred nanometers to a few microns. Meanwhile oxidation reactions also occur, burning the carbon atoms, therefore reducing the number density, size and mass of the particles. O2 and OH are reported as the main agents contributing to the oxidation of soot through reactions on the surface of soot and with its gaseous precursors [23,77]. Many approaches have been taken to simplify soot evolution to model these steps, such as linking the inception of soot nuclei to the incomplete pyrolysis of acetylene, Eq.21 [76]; further refinements are also applied to include benzene in the nucleation process [77]. Surface growth is modeled by HACA, Eq.20; while surface oxidation is modeled using Strickland-Constable’s [78] Figure 10- Phenomenological picture of soot formation [23] 28  O2 oxidation and Neoh’s OH oxidation [79]. These pathways for soot formation-oxidation are summarized in Eq.21-23.  𝐶2𝐻2 → 2𝐶𝑠 + 𝐻2 21 𝐶𝑠𝑜𝑜𝑡𝐻 + 𝑂2 → 𝐶𝑠𝑜𝑜𝑡−2𝐻 + 2𝐶𝑂  22 𝐶𝑠𝑜𝑜𝑡𝐻 + 𝑂𝐻 → 𝐶𝑠𝑜𝑜𝑡−1𝐻 + 𝐶𝑂 + 𝐻  23 To account for the agglomeration of the particles, balance equations for the number density of the soot primary and agglomerates are solved using the method of moments [80]. Even though different processes in the formation and oxidation of soot are modeled; these models are developed for specific flame configurations and are validated for the very same condition. Thus, applying them in environments other than those they were developed require lots  of tuning [81]. In empirical approaches, the formation and oxidation of soot are modeled by two global Arrhenius type reaction coefficients. Hiroyasu and Nishida [82,83] implemented this model for the simulation of diesel spray combustion in a direct injection configuration with the following formulation: 𝑑𝑚𝑠𝑓𝑑𝑡= 𝐴𝑓𝑚𝑓𝑔𝑃0.5 exp (−𝐸𝑆𝑓𝑅𝑇)  24 – a 𝑑𝑚𝑠𝑜𝑑𝑡= 𝐴𝑐𝑚𝑠𝑃𝑂2𝑃𝑃1.8 exp (−−𝐸𝑠𝑜𝑅𝑇)  24 – b 𝑑𝑚𝑠𝑑𝑡=𝑑𝑚𝑠𝑓𝑑𝑡−𝑑𝑚𝑠𝑜𝑑𝑡  24 – c Here 𝑚𝑠 is the soot mass, affected by a balance between formation (Eq.24 – a) and oxidation. Esf and Eso are the formation and oxidation activation energies and the pre-factors Af and Ac are typically used to tune the model. Due to its simplicity, the Hiroyasu model has gained popularity with its variants being used in many modeling studies, even those with detailed LES turbulence models [84]. The Hiroyasu model is also used for soot modeling in GOLD with the default values from [85]. While this results in reasonable predictions for a RANS turbulence model, the LES soot levels are two orders of magnitude smaller than the measured values, indicating that the pre-factor coefficients are heavily dependent on the turbulence model used. In this study a constant scaling factor Cs=100 is used in the post-processor to scale the final soot value. This constant is unchanged among all the simulation cases for proper comparison among different cases.  29  2.3 The Single-Cylinder Research Engine (SCRE)  The engine tests by Faghani and Mabson [48,49] that motivate this CFD study were performed on a 6-cylinder 4-stroke Cummins ISX 400 engine, modified to run with a single firing cylinder. Other 5 cylinders are motored, dummy injectors are installed, valves are bolted shut with the rocker arms removed, and the pistons are modified to run with single ring to lower the wall friction. The engine has 2.5 litre of displacement per cylinder, with a compression ratio of 17:1 and bore and stroke of 137 and 169 mm. The SCRE is connected to a water-cooled eddy-current dynamometer to absorb power by up to 150 kW. To overcome the friction load imposed by the reciprocation of other cylinders and other moving parts, an electric motor controlled with a vector drive is linked to the SCRE to provide additional power at low-load operations and dampen the torque pulses. The SCRE is designed to support independent control of intake and exhaust pressures and temperatures. A centrifugal compressor provides boosted air for the intake, modeling turbocharged engine conditions. The intake air enters an upstream surge tank where it mixes with a fraction of the recirculated exhaust gases (EGR) and its fluctuations are damped. Natural gas and diesel are provided at high pressures up to 30 MPa, with the diesel having slightly (1-2 MPa) higher pressure to avoided leaking gas in the diesel injector plenum and actuating the injector needles. The engine is fueled with Vancouver natural gas (typically 96% CH4, with 2% higher hydrocarbons and 2% N2 and CO2) and No.2 ULSD Diesel fuel. The engine exhaust enters a downstream surge tank at which point a fraction of the exhaust is passed to the intake stream passing an EGR regulator. More detailed specifications of the SCRE could be found at the references [47–49]. The SCRE uses a flush-mounted water-cooled Kistler 6067 piezoelectric pressure transducer for high-speed in-cylinder pressure measurement at every 0.5°. The data are then averaged over 45 cycles for mean pressure and rate of heat release calculations. The Kistler high-speed pressure transducer operates based pressure differential values; to obtain absolute pressure values, the sensor reading is pegged with a manifold piezoresistive pressure transducer at intake BDC of every cycle. 30  Before venting the engine exhaust gases, a small fraction of it is separated for gaseous and particulate emissions measurements. The gaseous emissions are analyzed by an AVL gaseous emission analyzer which measures wet (water vapour not removed) NOx and HC emissions and dry CO, CO2, and O2 emissions. The HC analyzer has the capability to measure both methane and non-methane HC emissions using Flame Ionization Detectors (FID). The NOx (NO + NO2) emissions are analyzed with a chemilumenescent detector, measuring the intensity of light emitted during the reaction of O3 with NO in a chamber. CO and CO2 are analyzed with the Non-Dispersive Infra-Red (NDIR) technique, using the infrared absorption band properties of CO and CO2. Finally O2 is measured using a paramagnetic oxygen sensor. Another portion of sampled emissions flows to the particulate measurement bench. Since the PM concentration is very high on the engine exhaust, the sample is diluted in an air-ejector mixing chamber to the ratio of 12:1. The dilution ratio is calculated by comparing a post-dilution CO2 reading with the AVL bench CO2 value. The PM mass is measured using multiple instruments with different techniques. A portion of the sample flows through a Tapered Element Oscillating Member (TEOM) for mass measurements. Another sample line passes through a denuder developed at UBC [47] to remove the PM volatile content for better light scattering measurement in the TSI Dusttrack DRX device. At certain operating points the PM population dynamics and detailed morphology analyses was also done with SMPS and TEM imaging technique.  The Coefficient of Variance (COV) of the SCRE outputs arises from errors in setting the inputs for each operating condition and instrumental drifts over time. Mabson [48] has investigated the COV of measurements over time at B-75 (or Mode-4, section1.1.2) operating conditions. The results from the work of Mabson are reported in Table 4. More detailed description of the SCRE facility, including the specification of fueling system, air handling utility, and emissions bench could be found in [47–49]. Table 4- B-75 Point over time coefficient of variance [47] RSTDs DRX TEOM NOx CO CH4 J36 19% 32% 5% 16% 9% 7-7-0 44% 38% 5% 24% 5%  31  The simulated injectors in this thesis include the production J-36 injector, the 7-7-0 injector, and the four paired-nozzle injectors. The arrangement of holes for the Reference 7-7-0 and the paired nozzle injectors are schematically shown in Figure 11 and the design specifications of all injectors are summarized in Table 5 [48]. The J-36 configuration is similar to the 7-7-0, except that there are 9 single gas holes and 7 diesel holes in the actual drawings; nevertheless, a 9-7 configuration cannot be modeled by a single-sector mesh. Therefore, a configuration with 9 gas and 9 aligned diesel holes is modeled instead. Given that the number of gas and diesel holes is different in the real configuration, the interlace angles between them vary between -20° to 20° and is different for each of the 9 gas jets. This might result in different time of the start of combustion in a real HPDI cylinder and a numerically simulated one; however knowing that, due to the sudden pressure and temperature rise, the ignition of one jet is rapidly followed by the ignition of the others in reality, a similar time of ignition for all jets can be assumed. This ignition time thought to be close to the ignition time for an aligned gas-diesel configuration; thus, the simulation error between a 9-7 and a 9-9 configuration is negligible.      Figure 11- Arrangement of holes in the baseline 7-7-0 (left) and paired-hole injectors (right). Figure adopted from Mabson [48]. 32  Table 5- Design specifications of the simulated nozzles Nozzle Production injector (J-36) Reference injector (7-7-0) S-H-S-A1 S-H-L-A2 L-H-S-A3 L-H-L-A4 Number of gas holes 9 7 14 14 14 14 Nominal gas hole D [mm] 0.65 0.73 0.52 0.52 0.57 0.57 Angle between gas jets [˚] - - 10 18 10 18 Gas flow [LPM] 61 61 61 61 73 73  2.4 Calculation of CFD inputs Data needed for initial and boundary conditions of the simulation are divided into three categories, summarized in Table 6. Table 6- Initial charge condition and injection specifications for GOLD 1. Direct from SCRE 2. Indirect from SCRE run files 3. Estimated from injector drawings and flow tests Speed 1500 T-90 440 [K] (CD) 0.8 P-90 4.5 [bar] YO2 0.205 Pilot injection rate shape, ramp (up, down) (0.22, 0.22) [ms] Mpilot 11 [mg/inj] YN2 0.7600 Gas injection rate shape, ramp (up, down) (0.6, 0.6) [ms] PPW 0.62 [ms] YCO2 0.0194 τpilot 0.35-0.65 [ms] Mgas 173.7 [mg/inj] YH2O 0.0156 τgas 0.6-0.9 [ms] Ppillot 27 [MPa]   Tpilot 320 [K] Pgas 25 [MPa]   Tgas (Sac, Nozzle) (320, 370) [K]  Inputs from the first category are used by GOLD without modifications. Inputs from the second category are calculated from the state of the charge at -90 aTDC and the EGR-EQR levels. The                                                           1 Small-Hole, Small-Angle 2 Small-Hole, Large-Angle 3 Large-Hole, Small-Angle 4 Large-Hole, Large-Angle 33  details of these calculations are provided in Appendix I. There is no direct method to measure the third category in the SCRE facility; so, most of these inputs are adjusted according to the data provided by Westport Innovations for the production injectors. However, the pilot and gas needle lift delays are further adjusted to match the phasing of combustion. 2.5 CFD post-processing methods The CFD results contain detailed information about the flow and emissions in different regions of the combustion chamber; further data reduction is needed to get meaningful results, and this is described below.  2.5.1 Heat release rates The rate of heat release in the cylinder is obtained from the cylinder pressure following the same method as in experimental analyses. The first law of thermodynamics is applied to the cylinder as a control volume and the rate of heat release is obtained to balance to boundary work due to pressure increase and cylinder volume change [86]. In Eq.25, γ is the poly-tropic exponent of the mixture with the constant value of 1.3 in experimental data processing and numerical post-processing routines. This poly-tropic exponent is slightly smaller than the average specific heat ratio of air for an isentropic compression mainly due to heat transfer to the cylinder wall and piston. 𝑑𝑄𝑛𝑑𝜃=𝛾𝛾−1𝑝𝑑𝑉𝑑𝜃+1𝛾−1𝑉𝑑𝑝𝑑𝜃  25  To reduce the pressure differentiation noise in the experimental post-processing, the HRR data for 45 consecutive cycles are averaged to reduce the pressure differentiation noise. On the other hand, only one cycle is simulated in the CFD and multi-cycle averaging is not possible. Thus in the CFD post-processor, the HHR data derived from the pressure differentiation are filtered using a moving 5-point averaging technique to reduce its noise. This procedure slightly left-shifts a HRR ramp-up curve, but as is shown in section 3.1 the error is less than 1° of crank angle and therefore can be tolerated when comparing CFD results with experiments. 34  2.5.2 Equivalence ratio – temperature maps The use of Equivalence ratio-Temperature maps is reported in many studies in diesel literature; especially, when new injection strategies are introduced [43,44]. Equivalence ration, EQR, or φ is the ratio of actual fuel in the mixture to the stoichiometric fuel needed; in terms of air-fuel ratio, this is expressed in Eq.26. The analyses using φ-T maps consist of a background image where relevant emission contours are marked in different regions of the map. The foreground of these figures shows the status of the cylinder fuel on the map with which reasonable qualitative comparisons can be made among different injection strategies. Φ =𝐴𝐹𝑠𝑡𝐴𝐹𝑎𝑐𝑡  26 In this study, the φ-T maps are generated using the GRI 3.0 mechanism implemented in CANTERA, the reactor network for the emissions contours was generated by Ehsan Faghani, and simulated a 3-ms combustion of natural gas in a perfectly stirred reactor at constant temperature to generate the background emission contours for C2H2, CO, and NO. The C2H2 is chosen as an indication of high PM formation zones. In the foreground, the mixture is binned into equivalence ratio and temperature bins and normalized to generate a 2-D distribution. A typical φ-T map used for the mixture analysis is shown in Figure 12. One novel feature used in this study is generating the differential maps such that the difference of the mass between two cases at each φ and T is color-coded in the map and studied. This is a convenient method to highlight the differences in cases where the maps look very similar otherwise. The post-processing code to generate the φ-T maps is provided in the appendix II. 35   Figure 12- Equivalence ratio-temperature maps generated using CANTERA. Regions of high C2H2, CO, and NOx are marked. 2.5.3 Equivalence ratio and mixture fraction distribution Complementing the previous analysis, the distribution of the fuel mass over a wide range of the equivalence ratios (0-10) and the entire range of the mixture fraction (0-1) is calculated at each crank angle. This distribution is particularly helpful when analyzing the paired-hole jets in which the amount of mixture and its distribution in too-rich and slightly-rich zones are of focal interest. The description is mathematically calculated as shown in Eq.27. Similar method for analyzing mixture properties for diesel group-holed injectors is also reported in literature [25], in which the diesel injection into a high-pressure chamber (no combustion) was analyzed using the Light Scattering-Absorption (LSA) method to obtain the mixture equivalence ratios. However, that study was done for a configuration outside the engine chamber, therefore neglecting the effect of pressure rise, charge flow and combustion. 𝑓(𝜁, 𝜃) =∑ (𝜌𝑖𝑉𝑖𝜁𝑖|𝜁≤𝜁𝑖<𝜁+𝑑𝜁)𝐶𝐹𝐷 𝑐𝑒𝑙𝑙𝑠𝐼𝑛𝑗𝑒𝑐𝑡𝑒 𝑔𝑎𝑠 𝑎𝑡 𝜃∘ 𝑎𝑇𝐷𝐶  27  Figure 13 shows a sample pdf plot for the baseline and a negative PSEP injection at 10° after Start of Injection (aSOI). In this sample, the negative PSEP injection produces a leaner mixture for the same time aSOI shown as a peak in the stoichiometric-slightly rich zone followed by a sharp 36  reduction and a long narrow tail in the rich side. Mixture pdf’s are particularly used for the study of paired jets in section 3.2.2.    2.5.4 Rich and lean mixture evolution, air entrainment The evolution of the rich mixture with equivalence ratio between 2 and 6 is important to the formation of soot. Particularly in the non-premixed combustion, this fuel has a sufficiently high temperature to promote the soot nucleation and surface growth reactions. The normalized mass of rich mixture in the chamber, Zrich in Eq.28, is compared for different injectors and combustion strategies. In this equation, 𝜁 and φ are the local mixture fraction and equivalence ratio respectively, and 𝑚𝑓 is the mass of total injected fuel. 𝑍𝑟𝑖𝑐ℎ =∑ (𝜌𝑖𝑉𝑖𝜁𝑖|2<𝜙<6)𝐶𝐹𝐷 𝑐𝑒𝑙𝑙𝑠𝑚𝑓  28  Furthermore, it is useful to know how the lean mixture evolves during the combustion event to assess the performance of the injectors in producing a lean mixture. This is of particular interest when comparing the Zlean values, Eq.29, at the ignition time to compare the level of premixing in partially premixed strategies. Zrich and Zlean take values between 0 and 1, depending on the amount of fuel existing in the region associated with each. Zlean is a monotonically increasing Figure 13- PDF of equivalence ratio for the baseline and negative PSEP injection. Both plots show the state of mixture at 10° after Start of Injection (aSOI). 0 2 4 6 8 1000.20.40.60.8EQRPDF  J-36, PSEP = 0.3 msJ-36, PSEP = -0.9 ms37  function approaching 1 by the end of the cycle as the injected fuel is completely mixed with the charge. On the other hand, Zrich rises as the injection starts and then drops to zero shortly after the end of injection. 𝑍𝑙𝑒𝑎𝑛 =∑ (𝜌𝑖𝑉𝑖𝜁𝑖|𝜙 < 1)𝐶𝐹𝐷 𝑐𝑒𝑙𝑙𝑠𝑚𝑓  29  To quantify the amount of air entrained into the jet, a normalized parameter called Zair is defined in a manner similar to Zrich and Zlean. That is, Zair is the mass of air inside the stoichiometric jet envelope normalized by the total injected fuel, or in mathematical terms, Eq.30. Zair is used in section 3.2.2 to study the early air entrainment into the base of paired jets. 𝑍𝑎𝑖𝑟 =∑ ((𝜌𝑖𝑉𝑖(1−𝜁𝑖))|𝜙>1)𝐶𝐹𝐷 𝑐𝑒𝑙𝑙𝑠𝑚𝑓  30  2.6 Numerical verification and sensitivity study Validation of LES CFD is an active area of research. While in laminar or RANS turbulent simulations the grid size mainly affects the numerical error, for LES simulations the error associated with the Sub-Grid Scale (SGS) motion is also dependent on the grid size, especially when the filtering operation is done in an implicit manner. In other words, resolving the mesh in a LES simulation results in more resolved eddies and thus more physical features, while the numerical error associated with the discretization is reduced as well. Therefore, the classical Richardson extrapolation approach cannot be simply used to obtain the numerical accuracy [87]. A key indicator of the quality of a LES simulation is the ratio of the resolved (grid level) to the total (resolved + SGS) kinetic energy, Eq.31. Pope has suggested 80% or more of the kinetic energy should be resolved for an LES simulation, while 95% or higher ratios could be deemed as DNS [53]. 𝐼𝑄𝑘 =𝑘𝑟𝑒𝑠𝑘𝑡𝑜𝑡=𝑘𝑟𝑒𝑠𝑘𝑟𝑒𝑠+𝑘𝑆𝐺𝑆  31  38  KSGS in the denominator is obtained by solving an extra equation for the TKE, but obtaining the resolved turbulent kinetic energy could be challenging. Ensemble or phase averaging has been used in the literature to obtain the resolved kinetic energy, with the former being extensively used in the stationary configurations; whereas the latter is recommended for non-stationary, in particular reciprocating engines. In phase-averaging method, velocity vectors are averaged at a certain crank angle for several consecutive cycles, then the variation about such phase-averaged velocity determines the resolved (grid-level) turbulent kinetic energy. The level of resolved turbulent kinetic energy also indicates the level of turbulence-induced cyclic variability. Unfortunately, both methods are impractical in the current simulations, wherein only the closed portion of a single cycle is simulated; thus neither ensemble-averaging nor phase-averaging are practical. In the current work, mesh sensitivity study is conducted to investigate the effect of grid resolution on the emissions, penetration and mixture evolution.   39  3 Results First, the accuracy of GOLD is tested by comparing pressure and heat release predictions to the experimental measurements. Successively refined meshes are used to examine the sensitivity of emissions, penetration, and mixture evolution to the grid size. The trends of emissions under incremental change of inputs are studied to investigate the variability in CFD-predicted emissions about the average trends over a finite interval. Finally, a compilation of CFD predictions are showed on PM-CO and PM-NOx maps and are compared with their experimental counterparts. CFD is then used to study the performance of Paired-Hole Injection, Late Post Injection (LPI), and Slightly Premixed Combustion (SPC). The mixing, penetration, and combustion behaviour are studied with the tools introduced in section 2.5, and the results are compared with baseline cases. 3.1 GOLD validation In this section, the pressure trace and heat release rate (HRR) curves are compared against their experimental counterparts to match the phasing of combustion with the baseline experiments. The Baseline operating condition is the ESC mode-4 (B-75) which was described in section1.1.2 and is fully quantified for HPDI operation on SCRE in Table 7. The performance at B-75 of the two single-hole nozzles, 7-7-0 and J-36 production injectors is simulated using successively refined CFD grids; paired-hole injectors are studied in detail in section 3.2, but for mesh-sensitivity study, some results are compared with the baseline points in this section. Table 7- Operating parameters at ESC mode B-75 (Mode-4) ESC Mode Speed (RPM) Load – GIMEP [bar] EQR (O2-based) EGR GRP [MPa] Diesel PPW [ms] 50% IHR PSEP B-75 1500 16.6 0.62 18% 25.4 0.62 11 0.3  3.1.1 Pressure, heat release rate, and emissions The inputs for CFD calculations are calculated based on Table 7 summarized in Table 8.  40  Table 8- CFD inputs for baseline simulations  Pilot Gas Charge state at -90° aTDC SOI PPW τ SOI GPW τ P [bar] T [K] YO2 YN2 YCO2 YH2O 7-7-0 -17 0.84 0.61 -8 2.39 0.74 4.48 440 0.205 0.76 0.0194 0.0156 J36 -17 0.84 0.67 -7 2.43 0.88 4.38 439 0.205 0.76 0.0194 0.0156  Figure 14 shows the pressure and the heat release rate (HRR) predictions for the 7-7-0 and J-36 injectors at mode B-75. As is shown in this figure, the apparent start of ignition for pilot and gas, as attained by adjusting injection delays, match the experimental SCRE curves. The combustion timing, defined as the time of 50% IHR and marked with black dots on both curves, is also the same for both cases indicating ignition and the overall combustion phasing are matched simultaneously.     Figure 14- Pressure and rate of heat release predictions at the baseline operating conditions for 7-7-0 injector (left) and J-36 injector (right).  -60 -40 -20 0 20 40 60 80050100150200CA []P [bar]  7-7-0, SCRE7-7-0, GOLD-60 -40 -20 0 20 40 60 80050100150200CA []P [bar]  J-36, SCREJ-36, GOLD-60 -40 -20 0 20 40 60 80-5005010015020010.8810.5CA []HRR [kJ/(m3.deg)]  7-7-0, SCRE7-7-0, GOLD-60 -40 -20 0 20 40-5005015020010.8210.5CA []HRR [kJ/(m3.deg)]  J-36, J-36, 41  An error of 10-15 bar (equivalently 7-10%) is observed between the peak of experimental and CFD pressure trace curves. First, the Low Heating Value (LHV) of n-heptane is slightly (about 5%) more than the Ultra Low-Sulfur Diesel (ULSD) fuel; thus slightly more heat is released for the same mass of pilot in the CFD. More released energy before the Top Dead Centre (TDC) causes higher rate of pressure increase observed for pilot combustion. Second, 5-point averaging method effectively advances the heat release curve by 0.5° to 1° which, although within the range of uncertainty of the experimental combustion timing, contributes to slightly more combustion occurring prior to TDC, thus causing higher peak pressure. Since the common criteria of matching the injected fuel mass, and combustion timing based on the 5-point averaged HRR curves are used throughout this thesis, the drift from the experimental pressure is similar among all cases.  The effect of injection timing on the heat release prediction is shown in Figure 15, wherein the previous cases are compared with cases of smaller injection delay. As observed, the start of injection (GSOIphys) has a direct impact on the phasing of HRR curve: earlier injection results in the same level of advancement in the HRR ramp-up.     SCRE operation using 7-7-0 and J-36 injectors is simulated for 3 EGR levels to investigate the ability of GOLD chemistry to capture the experimentally-observed emission trends. Experimental conditions are summarized in Table 9 and CFD inputs are summarized in Table 10. The first, second and third rows of each injector in Table 10 correspond to 0%, 18% and 25% EGR flow rates; and as such, the O2 mass fraction reduces when migrating from low to high EGR cases. To keep the same EQR level among all simulations, the mass of ingested O2 should be similar in all Figure 15- Effect of injection timing of HRR curves for 7-7-0 (left) and J-36 (right) injectors. -60 -40 -20 0 20 40 60 80-5005010015020010.910.5CA []HRR [kJ/(m3.deg)]  9.57-7-0, SCRE7-7-0, Timing1, GOLD7-7-0, Timing2, GOLD-60 -40 -20 0 20 40 60 80-5005010015020010.88CA []HRR [kJ/(m3.deg)]  10.5SCREJ-36, T ming 1, GOLDJ-36, T ming 2, GOLD42  cases resulting in higher intake charge mass for higher EGR levels. This causes an increased pressure and slightly reduced temperature for the charge at the beginning of the simulation. Table 9- SCRE commanded inputs for EGR sweep test Injector ESC mode PSOI PPW GSOI GPW EQR [O2 based] EGR [%] 7-7-0 B-75 -22.5 0.62 -14.6 1.62 0.62 0, 18, 25 J-36 B-75 -23 0.62 -15 1.8 0.61 0, 18, 25  Table 10- CFD inputs for EGR sweep simulations  Pilot Gas Charge state at -90° aTDC PSOI PPW τ GSOI GPW τ P [bar] T [K] YO2 YN2 YCO2 YH2O 7-7-0 -17 0.84 0.61 -8 2.39 0.74 3.82, 4.48, 4.5 431 441 428 0.233 0.205 0.198 0.7670 0.7600 0.7583 0.0 0.0194 0.0240 0.0 0.0156 0.0194 J-36 -17 0.84 0.67 -7 2.43 0.88 3.97, 4.38, 4.63 444 439 437 0.233 0.205 0.198 0.7670 0.7600 0.7583 0.0 0.0194 0.0240 0.0 0.0156 0.0194  Consistent with the experimental trends, increasing the EGR flow rate reduces O2 concentration, resulting in an overall richer combustion and hence higher rate of soot emission, as shown in Figure 16-top row. In spite of correct predictions of the soot emission trends, the percentage increase of soot from low to high EGR is under-predicted. As mentioned in section 2.2.4, the 2-equation Hiroyasu model is used in this study with model’s coefficients kept equal to the default values but a single calibration constant (C=100) is applied in the post processor. The results could be scaled for a baseline point to match the experimental data, but by maintaining constant model parameters throughout the whole work a fair ground is provided for cross-comparing the emission levels; indeed under-predictions are also observed in some cases.    43          High EGR is typically used to control the NOx emission by reducing the flame temperature as explained in section 1.3.2. The simulated and measured NOx values are remarkably close. Predicted CO does not follow the measurements, keeping an almost constant value throughout the simulations. CO is a major product in the final stages of the hydrocarbon oxidation, requiring the inclusion of the various oxidation pathways from the pure fuel to the final oxidized products. The natural gas combustion mechanism implemented in GOLD is validated based on experimental shock tube data [52]. However, as in many other studies, the ignition delay data and some intermediate species are used as the validation parameters. While correct ignition delay guarantees correct radical formation to initiate the combustion, it does not guarantee that 0 18 250.020.040.060.080.1EGR [%]PM [gr/kW.hr]  J-36, GOLDJ-36, SCRE0 18 250123456EGR [%]NOx [gr/kW.hr]  7-7-0, GOLD7-7-0, SCRE0 18 250123456EGR [%]NOx [gr/kW.hr]  J-36, GOLDJ-36, SCRE0 18 250246810EGR [%]CO [gr/kW.hr]  7-7-0, GOLD7-7-0, SCRE0 18 25024681EGR [%]CO [gr/kW.hr]  J-36, GOLDJ-36, SCREFigure 16- PM (top), NOx (middle), and CO (bottom) trends during an EGR sweep test for 7-7-0 (left) and J-36 (right) injector. Black dots show SCRE measurements and green dot the GOLD simulations. The experimental points are provided by Faghani and Mabson. [48, 49]. 0 18 250.020.040.060.080.1EGR [%]PM [gr/kW.hr]  7-7-0, GOLD7-7-0, SCRE44  all pathways leading to CO formation are considered throughout the combustion mechanism. Investigating the details of the combustion mechanism is beyond the scope of this work; however, the current results indicate that further improvements to the mechanism is needed for better prediction of CO. NOx and PM on the other hand are separated from the hydrocarbon mechanism; the former is modeled with an extended Zeldovich mechanism and the latter with an empirical two-equation Hiroyasu model. A compilation of the previous CFD results are summarized in the compact forms of PM-NOx and PM-CO emissions maps. Typical PM-NOx trade-off behaviour for EGR sweeps is observed for both J-36 and 7-7-0 injectors. Previous experimental research showed a strong correlation between the PM and CO emissions exist across a wide range of the operating conditions; more recent research showed this correlation is even observed for different injectors [48,49]. Therefore, a compilation of CO and PM emissions are plotted on a log-log scale in Figure 17 wherein the data are shown to be closely scattered about a line with a constant slope. Also shown in Figure 17, are the SCRE measurements corresponding to the simulations. Notably, a steeper PM-CO behaviour is observed; more importantly, PM and CO show a stronger correlation with the PM exponent of 0.57 from SCRE compared to 0.13 from GOLD.   45        3.1.2 GOLD sensitivity to mesh The phenomena occurring in a direct injection internal combustion engine involve time and length scales of different orders of magnitude. Namely, the charge motion inside the cylinder is regarded as a low Reynolds flow compared to the high-speed gas jet during injection. Therefore, to study the effect of mesh on the CFD-predicted performance and emissions, a finer mesh is used close to the jets, Figure 4. As explained in section 2.1.1, to keep a feasible computational time, the mesh is only further refined during injection, i.e. in -10° to 10° aTDC interval. As explained in section 2.6, LES numerical accuracy cannot be assessed with classical Richardson extrapolation method. Especially when emissions are concerned, highly non-linear reaction rates and excessive dependence on the flow mixing is involved. Nevertheless, since emissions are compared between different injection strategies, their variation with mesh resolution is of interest and thus is investigated in Figure 18. The circumferential mesh resolution on the injector 0 0.02 0.04 0.06 0.08 0.102468PM [gr/kW.hr]NOx [gr/kW.hr]  J-367-7-0GOLD0 0.02 0.04 0.06 0.08 0.102468PM [gr/kW.hr]NOx [gr/kW.hr]  J-367-7-0SCRE10-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J-367-7-0CO=PM0.13GOLD10-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J-367-7-0CO ~ 0.57SCREFigure 17- Compilation of the CFD data for J36 and 7-7-0 injectors on PM-NOx (top) and PM-CO (bottom) maps. Data from GOLD simulations (left) are compared with SCRE (right) experiments.  46  nozzle, dθ, is chosen as the mesh size indicator. The effect of successive mesh refinement on PM, NOx and CO emissions for 4 injectors, J-36, 7-7-0, L-H-S-A, S-H-S-A1, are investigated. The S-H-S-A and L-H-S-A emit higher PM on average and the variation between coarse and fine meshes is less than 10%, except for very fine L-H-S-A mesh where a sudden PM jump is observed. Not much variation in the PM predictions for J-36 injector is observed either, with PM emissions keeping an almost constant value of 0.04 gr/kW.hr for the three mesh resolutions. On the contrary, the 7-7-0 PM emission demonstrates much higher variations, by up to 30% compared to the baseline (the coarsest mesh).  It is critical to note here that only a single realization of the simulation is used in each case.  As in reality, one expects cyclic variability in flows and emissions. As the mesh resolution is changed, one expects a random change in results (related to the cyclic variability) independent of the mesh effect.    Other emissions, NOx and CO, show less variations with the mesh resolution, demonstrating a maximum of 20% change for the L-H-S-A injector. There is much more confidence in the NOx formation mechanisms than PM; particularly due to the simplicity of the thermal mechanism, its implementation poses less challenge in the CFD codes. As will be shown in the next sections, its readings also demonstrate close predictions of the measured readings. Less varying levels of NOx indicate some type of convergence for the quantities governing NOx formation, namely stoichiometric mixture and temperature in the reacting zone. CO also shows convergence; nevertheless, as demonstrated in PM-CO emission maps, less CFD-CO variation among different operating conditions is generally observed. The error bars are obtained using the Standard Error of Mean (SEM) calculation, Eq.32, based on the different emission predictions for each injection as shown in Figure 18. In Eq.32 𝑛 is the number of data available for each injector, and 𝑆𝑇𝐷 is their standard deviation. The S-H-S-A injector has the highest mass of soot, followed by the L-H-S-A and 7-7-0 injectors. As will be explained in more details in section 3.2, this is consistent with the SCRE tests as well; however, the ratio of S-H-S-A soot emission to the other injectors is under-predicted by CFD. The J-36                                                           1 Detailed information about the abbreviated notations and injectors’ specifications could be found in section2.3, Figure 11, and Table 5. 47  injector has the least variation of soot (its SEM is very small, the error bar is almost like a dot), but the other injector have very similar absolute errors.      𝑆𝐸𝑀 =𝑆𝑇𝐷√𝑛  32 0.5 1 1.500.050.10.15Circumferential mesh resolution [deg]PM [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AStandard error0.5 1 1.500.511.5Circumferential mesh resolution [deg]NOx [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AStandard error0.5 1 1.5012345Circumferential mesh resolution [deg]CO [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AStandard errorFigure 18- Sensitivity of CFD emissions to the mesh circumferential resolution for J36, 7-7-0, and two paired-hole injectors. Top: PM, Middle: NOx, Bottom: CO. 48  Due to the high level of variation in the CFD-PM predictions, alternative methods are developed to study the mixture formation in the cylinder, described in sections 2.5.3-2.5.2. As such, the penetration of the fuel jet is compared among different mesh resolutions in Figure 20. The penetration distances are measured on the gas injection plane, 10° below the fire deck, shown in Figure 19-left. The sudden jump observed in the penetration rate around 0.9 ms-aSOI is coincident with the impingement of the jet on the piston lip. The impingement happens on the edge of the piston bowl (piston lip), near the squish region; consequently, the jets split into four main directions: swirl-wise, counter swirl-wise, towards the piston bowl, and towards the squish region. In this work, the swirl-wise jet penetration is taken to be representative of the jet penetration after impingement as illustrated in Figure 19-right. For the most part, close similarity is observed between the penetration predictions, especially in the pre-impingement period. More uncertainty is observed in the penetration distance after the impingement; this is mostly due to the variation in the exact impingement spot and penetration progress which might slightly lay off-plane in penetration calculations. Nevertheless, the error is within 6-8% in the worst case for J-36 prior to impingement. The somewhat off penetration trend for very refined mesh in J-36 simulation could be attributed to the chaotic shape of the jet due to the more resolved, smaller eddies, Figure 21.     The standard error for the penetration of the single jets is less than 6% and takes its highest values just after the impingement, as indicated by Table 11. Also, according to the table, the J-36 jet penetrates slower than the 7-7-0 jet. Based on the theory of free jets this is justified given that J-36 has a smaller nozzle compared to the 7-7-0 injector. The multi-mesh error bars for Figure 19- Gas injection in HPDI engine. The gas jet penetration on the injection plane (left), and centreline vertical plane (middle) are marked. The penetration calculation after impingement is explained schematically (right). 49  penetration are used for a better comparison between the jet penetrations from different nozzles.          0 0.5 1 1.5 2020406080Time after start of injection [ms]Jet penetration [mm]  J-36, Theta resolution: 0.9 degJ-36, Theta resolution: 0.68 degJ-36, Theta resolution: 0.54 deg-8 -4 0 4 8CA aTDC [deg]0 0.5 1 1.5 2020406080Time after start of injection [ms]Jet penetration [mm]  7-7-0 l ti : 1.23 deg7-7-0, Theta resolution: 0.78 deg7-7-0, Theta resolution: 0.62 deg-8 -4 0 4 8CA aTDC [deg]Figure 20- Penetration sensitivity to the mesh resolution during -10° to 10° interval for J36 (top), 7-7-0 (middle) and paired-hole(bottom) injectors. 50   Table 11- Mean and standard error for the penetration of the jets based on multi-mesh simulations Time aSOI [ms] J-36 7-7-0 Mean Standard Error Mean Standard Error 0.3 10.84 0.12 11.21 0.21 0.6 26.38 1.05 31.16 0.39 0.9 35.41 1.12 37.27 0.26 1.2 47.35 2.65 56.1 2.28 1.5 56.26 2.21 59.81 3.31   Figure 21- Comparison between 7-7-0 jets with coarse (dθ = 0.78°, left) and fine (dθ = 0.62°, right) meshes. Some smaller eddies are marked in the right-hand side figure.  The effect of mesh refinement on mixture production of J36 and 7-7-0 injectors is shown in Figure 22. Very small difference is observed between the rich and lean mixtures for three successively refined J-36 meshes; the largest difference is observed at the beginning of the injection closing for Zrich parameter. Zlean varies slightly among the three meshes only late in the expansion stroke. For 7-7-0 injector, larger difference is observed especially in the evolution of Zlean. Intermediate mesh size, blue dashed line, shows the most desirable mixing behaviour, i.e. fastest rate of growth of the lean mixture. The difference supports the penetration and emissions observations; 51  that is, the intermediate mesh shows the fastest rate of penetration, lean mixture production, and soot emission. While Zrich is zero, the difference between the Zlean plots indicates that three different cases have different mass of fuel in the slightly rich region, 1<φ<2. Existence of such mixture even very late in the cycle demonstrates the sensitivity of late mixing to the mesh size. Overall, it can be said that Zrich is a more robust feature for mixture analysis.    Figure 22- Effect of mesh resolution on mixture production for J36 (top) and 7-7-0 (bottom) injectors. 3.1.3 GOLD sensitivity to input parameters and discussion of cyclic variability As noted above, only single simulations were done for each mesh, but one expects inherent variability from LES as one does from a real turbulent flow.  A simple way of running different realizations of the same conditions is test the sensitivity of the output to very small changes in input parameters, such as injection pressure and timing. For a baseline injection pressure of 25 -10 0 10 20 3000.20.40.60.81Crank angleZrich, Zlean  Resolution: 0.9 degResolution: 0.7 degResolution: 0.55 degZleanZrich-10 0 10 20 3000.20.40.60.81Crank angleZrich, Zlean  Resolution: 1 degResolution: 0.8 degResolution: 0.6 degZrichZlean52  MPa, the injection pressure is incrementally changed in both directions at steps of 4-bar, and the results are shown in Figure 23. The overall effect of pressure increase is a net reduction in the PM; however, large variations are observed around the trend-line, even with opposite trends in some short intervals. A decrease in PM for higher injection pressures is expected due to higher jet-induced kinetic energy and therefore enhanced mixing [13], and faster penetration with subsequently improved cylinder air utilization and shorter residence time in soot-forming zones. Modest reductions in NOx are observed which can be attributed to overall leaner combustion due to enhanced mixing as well. The CO shows no significant change throughout this study. Figure 24 shows change in the emissions when the injection timing is infinitesimally varied in steps of 0.2° (equivalent to 20 µs) in a 2° interval with the baseline GSOI being -7° aTDC. Like injection pressure, the overall effect of timing is observed with random variation about the trend-line. Retarding the combustion event has an overall soot reduction effect by shifting the high-temperature combustion around TDC further into the expansion stroke. Due to the same temperature reduction effects, the NOx emissions are also reduced in later timings. No meaningful change in the CO level is observed as before. It is worth noting in experimental measurements, the combustion timing is changed by retarding the injection, yet keeping a constant PSEP value; while in CFD study the GSOI is changed independently for a fixed Pilot timing.   53        Figure 23- Effect of Injection pressure on PM (top), NOx (middle) and CO (bottom) emissions. Injection pressure is increased in 4-bar increments. -20 -10 0 10 2000.020.040.060.080.1Injection pressure [+bar]PM [gr/kW.hr]-20 -10 0 10 2000.511.52Injection pressure [+bar]NOx [gr/kW.hr]-20 -10 0 10 200123456Injection pressure [+bar]CO [gr/kW.hr]54       Figure 24- Effect of combustion phasing (left) and injection timing (right) on PM (top), NOx (middle) and CO (bottom) emissions. The gas injection timing is increased by 0.2° CA in each step, but the 50% IHR timing is resolved in a 0.5° intervals. -1.5 -1 -0.5 0 0.5 1 1.500.020.040.060.080.1GSOI - GSOIbasePM [gr/kW.hr]-1.5 -1 -0.5 0 0.5 1 1.50.70.80.911.11.21.31.41.5GSOI - GSOIbaseNOx [gr/kW.hr]-1.5 -1 -0.5 0 0.5 1 1.5012345GSOI - GSOIbaseCO [gr/kW.hr]55  3.2 Analysis of the paired-nozzle injectors In this section, the performance of a series of paired-nozzle injectors for HPDI combustion is numerically investigated using GOLD. Firstly, the predictions of heat release rate for the four paired-nozzle injectors are presented and compared with a baseline single-hole injector under ESC-B75 operating condition. This is followed by CFD predictions reduced to 5 different metrics: 1- Penetration of the jets 2- Normalized entrainment rate (Zair) 3- Evolution of rich and lean mixtures (Zrich and Zlean) 4- Mixture distribution in φ-T maps 5- Distribution of fuel in mixture fraction and EQR space These metrics target the mixture preparation inside the cylinder from different perspectives to provide insight into the air-fuel mixing in paired-nozzle injection. Lastly, the effect of the attachment of the jets to the cylinder head is investigated by increasing the injector protrusion into the cylinder. 3.2.1 Heat release rates and final emissions at mode B-75 To compare the performance of different nozzle types, the phasing of combustion is adjusted to match the SCRE heat release curve; this is necessary to confirm the similarity of the charge state at the time of injection. The SCRE commanded parameters and the inputs for GOLD are summarized in Table 12 andTable 13. Table 12- SCRE commanded parameters for operation under mode B-75 for all paired nozzles. Injector ESC mode PSOI PPW GSOI GPW EQR [O2 based] EGR [%] S-H-S-A B-75 -23 0.62 -14.8 1.61 0.62 18 S-H-L-A -22.5 -14.4 1.45 L-H-S-A -22.5 -14.6 1.41 L-H-L-A -21.5 -13.6 1.38   56  Table 13- Inputs for CFD simulation of paired-nozzle injectors under mode B-75.  Pilot Gas Charge state at -90 aTDC SOI PPW τ SOI GPW τ P [bar] T [K] YO2 YN2 YCO2 YH2O S-H-S-A -18 0.84 0.56 -8.84 2.48 0.67 4.66 430 0.205 0.7600 0.0194 0.0156 S-H-L-A -17.5 0.84 0.67 -8.85 2.37 0.62 4.48 429 L-H-S-A -19.5 0.84 0.33 -8.57 2.12 0.67 4.60 434 L-H-L-A -18.5 0.84 0.33 -8.55 2.12 0.56 4.67 434        Figure 25- The rate of heat release curves for HPDI paired-nozzle injectors. Small holes (top), large holes (bottom), small angles (left), large angles (right) -60 -40 -20 0 20 40 60 80-50051001502002501112.5CA []HRR [kJ/(m3.deg)]  SCREGOLDS-H-S-A-60 -40 -20 0 20 40 60 80-5005010015020025011.112CA []HRR [kJ/(m3.deg)]  SCREGOLDS-H-L-A-60 -40 -20 0 20 40 60 80-5005011502002501111CA []HRR [kJ/(m3.deg)]  SCREGOLDL-H-S-A-60 -40 -20 0 20 40 60 80-50050015020025011.412.5CA []HRR [kJ/(m3.deg)]  SCREGOLDL-H-L-A57  The pilot and gas actuation delays are chosen to yield similar (within 1°) pilot and gas ignition times, as well as the 50% IHR, shown in Figure 25. GOLD predicts a much stronger premixed peak than observed in measurements for the L-H-S-A injector. In contrast, for the S-H-L-A injector, GOLD underestimates the peak heat release. Given that the injection timing is matched with the experimental HRR, more premixing stems from delayed ignition of gas, that is, more accumulated fuel before combustion. After careful investigation of the injectors under a microscope by Chris Mabson, it was noticed the two injectors with mismatched AHRR curves do not have aligned gas and diesel nozzles, Figure 26 [48]. Given this observation, the two S-H-L-A and L-H-S-A injectors were simulated with an 8˚ gas-diesel interlace offset. Expectedly, changing the interlace angle has removed the large premixed spike at the beginning of the L-H-S-A combustion and increased the premixing in S-H-L-A injector and resulted in better matching of peaks in AHRR plots, Figure 27. The current results indicate that gas-diesel interlace angle has a significant effect in defining the shape of AHRR and the level of premixing for these configuration of nozzles.     Figure 26- Misalignment between the gas and diesel holes in L-H-S-A injector. Figure adopted from Mabson [48]. 58      The L-H-S-A and S-H-S-A resulted in the lowest and highest PM emissions, respectively, and in section 3.2.2 onward, comparisons focus on these two injectors. The SCRE measurements indicate that even the best paired-nozzle injector emit more PM than the baseline single-hole injector, 7-7-0 (Figure 28).     The qualitative agreement of the predictions with measurements is encouraging, given the very simple soot model used here, and the uncertainty in the LES predictions due to cyclic variability and mesh sensitivity. In the following sections, the jet dynamics and the fuel mixture are studied using the penetration and entrainment plots, mixture evolution, phi-T maps and EQR and mixture fraction distributions. The CFD-predicted paired-nozzle emissions are added to the PM-NOx and PM-CO emission maps presented in section 3.1.1. No simulations were conducted to study the EGR sweep for paired-nozzles, so only the baseline B-75 points are presented. The results lie on the very right-hand--60 -40 -20 0 20 40 60 80-5005010015020025011.112CA []HRR [kJ/(m3.deg)]  SCREGOLDS-H-L-A-60 -40 -20 0 20 40 60 80-500501001502002501110.5CA []HRR [kJ/(m3.deg)]  SCREGOLDL-H-S-AFigure 27- HRR graphs of repeated S-H-L-A (left) and L-H-S-A (right) simulations with adjusted interlace angle. Figure 28- Measured (red) and simulated (blue) NOx (left) and PM (right) emissions for 7-7-0 and paired-nozzle injectors. 7-7-0 L-H-S-A S-H-S-A S-H-L-A L-H-L-A00.511.52NOx [gr/kW.hr]  GOLDSCRE7-7-0 L-H- -A S-H- -A S-H- A L-H-L-A00.050.10.150.2PM [gr/kW.hr]  GOLDSCRE, DRX59  side of the graphs showing generally higher levels of PM for all of these nozzles. The paired-jet results also closely scatter around the extrapolation of the PM-CO correlation generated for baseline injectors.       3.2.2 Penetration, air entrainment, and mixture formation Jets with higher penetration rates tend to utilize the oxygen more effectively and form less rich mixture, hence less soot [13,27]. The penetration of the single jet and the paired jets are compared on the jet centreline plane as discussed in section 3.1.2. As expected, the penetration rate of the single-nozzle injector exceeds the paired-nozzle both before and after impingement, Figure 30. Upon impingement, the jets suddenly spread on the piston, showing an abrupt increase of the nominal penetration lengths. The L-H-S-A and S-H-S-A jets show very close resemblance through the course of measurement, which is the instance of merging of the neighbouring jets. This happens at almost 25% IHR for the 7-7-0 injector and slightly before 50% IHR for the paired-nozzle injectors.  0 0.05 0.1 0.15 0.2 0.2502468PM [gr/kW.hr]NOx [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AGOLD0 0.05 0.1 0.15 0.2 0.2502468PM [gr/kW.hr]NOx [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-ASCRE10-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AGOLDCO~PM0.1310-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J-367-7-0L-H-S-A-H-S-ASCRECO ~ PM0.57Figure 29- Paired-nozzle emissions added to the baseline PM-NOx (top) and PM-CO (bottom) maps. Trend-lines are only fit to the baseline data. 60      Figure 30 also shows the evolution of Zair parameter (definition in section 0). The relative rate of air entrainment is similar for all three injectors, especially before the impingement point; this is counter-intuitive as it is expected that the paired-hole injectors entrain more air into the base of the jet. However, if the Zair values at the same penetration distances are compared, it is seen that paired-hole injectors do entrain more air (marked with dashed arrows in Figure 30). This suggests the plots of Zair versus the penetration distance should be examined. Figure 30-right also shows that the paired-hole jets entrain twice as much air at the same distance. This means that the mode of mixing is different between the injectors; while the 7-7-0 jet mixes with cylinder charge and produces the combustible mixture mostly through penetrating mechanism, the paired-hole injectors entrain in short distances and spread more radially. Figure 31 shows the evolution of the cylinder soot (left) and Zrich-Zlean parameters indicating the normalized amount of sooting and lean mixtures in the cylinder (right). From Figure 31-right, the two paired-hole injectors have different routes in producing more soot. The L-H-S-A injector forms more rich mixture early after start of injection because of its larger nozzle diameter and entrains less air from the periphery; while 7-7-0 and S-H-S-A injectors form almost equal levels of rich mixture. Nevertheless, the Zrich profile indicates that S-H-S-A suffers from poor late mixing which is also evident from the evolution of the cylinder lean mixture, Zlean, showing S-H-S-A injector has 10-15% less lean mixture compared to 7-7-0 and L-H-S-A. -10 -5 0 5 102060CA [deg]penetration [mm]-5 0 5 100510Zair  7-7-0L-H-S-AS-H-S-AFigure 30- Penetration and entrainment plots (left) for 7-7-0 (blue), L-H-S-A (green, x), and S-H-S-A (red, +). The two graphs on the left are combined into an entrainment-penetration plot (right). 0 10 20 30 40 50012345Penetration [mm]Zair  7-7-0L-H-S-AS-H-S-AImpingement61      The integrated effect of penetration and mixture evolution can be observed in the PM-crank angle graphs. According to the Hiroyasu soot model, Eq.24 section 2.2.4, the following factors contribute to higher soot formations  1- More rich fuel in the high pressure and temperature regions. 2- Longer residence time of soot in the high temperature regions. For paired-hole injectors both factors promote more soot formation; there is more rich mixture to begin with, increasing soot formation through mechanism 1. Also the jets penetrate less, therefore staying longer in the soot forming regions.  The slope and peak value of the formed soot is an indication of poor early-stage mixing for the large-hole nozzle. However, during the oxidation reactions, the L-H-S-A injector demonstrates close similarity to the single-hole nozzles; whereas the S-H-S-A injector shows a much slower oxidation rate later in the cycle. To find out the main reason for the poor performance of small-hole nozzles during late mixing, further in-cylinder investigation is needed. In the following section, the phi-T maps and EQR distribution plots are used for more in-depth study of these phenomena.   Figure 31- Evolution of soot (left) and rich and lean mixtures (right) in the cylinder. Comparison is made between 7-7-0, L-H-S-A, and S-H-S-A injectors. -50 0 50 1000123CA []PM [g/kg fuel]  7-7-0L-H-S-AS-H-S-A-10 10 20 3000.20.40.6CA [deg]Zrich,    Zlean  7-7-0L-H-S-AS-H-S-AZleanZrich62  3.2.3 Equivalence ratio - temperature maps The evolution of mixture formation is mapped on φ-T planes in Figure 32 for the 7-7-0 nozzle in three different timings, 25%, 50% and 75% of heat release. Details about the generation and interpretation of φ-T maps could be found in section 2.5.2. The CFD results are averaged over 1.5° crank angle (each crank angle averaged with its previous and next time steps), to smear the LES that might interfere with the interpretation of the graphs. Figure 32 shows that only a small part of the mixture is in the PM-forming region, even early in the process. At 75%-IHR, none of the mixture is in the PM-forming region, and only a little in the CO-forming zone. For a more direct comparison of each two injectors, the difference between fuel mixtures, instead, are shown on the phi-T planes in Figure 33. Dark red color indicates more mass for the first injector in a certain region and dark blue indicates otherwise. The difference of mixture mass between the L-H-S-A and 7-7-0 injectors at 25%, 50%, and 75% of IHR is shown in the top row of Figure 33. The sequence of figures strongly implies that the L-H-S-A injector produces rich mixture, especially during the early combustion event, confirming higher rate of soot formation early after start of injection. After 50% progress of the combustion, the difference between the two injectors fade away; however, there is still more fuel mass in the slightly rich mixture (1<EQR<2), resulting in more formation of CO. In the late combustion timing, 75% of IHR and afterwards, the two injectors barely show any difference. Even though it produces more rich mixture compared to the baseline injector, comparing S-H-S-A with 7-7-0 injector demonstrates more similarity at early timings. Nevertheless, the figures show that after the mid combustion event, the distribution of fuel for S-H-S-A injector is shifted to higher equivalence ratios and continues to remain the same until the end of combustion. This analysis shows that the higher PM emissions from the two paired jets originate from different periods of the combustion. While the large-hole injector forms a richer mixture in the beginning of the cycle and potentially forms more soot, it oxidizes the formed soot with the same rate as the 7-7-0 injector. On the other hand, the small nozzle forms less rich mixture at the beginning, instead has a very poor late mixing characteristics; therefore oxidizing less soot compared to the other two injectors. 63       Figure 32- Status of the mixture on the φ-T planes for the 7-7-0 injector at 25% (left), 50% (middle), and 75% (right) of IHR. 64          3.2.4 Mixture probability density function The plots of the mixture EQR distribution are presented in Figure 34 for a more quantitative study of the fuel distribution at early and mid-combustion timings. Paired-nozzle jets demonstrate similar pattern of EQR distribution in early injection timings, ramping up in near-stoichiometric regions and then ramping down with the same slope followed by a sudden decrease at EQRs between 2.5 to 3.5. Figure 34 reinforces the observations from the φ-T maps; i.e. the small nozzles form leaner mixture in the early-stage mixing compared to the large nozzles. Nevertheless, the paired-nozzles produce very similar distributions in the EQR range of 1.5-2.5. During mid-combustion, on the other hand, the small-hole nozzle contains significant amount of rich mixture in the EQR range of 2.5-4, corresponding to the temperatures of 1800-2000 K ideal Figure 33- Plots of fuel mass difference between L-H-S-A (left), and S-H-S-A (right) with the baseline 7-7-0 injector at the times of 25% (top), 50% (middle) and 75% (bottom) of IHR. 65  for fast soot formation. Observation of the pdf patterns for early and mid-combustion timing reveals the shifting of the early rich zone (φ between 1.5-3) to the leaner side. This indicates that the rich mixture formed in the early timings is transferred to the lean side as the jet is penetrating farther and mixing with more air. Nonetheless, there are dramatic differences between the two injectors in the very rich side of the plot (φ between 2.5-4.5).     This rich mixture for the S-H-S-A injector originates from the core of the jet, where the fuel stream has not yet mixed with enough oxidizer. Furthermore, this mixture is produced while the gas needle is closing the nozzle (between 8˚ to 14˚ aTDC); hence reducing the injection momentum Figure 34- The distribution function of the equivalence ratio (top) for 7-7-0 and paried-nozzle injectors at 25% (left), and 50% (right) IHR timing. Different φ intervals are marked up on the injection plane (bottom). 66  and decelerating the gas flow. Therefore, the S-H-S-A injector mostly suffers from the poor end-of-injection mixing because of the low-momentum fuel parcels from the nozzle. 3.2.5 Jet interaction with cylinder and fire-deck High emissions could arise from the attachment of the jets to the fire-deck. This deteriorates the mixing and results in the build-up of rich mass near the head, eventually leading to higher rates of soot. The fire-deck attachment is more severe for the paired-nozzle injectors, especially the small-nozzle group, due to the blockage of air access to the top part of the cylinder, between the gas jet and the cylinder head. To investigate the effect of fire-deck attachment on the mixture formation of paired-nozzle injectors, CFD simulations with increased injector tip protrusion inside the cylinder were carried out. The measured and modeled gas nozzle protrusions inside the cylinder are summarized in Table 14. Figure 35 shows the shape of jets in CFD simulations with small (baseline) and large protrusions. After increasing the gas jet protrusion inside the cylinder, the fire-deck-attachment issue is mostly mitigated, but the targeting of the jet on the piston is changed. In the former case the jet is targeted to impinge on the piston lip; whereas in the latter it hits the piston bowl and do not get the extra mixing benefit after impingement. This phenomenon adversely affects the post-impingement mixing and pollutant formation; therefore, to avoid further complexities only the penetration and mixture evolution of the jets are compared to directly assess the effect of the gas nozzle protrusion on the jet dynamics. Table 14- The gap between the gas nozzle centreline and the cylinder head Injector Gas nozzle protrusion in the cylinder [mm] Measured Baseline CFD Increased protrusion CFD 7-7-0 1.2 1.2 3.2 L-H-S-A 1.6 1.2 3.2 S-H-S-A 1.5 1.2 3.2    67      As Figure 36 suggests, increasing the gap between the nozzle and the cylinder head does not significantly affect the pattern of mixture formation and the penetration of jets. The jets with higher protrusion (thinner dashed lines) show similar rate of penetration as the baseline cases and similar path of rich mixture formation. Given these results, the likelihood of the fire-deck interference with the mixture formation of the paired-jets is considered to be insignificant. Distillation of the results from previous sections suggests the sooting tendency of the paired-nozzles should be attributed to the jet dynamics characteristics of the jets, not the cylinder/charge motion interactions.       0 0.5 1 1.5 2 2.5 3020406080Time after start of injection [ms]Jet penetration [mm]  7-7-0L-H-S-AS-H-S-A7-7-0: Large protrusionL-H-S-A: Large protrusionSHSA: Large protursion-5 0 5 10 15CA aTDC [deg]Impingement-10 0 10 20 3000.050.10.150.20.250.3CA [deg]Zrich  7-7-0L-H-S-AS-H-S-A7-7-0, Large protrusionL-H-S-A, Large protrusionS-H-S-A. Large protrusionFigure 35- Fuel jet interaction with the cylinder wall for 1.2-mm (left) and 3.2-mm (right) gaps between nozzle centre-line and cylinder head. Figure 36- Penetration (left) and rich-mixture evolution (right) for small and large protrusion injections. 68  3.3 Analysis of the LPI strategy In SCRE tests, Late Post Injection (LPI) reduced the PM emissions by up to 75% without negatively affecting other emissions with only slight fuel penalties. Table 15 summarizes the results of the optimized LPI engine tests on the SCRE facility, borrowed from the work of Ehsan Faghani [49]. Table 15- Performance of the optimized LPI strategy, SCRE measurements and CFD.  CFD Experiments Baseline LPI Baseline LPI PM (g/kw-hr) 0.055 0.028 0.02 0.005 NOx (g/kw-hr) 1.08 0.99 1.4 1.4 CH4 (g/kw-hr) 0.43 0.44 0.55 0.45 GISFC (g/kw-hr) 171 175 175 176  Understanding the mechanism of soot reduction in LPI in HPID engines is pursued in this study. There is no unanimous agreement in the diesel literature about the effectiveness of this strategy. Enhanced second-pulse mixing, soot oxidation by the combustion of the second pulse, and the short-pulse effects are suggested as the effective mechanism in diesel literature [30]. SCRE tests shows that in an LPI strategy, the second pulse forms negligible amounts of soot [49]. This finding reinforces the third suggested mechanism as being the main cause of soot reduction in the SCRE facility. These hypotheses are tested with GOLD by simulating SCRE under 5 operating conditions, outlined in Table 16. The purpose of these tests is to show the evolution of the LPI strategy as the pulse separation is increased up to the optimum GSEP value and compare the results with the injection of single pulses with 100% and 85% of the nominal fuel. It is desired to see how the strategy improves as the flames are separated. 3.3.1 Heat release rates and final emissions First, Figure 37 shows the comparison the CFD heat release curves with their experimental counterparts to match the phasing of the three combustion events. The figures demonstrate accurate matching between the pilot and first gas combustion timings. The second gas jet shows a more diffuse and slower combustion, marked by shorter second peak and wider duration, 69  compared to the SCRE. Different combustion behaviour of the second pulse is due to the different state of the charge at the beginning of the second injection. A major portion of the cylinder O2 is already consumed; therefore the second pulse burns in high-EGR like charge resulting in a slower combustion. This behaviour is more accentuated in GOLD, because of the chemistry trajectories generated for only limited realizations of charge O2 concentrations. Table 16- SCRE and GOLD injection specifications for LPI tests. Number of injections 1st GSOI Fuel mass in the 1st pulse Fuel share of the 1st pulse GSEP 2nd GSOI Fuel mass in the 2nd pulse SCRE GOLD SCRE GOLD SCRE GOLD 1 -14.7 -6.7 167 100% - - - - - 1 -14.7 -6.7 141.5 100% - - - - - 2 -15 -6.8 142 85% 1.0 1.0 7 14.5 25 2 -15 -6.8 142 85% 1.5 1.6 11.5 20.5 25 2 -15 -6.8 142 85% 2.0 2.0 14.5 24.5 25       -60 -40 -20 0 20 40 60 80-5005010015020025011.211CA []HRR [kJ/(m3.deg)]  SCREGSEP = 1.0 ms-60 -40 -20 0 20 40 60 80-5005010015020025011.211CA []HRR [kJ/(m3.deg)]  SCREGSEP = 1.6 ms-60 -40 -20 0 20 40 60 80-50050105020025011.611CA []HRR [kJ/(m3.deg)]  SCREGSEP = 2.0 msFigure 37- Heat release rates for LPI simulations with GSEP’s of 1.0 ms (top-left), 1.5 ms (top-right), and 2.0 ms (middle-left). Plots are compared with experimental curves. 70  Figure 38 shows the difference between the emissions of the five simulated cases. Much higher predictions for the absolute simulated soot values are predicted; this pattern, also observed in the previous sections, could be mitigated by adjusting the model scaling factor to match at either of the 7-7-0 or J-36 baseline cases. The current scaling factor is unchanged however, to keep consistency in comparing various simulation cases.     LPI emissions are also plotted on the PM-NOx and PM-CO emissions maps alongside previous baseline and paired-hole emissions. LPI points lay in the same region in the PM-NOx maps as the SCRE data; i.e. in the same NOx level with reduced PM (left-shifted). Location of LPI points on the LPI strategy, on the other hand, is slightly above the baseline trend-line; while SCRE points lay below the trend-line, indicating more CO reductions for LPI strategy in reality.   Figure 38- NOx (left) and PM (right) emissions for single and late-post injections. SI(100%) SI(85%) GSEP=1.0ms GSEP=1.6ms GSEP=2.0ms00.511.52NOx [gr/kW.hr]  GOLDSCRESI(100%) SI(85%) GSEP=1.0ms GSEP=1.6ms GSEP=2.0ms00.020.040.060.080.1PM [gr/kW.hr]  GOLDSCRE, DRX71        3.3.2 PM formation in LPI strategy The evolution of the absolute mass of formed PM in the cycle is compared among the five cases in Figure 40. All the LPI cases form the same mass of soot as a single injected pulse with 85% of the nominal fuel. Therefore as the figure suggests, a main benefit in using an LPI strategy is the reduced mass of soot formed (the peak in the graph) due to the reduced injected mass in the pulse. Furthermore by increasing GSEP, the final level of soot emissions is reduced to a level very close to that of a single small pulse. This implies that the close-coupled pulses are not as effective as split pulses. In other words, close interaction between the pulses generates more soot in comparison with the isolated pulses, i.e. split flames. 0 0.05 0.1 0.15 0.2 0.2502468PM [gr/kW.hr]NOx [gr/kW.hr]  J367-7-0L-H-S-AS-H-S-AJ-36, LPIGOLD0 0.05 0.1 0.15 0.2 0.2502468PM [gr/kW.hr]NOx [gr/kW.hr]  J367-7-0L-H-S-AS-H-S-AJ-36, LPISCRE10-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J-367-7-0L-H-S-AS-H-S-AJ-36, LPICO=PM0.13GOLD10-310-210-110010-1100101102PM [gr/kW.hr]CO [gr/kW.hr]  J367-7-0L-H-S-AS-H-S-AJ-36, LPICO ~ PM0.57SCREFigure 39- PM-NOx (top) and PM-CO (bottom) emission maps. 72     The split-flame concept is demonstrated in Figure 41; in the right hand side, the contours of formed soot in the cylinder are shown (view from top), the region in the black box highlights the formed soot in the second pulse. As is evident from the figures, the soot formation-oxidation reactions from the two pulses heavily interact with each other causing higher soot yields by the end of the cycle. On the contrary, by separating the sooty reactions from the two flames, the high flame temperature of the second pulse do not accelerate the formation reactions in the first pulse. The difference between the formed soot from the first and second pulses are also marked in Figure 41, showing the effect of pulse interaction in soot formation. Comparison of LPI with other strategies on emission maps is conducted in the next section after discussing SPC strategy. 01234x 10-4In-cylinder PM [gr]  -80 -40 0 40 80 12001CA []Fuling [a.u]SI, 100% FISI, 85% FILPI, GSEP = 1.0 msLPI, GSEP = 1.6 msLPI, GSEP = 2.0 msFigure 40- Evolution of in-cylinder soot for single injection cases with nominal and 85% injected fuel, and three LPI cases with GSEP’s of 1.0, 1.5, and 2.0 ms. 73     3.4 Analysis of SPC injection strategy  Engine experiments on SCRE showed that Slightly Premixed Charge (SPC) combustion offers great reduction in soot with less sensitivity to the EGR and EQR levels. The challenges in simulating the SPC strategy with GOLD are discussed and the advantages in mixture formation and PM evolution are discussed. 3.4.1 Heat release rates and final emissions The SCRE and GOLD inputs for SPC operating conditions are summarized in Table 17 which are attempts to simulate the SPC regime at negative PSEPs of -0.9 and -2.1 ms. Figure 42 shows the simulated HRR curves compared with the experimental HRR for a PSEP of -0.9 ms. For these simulations the standard J36 injection delay for previous cases were used to maintain a consistent conversion of SCRE to GOLD data. Different parameters are changed to delay the combustion timing of the gas jet but none was successful. The main difference between the simulated combustion behaviour versus the experimental trends is the onset of gas ignition. GOLD predicts instantaneous gas ignition upon the ignition of diesel fuel; whereas, the experimental HRR suggests a transition phase exists with low heat-release rates prior to main gaseous combustion. Figure 41- Regions of high PM concentration in a close-coupled (left) versus split-flame (right) regimes. 74  Table 17- SCRE operating parameters in SPC strategy PSEP PSOI GSOI EGR EQR SCRE GOLD SCRE GOLD 0.3 (Baseline) -23 -17 -14.7 -7 18% 0.62 -0.9 -17 -11, -5 -19.6 -12, -12 18% 0.62 -2.1 -14 -11, -8 -27.4 -21.4, -19 18% 0.62  The issue of combustion phasing mismatch lies in the method of handling natural gas and diesel combustion in the CFD tool. Separate sets of trajectories handle the reaction progress of natural gas and diesel while the chemical effect of diesel-gas mixing is not considered. In other words, the diesel flame ignites the natural gas jet through temperature effects without considering the chemical interaction of the two fuels. This leads to problems in correctly capturing the ignition delay of the gas when operating in the SPC regime wherein considerable diesel-gas mixing happens prior to ignition. In order to compensate for this effect, changes in the gas-diesel interlace angle and injection temperature are made to delay the gas ignition through physical separation of the flames or temperature-induced chemical delay. In spite of these efforts, no significant changes are induced in the phasing of combustion pointing to the significance of diesel-gas interaction in this slightly premixed regime.    75     To simulate the SCRE combustion chamber to the best approximation, the diesel pulse is retarded to match the phasing of combustion. Increasing the negative gap between the gas and diesel pulses results in more gas premixing manifested as much higher HRR peaks in Figure 43. Therefore, the pilot pulse adjustment is made only to the degree to match the beginning of the gas combustion (start of the HRR ramp up).     3.4.2 Mixture formation The evolution of the rich mixture inside the cylinder is shown for the two simulated cases in Figure 44. SPC demonstrates two beneficial features. First, it shows significant increase in the level of the premixed fuel at the time of ignition; in particular, the premixed fraction in the case with PSEP=-0.9 ms is 35% and for PSEP=-2.1 is 57% compared to the 3% premixed fuel for the baseline case. Second, the profile of the rich mixture (2<EQR<6) has changed when switching to -60 -40 -20 0 20 40 60 80020040060080011.2CA []HRR [kJ/(m3.deg)]  SCREPSEP = -0.9 [ms], BaselinePSEP = -0.9 [ms], Adjusted timing-60 -40 -2 0 20 40 60 80020040060080011.3CA []HRR [kJ/(m3.deg)]  SCREPSEP = -2.1, BaselinePSEP = -2.1, adjusted timingFigure 42- Simulated (coloured) and SCRE (black) HRR curves for the standard J-36 actuation delay. Attempted methods to delay the gas ignition are listed in the legend box. Figure 43- Simulation of SCRE at SPC regime with varied PSOI for closest matching of HRR phasing. -60 -40 -20 0 20 40 60 800100200300400CA []HRR [kJ/(m3.deg)]  SCRESPC, PSEP = -0.9 [ms]Interlace + 20 [deg]Interlace - 20 [deg]Phenomenologicalgas-diesel chemistryTinj - 70 [K]76  the SPC strategy. The maximum height of the normalized rich mixture is reduced to 12% from a baseline value of 25%. The second difference, given that same injector was used and constant rail pressure was maintained, is counterintuitive in the first glance. However, by retarding the ignition during the gas injection cylinder pressure is kept low; hence the sac-to-chamber pressure ratio stays higher than the baseline case. As a result, the fuel is injected in a lower pressure, lower density charge; thus penetrates and mixes with more air. This causes less formed rich mixture and faster rate of lean mixture evolution.     Figure 45 shows the comparison between the jet penetrations in a J-36 baseline case compared to a J-36 jet operating at SPC regime. In-line with the observations from Figure 44, the jet penetrates faster under SPC regime than in the conventional baseline condition.   -10 -5 0 5 10 15 20 2500.20.40.60.81CA [deg]Zrich , Zlean  J36, PSEP = 0.3 [ms]J36, PSEP = -0.9 [ms]-20 -10 0 10 20 3000.20.40.60.81CA [deg]Zrich, Zlean  J36, PSEP = 0.3 [ms]J36, PSEP = -2.1 [ms]Figure 44- Evolution of the rich and lean mixtures for the simulated SPC cases with PSEP=-0.9 ms (left) and PSEP=-2.1 ms (right). 77     The non-premixed combustion at PSEP = 0.3 ms and SPS with PSEP = -2.1 ms are compared on φ-T maps in Figure 46. SPC combustion almost entirely bypasses the soot formation region by allowing far more premixing than the conventional non-premixed case. The NOx emissions are expected to be higher since, more mixture enters the NOx formation peninsula. The reason for such high-temperature stoichiometric flame is the advancement of injection in the PSEP = -2.1 ms and therefore an almost instantaneous combustion at TDC. Figure 47 shows the evolution of cylinder PM for the baseline case and the two SPC cases with PSEPs of -0.9 and -2.1 ms. In-line with previous observations, minute amounts of soot is formed, of which almost all is oxidized for PSEP = -2.1 ms.   Figure 45- fuel penetration under baseline and SPC J36 injection. 0 0.5 1 1.5 20102030405060Time after start of injection [ms]Jet penetration [mm]  7-7-0, B75J-36, B75J36, Neg-PSEP78        Figure 46- Mixture difference between positive (blue) and negative (red) PSEPs on Equivalence ratio-Temperature maps at 11° (top), 20°, and 23° (bottom) after start of injection (aSOI). 79      Figure 47- PM evolution for J36 injection at PSEP=0.3 ms (baseline-black), PSEP=-0.9 ms (blue), and PSEP=-2.1 ms (green). -50 0 50 1000123CA []PM [g/kg fuel]  J36, PSEP = 0.3 [ms]J36, PSEP = -2.1 [ms], Adjusted timingJ36, PSEP = -0.9 [ms], Adjusted timing80  4 Conclusions The performance of two injection strategies and novel paired-nozzles for HPDI combustion were investigated numerically. The numerical tool employed an LES turbulence model to resolve the turbulent charge and injection motion. Fuel oxidation chemistry was handled through pre-tabulated TGLDM method, NOx through an extended Zeldovich mechanism and soot with a phenomenological 2-equation model. To assess the validity of this LES study the following simulations were carried out: a) The sensitivity of the results to the CFD grid size was investigated. Reducing the mesh size was carried out, only for the duration of injection. Refining the mesh size resulted in the manifestation of more turbulence structure. J-36 showed the least sensitivity to the mesh size with less than 4% change in all the outputs after two stages of mesh refinement. In general, the mesh sensitivity was most extreme for soot predictions compared to other emissions. The highest mesh sensitivity was seen in 7-7-0 and fine-mesh L-H-S-A simulations; yet on average higher soot levels were reported for paired-nozzle injectors. b) Acceptable predictions for emissions trends were observed for input parameter sweep tests. Limited EGR sweep simulations reproduced experimental trends for NOx and PM. More tests were simulated for injection timing and pressure effects; overall trends were satisfactory, but fluctuations of up to 100% about the trend-line were observed. Mixture evolution and penetration were not substantially different, thus the difference in soot is mostly attributed to local phenomena, notably near-piston rich mixture. c) Emissions demonstrated trade-off behaviour when plotted in a PM-NOx space and showed a linear correlation when plotted in log scale on the PM-CO maps. The exponential dependence of the PM and CO was much weaker than the experimental observations. 1- Paired-nozzles The penetration of paired jets was less than a single jet. This is however, accompanied with higher rate of air entrainment during the injection. Therefore, the concept of early mixing is achieved to some degree; however, the effect of penetration is more crucial in mixture 81  preparations and ought not to be compromised. Consistent with the experiments, L-H-S-A and S-H-S-A were the best and worst paired-nozzle injectors. S-H-S-A injector particularly suffered from end-of-injection mixture leaning. More simulations demonstrated similar patterns even with larger gap between the nozzle and cylinder head, ruling out the effect of fire-deck attachment to such high soot formations. Mixture distribution plots and in-cylinder contours showed thicker rich jet core especially for the S-H-S-A nozzle. 2- LPI The benefit of post injection for natural gas injection arises from the flame splitting. Mass of formed soot from the first pulse was equal to an isolated single small injection (85% of fuel) which formed 20% less soot compared to the baseline pulse. However, after soot oxidation the amount of net soot from a single pulse was only 40% of the baseline case, demonstrating non-linear correlation between pulse width/mass and net soot. The net soot of LPI points were higher than the short single-pulse case, but the gap was less with higher GSEPs. Providing enough splitting for the injection events, the soot reactions do not interact; larger GSEPs are particularly effective reducing the formed soot from the second pulse as shown in the PM evolution graphs. 3- SPC Two benefits are gained when gas injection is advanced relative to the diesel injection. First, the level of gas premixing from 1-4% is increased to 30-50%; therefore less rich fuel is available for soot formation. Likewise, the ignition point is retarded in the Zrich plot, indicating the rich mixture is only available during a short period at the end of injection. Second, due to higher sac-to-cylinder pressure ratio during injection higher penetration rates and better air utilization is experienced during SPC operation. Phi-T maps of SPC showed that by allowing far more premixing in SPC, the entire PM region is bypassed, justifying the minute amount of formed soot.    82  5 Recommendations There are further questions regarding the simulation of HPDI using numerical tools, in particular: 1- Full-cylinder multi-cycle simulation is needed to go along with the LES model. Given the substantial sensitivity of soot to the mesh and simulation inputs, an averaging procedure is needed to obtain statistically converged results. This demands high computational times, but could be achieved, by parallelizing the simulation on multi-core machines. Cyclic engine emission studies could facilitate this, providing experimental inputs for bridging between cyclic fluctuations and multi-cycle mean measurements. 2- To study the soot evolution better, the distribution of the soot mass, instead of the fuel mass, should be studied on φ-T planes. Using this approach, the regions of high soot formation and oxidation could be visualized directly, instead of being inferred based on the mixture location relative to high C2H2 regions. 3- CO emissions showed less sensitivity to the HPDI operating conditions than soot. This could be a result of manifold parameterisation method in TGLDM, inaccuracy of current natural gas combustion mechanism, or tail-pipe oxidation. Optical studies could help identifying the in-cylinder CO formation phenomenology to back up potential combustion mechanism improvements and implementation in CFD. 4-  Paired-hole injection showed much difference from the baseline single-hole injection. To confirm this experimentally, optical assessment of mixture properties and penetration of these injectors should be undertaken. 5- Soot formation from the second pulse is minimized by providing enough gap for the post injection. Optical in-cylinder measurements or more advanced soot chemistry in CFD can provided more insight into this mechanism. 6- Gas-diesel kinetic interaction is blamed for the poor combustion predictions in SPC regime. To improve GOLD predictions for this regime, chemistry of natural gas mixed with trace diesel (heptane) should be modeled either kinetically or phenomenologically. 83  Bibliography [1] Davis, S. C., Diegel, S. W., and Boundy, R. G., 2014, "Transportation Energy Data Book," Oak Ridge National Laboratory, Oak Ridge [2] Crude Oil and Petroleum Products Market. (n.d.). Retrieved February 10, 2015, from www.nrcan.gc.ca/energy/crude-petroleum/4541. [3] Birol, F., 2009, "World Energy Outlook," International Energy Agency (IEA), Paris. [4] Andres, R. J., Boden, T. a., Bréon, F. M., Ciais, P., Davis, S., Erickson, D., Gregg, J. S., Jacobson, A., Marland, G., Miller, J., Oda, T., Olivier, J. G. J., Raupach, M. R., Rayner, P., and Treanton, K., 2012, “A Synthesis of Carbon Dioxide Emissions from Fossil-Fuel Combustion,” Biogeosciences Discuss., 9, pp. 1299–1376. [5] Lloyd, A. C., and Cackette, T. a., 2001, “Diesel Engines: Environmental Impact and Control,” J. Air & Waste Manage. Assoc., 51(6), pp. 809–847. [6] Sydbom, A., Blomberg, A., Parnia, S., Stenfors, N., Sandström, T., and Dahlén, S. E., 2001, “Health Effects of Diesel Exhaust Emissions.,” Eur. Respir. J.  Off. J. Eur. Soc. Clin. Respir. Physiol., 17, pp. 733–746. [7] Mauderly, J. L., 1994, “Toxicological and Epidemiological Evidence for Health Risks from Inhaled Engine Emissions,” Env. Heal. Perspect, 102(4), pp. 165–171. [8] Bünger, J., Krahl, J., Baum, K., Schröder, O., Müller, M., Westphal, G., Ruhnau, P., Schulz, T. G., and Hallier, E., 2000, “Cytotoxic and Mutagenic Effects, Particle Size and Concentration Analysis of Diesel Engine Emissions Using Biodiesel and Petrol Diesel as Fuel,” Arch. Toxicol., 74, pp. 490–498. [9] Benbrahim-Tallaa, L., Baan, R., Grosse, Y., Secretan-Lauby, B., El Ghissassi, F., Bouvard, V., Guha, N., Loomis, D., and Straif, K., 2012, “Carcinogenicity of Diesel-Engine and Gasoline-Engine Exhausts and Some Nitroarènes,” Pollut. Atmos., 2045(12), pp. 43–44. [10] Bernard, S. M., Samet, J. M., Grambsch, A., Ebi, K. L., and Romieu, I., 2001, “The Potential Impacts of Climate Variability and Change on Air Pollution-Related Health Effects in the United States,” Environ. Health Perspect., 109 (2), pp. 199–209. [11] DieselNet: Diesel Engine Emissions Online. Retreived January 20, 2015, from www.dieselnet.com/standards/eu/hd.php. [12] DieselNet: Diesel Engine Emissions Online. Retreived January 20 2015, from www.dieselnet.com/standards/cycles/#eu-hd. [13] McTaggart-Cowan, G., 2006, “Pollutant Formation in a Gaseous-Fuelled, Direct Injection Engine,” University of British Columbia, Vancouver. 84  [14] Korakianitis, T., Namasivayam, A. M., and Crookes, R. J., 2011, “Natural-Gas Fueled Spark-Ignition (SI) and Compression-Ignition (CI) Engine Performance and Emissions,” Prog. Energy Combust. Sci., 37(1), pp. 89–112. [15] Reitz, R. D., 2013, “Directions in Internal Combustion Engine Research,” Combust. Flame, 160(1), pp. 1–8. [16] Papagiannakis, R. G., Rakopoulos, C. D., Hountalas, D. T., and Rakopoulos, D. C., 2010, “Emission Characteristics of High Speed, Dual Fuel, Compression Ignition Engine Operating in a Wide Range of Natural Gas/Diesel Fuel Proportions,” Fuel, 89(7), pp. 1397–1406. [17] Reitz, R. D., and Duraisamy, G., 2014, “Review of High Efficiency and Clean Reactivity Controlled Compression Ignition (RCCI) Combustion in Internal Combustion Engines,” Prog. Energy Combust. Sci., 46, pp. 12–71. [18] Westport HPDI 2.0 (n.d.). Retreived December 20, 2014, from http://www.westport.com/is/core-technologies/hpdi-2. [19] Johnson, T., 2008, “Diesel Emission Control in Review,” SAE Int. J. Fuels Lubr., 1(1), pp. 68–81. [20] Johnson, T. V, 2009, “Review of Diesel Emissions and Control,” Int. J. Engine Res., 10(5), pp. 275–285. [21] Warnatz, J., Maas, U., and Dibble, R., 2001, "Combustion: physical and chemical fundamentals, modeling and simulation, experiments, pollutant formation," Springer, Berlin. [22] Khan, I., Greeves, G., and Wang, C., 1973, “Factors Affecting Smoke and Gaseous Emissions from Direct Injection Engines and a Method of Calculation,” SAE Tech. Pap., No. 730169. [23] Tree, D. R., and Svensson, K. I., 2007, “Soot Processes in Compression Ignition Engines,” Prog. Energy Combust. Sci., 33(3), pp. 272–309. [24] Mohan, B., Yang, W., and Chou, S. K., 2013, “Fuel Injection Strategies for Performance Improvement and Emissions Reduction in Compression Ignition Engines—A Review,” Renew. Sustain. Energy Rev., 28, pp. 664–676. [25] Gao, J., Matsumoto, Y., and Nishida, K., 2009, “Experimental Study on Spray and Mixture Properties of the Group-Hole Nozzle for Direct-Injection Diesel Engines, Part I: a Comparative Analysis with the Single-Hole Nozzle,” At. Sprays, 19(4), pp. 321–337. [26] Nishida, K., Gao, J., and Matsumoto, Y., 2009, “Experimental Study on Spray and Mixture Properties of the Group-Hole Nozzle for Direct-Injection Diesel Engines, Part Ii: Effects of Included Angle and Interval Between Orifices,” At. Sprays, 19(4), pp. 339–355. [27] Gao, J., Matsumoto, Y., Namba, M., and Nishida, K., 2009, “An Investigation of Mixture Formation and In-Cylinder Combustion Processes in Direct Injection Diesel Engines Using Group-Hole Nozzles,” Int. J. Engine Res., 10(1), pp. 27–44. 85  [28] Gao, J., Park, S. W., Wang, Y., Reitz, R. D., Moon, S., and Nishida, K., 2010, “Simulation and Analysis of Group-Hole Nozzle Sprays Using a Gas Jet Superposition Model,” Fuel, 89(12), pp. 3758–3772. [29] Zhang, L., 1999, “A Study of Pilot Injection in a DI Diesel Engine,” SAE Tech. Pap., No. 1999-01-3493. [30] O’Connor, J., and Musculus, M., 2013, “Post Injections for Soot Reduction in Diesel Engines: A Review of Current Understanding,” SAE Tech. Pap., No. 2013-01-0917. [31] Han, Z., Uludogan, A., Hampson, G., and Reitz, R., 1996, “Mechanism of Soot and NOx Emission Reduction Using Multiple-Injection in a Diesel Engine,” SAE Tech. Pap., No. 960633. [32] Faghani, E., Patychuk, B., McTaggart-Cowan, G., and Rogak, S., 2013, “Soot Emission Reduction from Post Injection Strategies in a High Pressure Direct-Injection Natural Gas Engine,” SAE Tech. Pap., No. 2013-24-0114. [33] Stanglmaier, R. H., and Roberts, C. E., 1999, “Homogeneous Charge Compression Ignition (HCCI): Benefits, Compromises, and Future Engine Applications,” SAE Tech. Pap., No. 1999-01-3682. [34] Gan, S., Ng, H. K., and Pang, K. M., 2011, “Homogeneous Charge Compression Ignition (HCCI) Combustion: Implementation and Effects on Pollutants in Direct Injection Diesel Engines,” Appl. Energy, 88(3), pp. 559–567. [35] Musculus, M. P. B., Miles, P. C., and Pickett, L. M., 2013, "Conceptual Models for Partially Premixed Low-Temperature Diesel Combustion," Prog. Energy Combust. Sci., 39(2), pp. 246-283. [36] Zhang, F., Xu, H., Zhang, J., Tian, G., and Kalghatgi, G., 2011, “Investigation into Light Duty Dieseline Fuelled Partially-Premixed Compression Ignition Engine,” SAE Tech. Pap., No. 2011-01-1411. [37] Munshi, S., McTaggart-Cowan, G., Huang, J., and Hill, P. G., 2011, “Development of a Partially-Premixed Combustion Strategy for a Low-Emission, Direct Inject High Efficiency Natural Gas Engine,” Proceedings of the ASME ICE Division Fall Technical Conference, ICEF2011-60181, pp. 1–14. [38] McTaggart-Cowan, G. P., Mann, K., Wu, N., Huang, J., and Munshi, S. R., 2012, “Particulate Matter Reduction from a Pilot-Ignited, Direct Injection of Natural Gas Engine,” Proceedings of the ASME Internal Combustion Engine Division’s 2012 Fall Technical Conference, ICEF2012, pp. 1–11. [39] Schaller, R. R., 1997, “Moore’s Law: Past, Present and Future,” IEEE Spectr., 34(6), pp. 52-59. [40] Thompson, S. E., and Parthasarathy, S., 2006, “Moore’s Law: The Future of Si Microelectronics,” Mater. Today, 9(6), pp. 20–25. 86  [41] Davis, R. S., Mandrusiak, G. D., and Landenfeld, T., 2009, “Development of the Combustion System for General Motors’ 3.6 L DOHC 4V V6 Engine with Direct Injection,” SAE Tech. Pap., No. 2008-01-0132. [42] Shrivastava, R., Hessel, R., and Reitz, R. D., 2002, “CFD Optimization of DI Diesel Engine Performance and Emissions Using Variable Intake Valve Actuation with Boost Pressure , EGR and Multiple Injections,” SAE Tech. Pap., No. 2002-01-0959. [43] Akihama, K., Takatori, Y., Inagaki, K., Dean, A. M., and Sasaki, S., 2001, “Mechanism of the smokeless Rich Diesel Combustion by Reducing Temperature,” SAE Tech. Pap., No. 2001–01–0655. [44] Kaario, O., Brink, A., Lehto, K., Keskinen, K., and Larmi, M., 2011, “Studying Local Conditions in a Heavy-Duty Diesel Engine by Creating Phi-T Maps,” SAE Tech. Pap., No. 2011-01-0819. [45] Rutland, C. J., 2011, “Large-Eddy Simulations for Internal Combustion Engines - a Review,” Int. J. Engine Res., 12, pp. 421–451. [46] Goryntsev, D., Sadiki, A., Klein, M., and Janicka, J., 2009, “Large Eddy Simulation Based Analysis of the Effects of Cycle-to-Cycle Variations on Air-Fuel Mixing in Realistic DISI IC-Engines,” Proc. Combust. Inst., 32(2), pp. 2759–2766. [47] Patychuk, B. D., 2013, “Particulate Matter Emission Characterization From a Natural Gas High-Pressure Direct-Injection Engine,” University of British Columbia, Vancouver. [48] Mabson, C., 2015, “Emissions Characterisation of Paired Gaseous Jets in a Pilot-Ignited Natural-Gas Compression-Ignition Engine,” University of British Columbia, Vancouver. [49] Faghani, E., 2015, “Effect of Injection Strategies on Particulate Matter Emissions in HPDI Natural-gas Engines,” University of British Columbia, Vancouver. [50] Faghani, E., and Rogak, S. N., 2012, “Application of CFD and Phenomenological Models in Studying Interaction of Two Turbulent Plane Jets,” Int. J. Mech. Eng. Mechatronics, 1(1). pp. 36-49. [51] Jasak, H., Jemcov, A., and Tukovic, Z., 2007, “OpenFOAM : A C ++ Library for Complex Physics Simulations,” International Workshop on Coupled Methods in Numerical Dynamics, IUC, Dubrovnik, pp. 1–20. [52] Huang, J., 2006, “Natural Gas Combustion under Engine-Relevant Conditions,” University of British Columbia, Vancouver. [53] Pope, S. B., 2000, “Turbulent Flows,” Cambridge university press, Cambridge. [54] Sone, K., and Menon, S., 2003, “Effect of Subgrid Modeling on the In-Cylinder Unsteady Mixing Process in a Direct Injection Engine,” J. Eng. Gas Turbines Power, 125(2), p. 435-443. 87  [55] Klimenko, A. Y., and Bilger, R. W., 1999, “Conditional Moment Closure for Turbulent Combustion,” Prog. Energy Combust. Sci., 25(6), pp. 595–687. [56] Bushe, W. K., and Steiner, H., 1999, “Conditional Moment Closure for Large Eddy Simulation of Nonpremixed Turbulent Reacting Flows,” 11(7), pp. 1896–1906. [57] Steiner, H., and Bushe, W. K., 2001, “Large Eddy Simulation of a Turbulent Reacting Jet with Conditional Source-Term Estimation,” Phys. Fluids, 13(3), p. 754-769. [58] Spalding, D. B., 1961, “A Single Formula for the "Law of the Wall",” J. Appl. Mech., 28(3), p. 455-458. [59] DeVilliers, E., 2006, “The Potential of Large Eddy Simulation for the Modeling of Wall Bounded Flows,” Imperial College of Science, Technology and Medicine. [60] Beale, J. C., and Reitz, R. D., 1999, “Modeling Spray Atomization with the Kelvin-Helmholtz/Rayleigh-Taylor Hybrid Model,” Atom. Spray, 9(6), pp. 623–650. [61] OpenFOAM User Guide | CFD Direct., Retreived February 4, 2015, from https://cfd.direct/openfoam/user-guide/. [62] Li, Y., Kirkpatrick, A., Mitchell, C., and Willson, B., 2004, “Characteristic and Computational Fluid Dynamics Modeling of High-Pressure Gas Jet Injection,” J. Eng. Gas Turbines Power, 126(1), pp. 192-197. [63] Ouellette, P., and Hill, P., 2000, “Turbulent Transient Gas Injections,” J. Fluids Eng., 122(4), pp. 743–752. [64] Foss, M. M., 2007, "Introduction to LNG: An Overwiew on Liquified Natural Gas (LNG), Its Properties, Organization of the LNG Industry and Safety Considerations", Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin, Huston. [65] Lifshitz, A., Scheller, K., and Burcat, A., 1971, “Shock-Tube Investigation of Ignition in Methane-Oxygen-Argon Mixtures,” Combust. Flame, 16(3), pp. 311–321. [66] Westbrook, C., 1979, “An Analytical Study of the Shock Tube Ignition of Mixtures of Methane and Ethane,” Combust. Sci. Technol., 20(1-2), pp. 5–17. [67] Asaba, T., Yoneda, K., Kakihara, N., and Hikita, T., 1963, “A Shock Tube Study of Ignition of Methane-Oxygen Mixtures,” Symp. Combust., 9(1), pp. 193–200. [68] Yamamoto, T., Kobayashi, N., Arai, N., and Tanaka, T., 1997, “Effects of Pressure on Fuel-Rich Combustion of Methane-Air under High Pressure,” Energy Convers. Manag., 38(10), pp. 1093–1100. 88  [69] Huang, J., Hill, P. G., Bushe, W. K., and Munshi, S. R., 2004, “Shock-Tube Study of Methane Ignition under Engine-Relevant Conditions: Experiments and Modeling,” Combust. Flame, 136(1-2), pp. 25–42. [70] Huang, J., and Bushe, W. K., 2006, “Experimental and Kinetic Study of Autoignition in Methane/Ethane/Air and Methane/Propane/Air Mixtures under Engine-Relevant Conditions,” Combust. Flame, 144(1-2), pp. 74–88. [71] Curran, H. J., Gaffuri, P., Pitz, W. J., and Westbrook, C. K., 1998, “A Comprehensive Modeling Study of n-Heptane Oxidation,” Combust. Flame, 114(97), pp. 149–177. [72] Mallampalli, H. P., Fletcher, T. H., and Chen, J. Y., 1998, “Evaluation of CH4/NOx Reduced Mechanisms Used for Modeling Lean Premixed Turbulent Combustion of Natural Gas,” J. Eng. Gas Turbines Power, 120(4), p. 703-712. [73] Maas, U., and Pope, S. B., 1992, “Simplifying Chemical Kinetics: Intrinsic Low-Dimensional Manifolds in Composition Space,” Combust. Flame, 88(3-4), pp. 239–264. [74] Pope, S. B., and Maas, U., 1993, “Simplifying Chemical Kinetics: Trajectory-Generated Low-Dimensional Manifolds,” Mechanical and Aerospace Engineering Report: FDA 11. [75] Frenklach, M., 2002, “Reaction Mechanism of Soot Formation in Flames,” Phys. Chem. Chem. Phys., 4(11), pp. 2028–2037. [76] Leung, K. M., Lindstedt, R. P., and Jones, W. P., 1991, “A Simplified Reaction Mechanism for Soot Formation in Nonpremixed Flames,” Combust. Flame, 87(3-4), pp. 289–305. [77] Tao, F., Golovitchev, V. I., and Chomiak, J., 2004, “A Phenomenological Model for the Prediction of Soot Formation in Diesel Spray Combustion,” Combust. Flame, 136(3), pp. 270–282. [78] Nagle, J., and Strickland-Constable, R. F., 1962, “Oxidation of Carbon between 1000-2000 C,” Proc. fifth carbon Conf., 1, New York: Pergamon. [79] Neoh, K. G., Howard, J. B., and Sarofim, A. F., 1981, “Soot Oxidation in Flames,” Particulate Carbon, pp. 162–282. [80] Hong, S., Wooldridge, M. S., Im, H. G., Assanis, D. N., and Pitsch, H., 2005, “Development and Application of a Comprehensive Soot Model for 3D CFD Reacting Flow Studies in a Diesel Engine,” Combust. Flame, 143(1), pp. 11–26. [81] Xi, J., and Zhong, B. J., 2006, “Soot in Diesel Combustion Systems,” Chem. Eng. Technol., 29(6), pp. 665–673. [82] Nishida, K., and Hiroyasu, H., 1989, “Simplified Three-Dimensional Modeling of Mixture Formation and Combustion in a D.I. Diesel Engine,” SAE Tech. Pap., No. 890269. 89  [83] Hiroyasu, H., and Toshikazu, K., 1983, “Development and Use of a Spray Combustion Modeling to Predict Diesel Engine Efficiency and Pollutant Emissions,” Bull. JSME, 26(214). [84] Li, Y. H., and Kong, S. C., 2008, “Diesel Combustion Modelling Using LES Turbulence Model with Detailed Chemistry,” Combust. Theory Model., 12(2), pp. 205–219. [85] Boussouara, K., and Kadja, M., 2009, “Numerical Investigation of Soot Formation in Diesel Jet Flame with KIVA-3V,” Revue Eng. Renouv., 12(1), pp. 55 – 62. [86] Heywood, J. B., 1988, "Internal combustion engine fundamentals," Mcgraw-hill, New York. [87] Celik, I. B., Cehreli, Z. N., and Yavuz, I., 2005, “Index of Resolution Quality for Large Eddy Simulations,” J. Fluids Eng., 127(5), pp. 949-958.     90  Appendices Appendix A: Initial engine mixture calculation Cylinder charge composition is calculated based on the overall fuel-oxygen equivalence ratio and the EGR flow rate. The flowchart for the calculation of mixture composition is shown in Figure 48.  Figure 48- Flowchart for cylinder (EGR-diluted) charge composition  For simplification, CH4 is considered as the only fuel constituent; also, only the major combustion products (CO2, H2O, O2, N2) are considered for the exhaust composition calculations, Eq.33. 𝐶𝐻4 +𝑎𝑠𝜙(𝑂2 + 3.76𝑁2) → 𝛼𝐶𝑂2 + 𝛽𝐻2𝑂 + 𝛾𝑂2 + 𝜃𝑁2  33  The stoichiometric coefficient, as, is calculated by substituting 1 for φ; for a given φ, the exhaust composition is calculated by balancing the mass of four atoms involved in the reaction, C, H, O, N. Assuming a standard 23% O2, 77% N2 as the constituents of the fresh air, the composition of the diluted charge is calculated using Eq.34. In this equation and Figure 48, the subscripts d, e, 91  and I stand for diluted, exhaust, and intake. Eq.34 is applicable to any generic species α in the mixture. 𝑦𝑑,𝛼 = (𝐸𝐺𝑅) × 𝑦𝑒,𝛼 + (1 − 𝐸𝐺𝑅) × 𝑦𝑖,𝛼  34 Using above calculations, the mixture composition for the three EGR levels used in this study are summarized in Table 18. The values of calculated oxygen mass fraction are compared with their measured counterparts at the intake manifold. The current calculations over-predict the initial oxygen mass fraction, but the error is less than 0.6% in all cases. This source of this error is neglecting the combustion products brought by EGR in the LHS of Eq.33. To take the effect of EGR products in the combustion reaction, the calculations should be done in an iterative manner which is not carried in this work given the negligible amount of error. Table 18- Sample initial charge composition for different working conditions. EGR EQR O2 N2 CO2 H2O Calculated Measured 0% 0.62 0.23 0.23 0.77 0.0 0.0 18% 0.62 0.2039 0.2026 0.7646 0.0173 0.0142 25% 0.62 0.1937 0.1922 0.7625 0.0240 0.0198     92  Appendix B: Φ-T map calculations Constructing φ-T maps includes two steps; firstly, a background emission map is generated using CANTERA with GRI-3.0 chemical mechanism. Products’ concentrations are calculated for 3-ms of reaction time at each equivalence ratio and temperature; emission contours are then generated given the final concentrations of C2H2, CO, and NOx. To maintain a constant temperature during the reaction, a large heat transfer coefficient is assumed, such that the reacting system is in equilibrium with the surrounding during the time step integration. To map the cylinder mixture on the emission contours, the CFD mixture is divided into φ and T bins. Figure 49 shows the flowchart for generating φ-T maps. The mixture equivalence ratio at each CFD cell is calculated given the mixture fraction value, Eq.35. The stoichiometric air-fuel ratio is calculated based on the complete combustion of methane with air, Eq.36; the subsequent value for the AFst is 17.16.     𝜙 =𝐴𝐹𝑠𝑡1𝜁−1  35 Figure 49- Flowchart for binning the mixture on a φ-T plane. 93  𝐶𝐻4 + 2(𝑂2 + 3.76𝑁2) → 𝐶𝑂2 + 2𝐻2𝑂 + 7.52𝑁2  36  The binned mixture composition is normalized with φ and T bin sizes, to rule out the dependence on the bins’ size or numbers. This is especially necessary when mapping the difference of the mixture between two cases on the φ-T maps. 

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:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0166233/manifest

Comment

Related Items