UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Combustion modelling in spark-ignition engines using conditional source-term estimation Nivarti, Girish Venkata 2013

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

Item Metadata


24-ubc_2013_fall_nivarti_girish.pdf [ 1.99MB ]
JSON: 24-1.0074086.json
JSON-LD: 24-1.0074086-ld.json
RDF/XML (Pretty): 24-1.0074086-rdf.xml
RDF/JSON: 24-1.0074086-rdf.json
Turtle: 24-1.0074086-turtle.txt
N-Triples: 24-1.0074086-rdf-ntriples.txt
Original Record: 24-1.0074086-source.json
Full Text

Full Text

Combustion Modelling in Spark-Ignition Engines usingConditional Source-term EstimationbyGirish Venkata NivartiB. Tech. (Honours), Indian Institute of Technology Kharagpur, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2013c? Girish Venkata Nivarti, 2013AbstractConditional Source-term Estimation (CSE) is a chemical closure model for thesimulation of turbulent combustion. In this work, CSE has been explored for mod-elling combustion phenomena in a spark-ignition (SI) engine. In the arbitrarilycomplex geometries imposed by industrial design, estimation of conditionally av-eraged scalars is challenging. The key underlying requirement of CSE is that con-ditionally averaged scalars be calculated within spatially localized sub-domains.A domain partitioning algorithm based on space-filling curves has been developedto construct localized ensembles of points necessary to retain the validity of CSE.Algorithms have been developed to evenly distribute points to the maximum extentpossible while maintaining spatial locality. A metric has been defined to estimaterelative inter-partition contact as an indicator of communication in parallel com-puting architectures. Domain partitioning tests conducted on relevant geometrieshighlight the performance of the method as an unsupervised and computationallyinexpensive domain partitioning tool.In addition to involving complex geometries, SI engines pose the challenge ofaccurately modelling the transient ignition process. Combustion in a homogeneous-charge natural gas fuelled SI engine with a relatively simple chamber geometryhas been simulated using an empirical model for ignition. An oxygen based reac-tion progress variable is employed as the conditioning variable and its stochasticbehaviour is approximated by a presumed probability density function (PDF). Atrajectory generated low-dimensional manifold has been used to tabulate chemistryin a hyper-dimensional space described by the reaction progress variable, temper-ature and pressure. The estimates of pressure trace and pollutant emission trendsobtained using CSE accurately match experimental measurements.iiPrefaceThis thesis has been the result of collaboration with researchers from the Univer-sity of British Columbia and Westport Innovations Inc. in Vancouver, Canada.The work has been funded by Mitacs Inc., Canada and Westport Innovations Inc.,Canada. I am responsible for conducting all parts of the research and preparing thecorresponding manuscripts.? Chapter Three: contains work that was accomplished at the University ofBritish Columbia under the supervision of Dr. Kendal Bushe. Part of thework has been published as conference proceedings under the title:? Girish V. Nivarti, M. Mahdi Salehi, W. Kendal Bushe, Conditionalsource-term estimation for large eddy simulation of turbulent premixedflames in complex geometries, Proceedings of the Combustion InstituteCanadian Section, 03-001, 2013.Ms. Nasim Shahbazian and Dr. Jim Huang contributed geometry data forconducting validation tests of the developed method. A licensed copy ofTecplot360 software was used for plotting domain geometries. I am respon-sible for conducting all parts of the research and preparing the correspondingmanuscripts. This work shall soon be submitted for journal publication.? Chapter Four: involves research work conducted at Westport InnovationsInc. under the supervision of Dr. Jim Huang and Dr. W. Kendal Bushe undera non-confidentiality agreement signed through Mitacs Accelerate and Mi-tacs Globalink. The open-source software package, OpenFOAM was usedfor performing relevant simulations. A licensed MATLAB software copywas used for performing low-pass filter calculations. I am responsible foriiiPrefacewriting the entire C++ code pertinent to the combustion model. I am respon-sible for conducting all parts of the research and preparing the correspondingmanuscripts. This work shall soon be submitted for journal publication.The entire manuscript was prepared using the LATEX package.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.1 Reynolds and Favre Averaging . . . . . . . . . . . . . . . 82.1.2 Filtering and Spatial Averaging . . . . . . . . . . . . . . 102.1.3 Probability Density Function Methods . . . . . . . . . . . 11vTable of Contents2.2 Turbulent Combustion . . . . . . . . . . . . . . . . . . . . . . . 122.2.1 Combustion Chemistry . . . . . . . . . . . . . . . . . . . 132.2.2 Regimes of Combustion . . . . . . . . . . . . . . . . . . 152.2.3 Chemical Closure Methods . . . . . . . . . . . . . . . . . 183 Spatial Partitioning Algorithm . . . . . . . . . . . . . . . . . . . . . 253.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.1.1 Relevant Background . . . . . . . . . . . . . . . . . . . . 263.1.2 Space-filling Curves . . . . . . . . . . . . . . . . . . . . 283.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.3 Test Geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494 Engine Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.1.1 Relevant Background . . . . . . . . . . . . . . . . . . . . 514.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.2.1 Conditional Source-term Estimation . . . . . . . . . . . . 544.2.2 Presumed Probability Density Function . . . . . . . . . . 564.2.3 Chemistry Reduction . . . . . . . . . . . . . . . . . . . . 584.2.4 Governing Equations . . . . . . . . . . . . . . . . . . . . 604.2.5 Spark-Ignition Model . . . . . . . . . . . . . . . . . . . . 624.3 Problem Specifications . . . . . . . . . . . . . . . . . . . . . . . 644.3.1 Experimental Validation Case . . . . . . . . . . . . . . . 644.3.2 Numerical Solver . . . . . . . . . . . . . . . . . . . . . . 654.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.4.1 Flame Kernel Development . . . . . . . . . . . . . . . . 684.4.2 Pressure Trace . . . . . . . . . . . . . . . . . . . . . . . 714.4.3 Emissions . . . . . . . . . . . . . . . . . . . . . . . . . . 754.4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 774.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77viTable of Contents5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.1 Inverse Problem Implementation . . . . . . . . . . . . . . . . . . 795.2 Industrial Modelling Application . . . . . . . . . . . . . . . . . . 805.3 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . 81References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82viiList of Tables3.1 Specification and characteristic dimensions of test cases. . . . . 353.2 Mesh type and features for different test geometries. . . . . . . 363.3 Performance of spatial partitioning algorithm in different ge-ometries and inter-point distance parameters defining the corre-sponding Morton order curve. . . . . . . . . . . . . . . . . . . 484.1 Specifications of the Ricardo Hydra research engine. . . . . . . 65viiiList of Figures1.1 Diversity of time scales in a turbulent flame. . . . . . . . . . . 22.1 Kolmogorov?s energy cascade hypothesis. . . . . . . . . . . . 62.2 Classical turbulent premixed combustion regime diagram. . . 162.3 Turbulence-chemistry interaction in premixed flames. . . . . . 172.4 Conditional averaging of a scalar in a turbulent reactive flow. . 202.5 Chemical closure in conditional source-term estimation. . . . 233.1 Contiguity of linear mapping within partitions using the (a)Hilbert curve, and the (b) Morton order curve. Leaps betweenconsecutive points can be significantly larger in the Mortonorder curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.2 Equi-sized partitioning of complex or refined mesh structuresresults in poor spatial locality within clusters. Regions of sim-ilar shading are assigned the same cluster. . . . . . . . . . . . 303.3 Partitioning of the Morton order curve between points locatedacross leaps increases cluster locality at the expense of losingcluster size balance. . . . . . . . . . . . . . . . . . . . . . . . 333.4 Distances between consecutive points (block-centres) on theMorton order curve displayed for the first two steps of the par-titioning algorithm. The normalized distance between points(r?i ? ri/rmax) has been employed for the case with nb = 3200,and k = 128. . . . . . . . . . . . . . . . . . . . . . . . . . . 37ixList of Figures3.5 Distances between consecutive points (block-centres) on theMorton order curve displayed for final steps of the partition-ing algorithm. The normalized distance between points (r?i ?ri/rmax) has been employed for the case with nb = 3200, andk = 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.6 Degree of computational load balance for the Bunsen burnergeometry gauged using relative cluster (c?i) for strictly equalpartitioning. Inter-processor communication is indicated bythe locality metric (LR,c). . . . . . . . . . . . . . . . . . . . . 403.7 Degree of computational load balance for the Bunsen burnergeometry gauged using relative cluster (c?i) for leap-based par-titioning. Inter-processor communication is indicated by thelocality metric (LR,c). . . . . . . . . . . . . . . . . . . . . . . 413.8 Comparison of intra-cluster (LR) and inter-cluster (LR,c) local-ity measures for manual partitioning (MP), equi-sized parti-tioning (EP) and using the spatial partitioning algorithm (SPA)on the Morton order curve in the Bunsen burner geometry fork = 128 divisions. . . . . . . . . . . . . . . . . . . . . . . . 423.9 Clusters produced from the equi-sized partitioning of Mor-ton order curve often contain leaps that reduce spatial localitywithin. Even for a simple cylindrical Bunsen burner geometry,leaps occur at frequent and irregular intervals. . . . . . . . . . 433.10 A qualitative picture of the difference between regular equi-sized partitioning and our spatial partitioning algorithm ap-plied to Bunsen burner geometry shown above. . . . . . . . . 443.11 Computational load balance measures (c?i and Lrms,c) estimatedfor the bluff-body burner for different cluster divisions, k =32,64 and 128. . . . . . . . . . . . . . . . . . . . . . . . . . 453.12 Computational load balance measures (c?i and Lrms,c) estimatedfor the Ricardo Hydra engine for cluster divisions, k = 32,64and 128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.13 Sections of complex geometries clustered using the spatial par-titioning algorithm along with appropriate sectional views. . . 47xList of Figures4.1 Laminar flame solutions of mixtures with varying relative air-fuel ratio (? = 1,1.1, . . . ,1.5) used for the regularisation ofCSE solutions. . . . . . . . . . . . . . . . . . . . . . . . . . 564.2 The ? -probability density function (PDF) has been tabulated apriori for different means and variances. A sample distributionof the two presumed PDFs employed has also been illustrated. 574.3 Conditional reaction rates at 25 bar obtained directly comparedwith an interpolation between tabulated values at 10 bar and 30bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.4 Iterative algorithm of the CSE-TGLDM chemical closure model. 624.5 Bowl-in piston geometry of the Ricardo Hydra engine. . . . . 644.6 Axi-symmetric mesh constructed at 270? crank angle. . . . . 664.7 Development of flame kernel shown as the degree of MFB until5% MFB and corresponding conditionally averaged tempera-tures obtained for a lean mixture with ? = 1.5 at 2500 rpm. . 674.8 Mean and variance contours of progress variable at 5% MFBfor a mixture with ? = 1.5 at 2500 rpm. . . . . . . . . . . . . 684.9 Initiation timing of high-reactivity flame kernel with ? 2.5%chamber volume for varying ? and their combustion durations. 704.10 Pressure map for a lean mixture (? = 1.5) at 2500 rpm. . . . . 714.11 Heat release rate for a lean mixture (? = 1.5) at 2500 rpm. . . 724.12 Variation of MFB for a lean mixture (? = 1.5) at 2500 rpm. . 734.13 Peak pressure and associated crank angle for varying ? . . . . 744.14 Sensitivity of maximum cylinder pressure and the correspond-ing crank angle to the pressure at IVC in the simulations. . . . 764.15 Pollutant emission trends with varying ? . Numerical estimateshave been scaled appropriately. . . . . . . . . . . . . . . . . . 77xiList of Algorithms3.1 Spatial partitioning algorithm for multi-block meshes. . . . . 313.2 Threshold consecutive point distance for partitioning. . . . . . 323.3 Cluster coarsening for reviving computational load balance. . 34xiiList of Symbols? Scalar conditional average? Laminar flame thickness?i j Kronecker delta function?? Chemical production rate?0 Integral length scale? Dissipation rate of energy? Kolmogorov length scale? Relative air-fuel ratioD Molecular diffusion constantJ Molecular diffusion flux? Coefficient of laminar kine-matic viscosity?T Coefficient of turbulent kine-matic viscosity? Fluid density~x Spatial coordinatec? Favre-averaged value of cc???2 Favre-averaged variance of cc Reaction progress variableDa Damk?hler numberk Kinetic energyKa Karlovitz numberP PressureRe Reynolds numberRet Turbulent Reynolds numberSL Laminar flame speedSc Schmidt numberT Temperatureu? Turbulence intensityui ith component of velocityYk kth species mass fractionxiiiList of AcronymsBML Bray-Moss-LibbyCMC Conditional Moment ClosureCSE Conditional Source-term EstimationDNS Direct Numerical SimulationEBU Eddy Break-upFSD Flame Surface DensityLES Large-eddy SimulationLEM Linear-eddy ModelMMC Multiple Mapping ConditioningODT One-dimensional TurbulencePCM Presumed Conditional MomentPDF Probability Density FunctionRANS Reynolds Averaged Navier StokesTGLDM Trajectory Generated Low-dimensional ManifoldAMR Adaptive Mesh RefinementNP Non-deterministic Polynomial-timexivList of AcronymsRSB Recursive Spectrum BisectionRCB Recursive Coordinate BisectionSFC Space-filling CurveATDC After Top Dead CenterBTDC Before Top Dead CenterCA Crank AngleEVO Exhaust Valve OpensGIW Gross Indicated WorkGRI Gas Research InstituteHRR Heat Release RateIVC Intake Valve ClosesMFB Mass Fraction BurnedSI Spark-IgnitionTDC Top Dead CenterxvAcknowledgementsTurbulent combustion is a formidable subject of inquiry. Over the past few years ofresearch, I have evolved as a person through learning from my mistakes. Through-out my errors, doubts, and fears, my supervisor, Prof. W. Kendal Bushe, has shownme courage and calm. In his sheltering aura of knowledge and sheer cool, I haverisen to a stage where I can comprehend the basic questions within our discipline.I owe this work and the related successes in my life entirely to him.My academic growth has been accelerated by the impeccable mentorship ofmy professors and collaborating researchers. Foremost amongst them, Dr. MahdiSalehi has been a friend, philosopher and guide to me. I humbly thank him for hiscompany and his benefactions, academic and otherwise. I sincerely thank my pro-fessors, Prof. Neil Balmforth, Prof. Mona Berciu and Prof. Carl Ollivier-Gooch,whose passion towards their disciplines has been inspirational. I am grateful to Dr.Jim Huang for introducing me to the world of industrial problems which was oncebizarre to me. I also thank my colleagues Hong Tsui and Graham Hendra for theirpositive company.I owe my relentlessly optimistic attitude to my mother, my father and my sisterwho, through their love and support, allowed me to grow healthily and live pas-sionately in this continent half the world away from home. I am deeply indebtedto Lorenzo Lane, Dr. Artan Sheshmani, Amanda Grochowich, Razan Hamdi andTim H?llering for their true friendship, inspiration, love, guidance and science! Idedicate this thesis to the land where I first delved into research, where some of themost memorable moments of my life were experienced and where, today, I standbeckoned by the vast ocean of learning. . .xviTo Canada.xvii1 | IntroductionEnergy has remained vital to the growth of our civilisation for centuries. Butover the past few decades, the rising population and rapid urban development havecaused an acute inflation of energy demand, encouraging a globally pervasive de-bate on energy policy. A consensus has emerged mandating the alleviation of do-mestic and industrial carbon footprints on the environment. Amidst the ongoingpolemic on current power generation technologies, the Stern review [1] categori-cally decries fossil fuel combustion as a source of harmful pollutants and green-house gases. Today, low-polluting renewable sources of energy form a topic ofsignificant research. Moreover, a growing conscience for our future generationshas initiated a shift towards sustainable lifestyles.The ground reality, however, is that building an infrastructure of alternativeenergy hinges on existing resources. Burgeoning demands of our society cannotbe fulfilled by renewable energy, particularly with the current pace of develop-ment. The significance of combustion has been clearly illustrated by Tollefson andMonastersky [2] among others. While fossil fuel reserves are said to be declining,new repositories of combustible fuels are being discovered and invented in formshitherto unknown. All roads to the immediate future of our civilisation dependinevitably on combustion ? the challenge is to use every drop of fuel judiciously,with diminishing environmental effects.Scientific interest in combustion research has been fuelled further by the longadvancing strides of technology with each finding. In contrast to their counterpartsfrom early twentieth century, combustion systems today are outcomes of prolifictheoretical and experimental research. Contemporary aircraft engines are muchmore efficient, cleaner and safer than their immediate precursors. Although ex-11. Introductionperiments have been crucial to these drastic developments in engine design, theyhave limited capability to probe the intricate geometries of modern devices. Theadvent of high performance computing technology has not only aided experimentalresearch in this regard, but it has become a necessary predictive tool in the analysisof alternative fuels. In fact, Moin and Kim [3] describe supercomputers as pow-erful assets which might uncover the secrets long held behind natural phenomenathat characterise the behaviour of flames.101 10?1 10?3 10?5 10?7 10?910?2 10?4 10?6Unit: sTurbulenceCombustion ChemistryFigure 1.1: Diversity of time scales in a turbulent flame.In combustion devices of engineering interest, fuel mixtures burn as turbulentreactive flows. Chaotic properties of turbulence impart a disparity of scales to theflow physics, while the chemistry in flames comprises an even broader spectrumof time and length scales. Nitrogen oxides (collectively termed, NOx) are prod-ucts of gradual transformation at high activation temperatures, whereas myriadsof radicals such as hydroxyl ions (OH?) survive fleeting nanoseconds. As dis-cussed by Pope [4], small scales and many species are but among the manifoldother formidable challenges that render turbulent combustion to be mathematicallyand computationally prohibitive for decades to come. Even the existence of an an-alytical solution to the non-linear governing equations of turbulent flows is an openMillenium Prize Problem, formally stated by Fefferman [5].21. IntroductionIn the absence of an exact solution, computational fluid dynamics (CFD) hasemerged as an approach to numerically estimate solutions of time sensitive industryscale problems. Fundamental aspects of flame behaviour, including ignition, ex-tinction and flame propagation, are modelled based on experimental knowledge,asymptotic analysis, theoretical insight and often, scientific faith. Combustionmodelling seeks to resolve the interaction between chemistry and turbulence atsmall scales; this is infamously known as the chemical closure problem. Whenthe fuel and oxidiser are mixed prior to combustion, the chemical reactions andturbulent flow exhibit a stronger coupling, particularly at the high turbulence in-tensity observed in engines. Premixed flames, however, can propagate in lean fuelmixtures and demand relatively low ignition temperatures. These conditions pro-vide an unsuitable environment for the emission of pollutants such as unburnedhydrocarbons, carbon monoxide (CO) and NOx. As modern engines increasinglyexploit a greater fraction of the premixed burn of fuel, there is an urgent need fordeveloping precise modelling tools.Conditional source-term estimation (CSE) is a combustion model formulatedby Bushe and Steiner [6] in 1999 at the Center for Turbulence Research, Stan-ford. In contrast to previously developed models for chemical closure, CSE is notrestricted to any specific regime of flame behaviour. CSE is poised to accuratelydescribe turbulence-chemistry interactions at high turbulence intensities due to thelack of common binding assumptions. During the past decade, this model has beensuccessful in calculating quantities of practical interest for canonical non-premixedflame problems; the application to premixed flame simulations has been relativelyrecent. The subject of this thesis is to extend the applicability of CSE to problemsof industrial context. The work has been presented as a part on the theoreticaldevelopment of CSE and another on its application to problems in contemporarynatural gas engines.32 | BackgroundAccurate modelling of turbulent reactive flows hinges on a sound understandingof the complex constitutive processes. Knowing the fluid mechanical propertiesand their variation with the temperatures exhibited by a combustion system is ab-solutely essential for describing transport phenomena and other mechanisms suchas heat transfer, molecular diffusion and convective turbulent transport [7]. Uponmixing, reacting species undergo chemical transformations that follow non-linear,inter-dependent mechanisms. Pollutant formation and fuel consumption are gov-erned by these mechanisms; in particular, nitrogen undergoes transformations toproduce environmentally dangerous oxides, NOx [8]. Even methane, the simplestof simplest of hydrocarbons, may undergo over 340 elementary reactions involv-ing up to 50 intermediate species upon combustion in the presence of air [9]. De-tailed mechanisms are often simplified and tabulated prior to simulation in orderto improve computational efficiency [10]. Tabulation of reaction constants andchemical reaction mechanisms has not been accomplished for higher hydrocar-bons [11, 12]. In gas phase combustion, the energy released from reactions affectsthe transport properties of the system and enhances turbulent micro-mixing, whilstthe increased micro-mixing helps escalate reaction rates by bringing reactants to-gether faster. Resolution of these non-linear turbulence-chemistry interactions is apersistent challenge in gas-phase combustion. This chapter provides a conceptualbackground of mathematical investigative tools necessary for developing modelsto resolve this coupling.42. Background2.1 TurbulenceTurbulent flows present a non-linear, multi-scale problem described in completedetail by the inherently coupled system of Navier-Stokes equations. Originallywritten as integral conservation laws by White [13], these differential equationsdefine conservation of mass and momentum assuming the continuum hypothesis,wherein flow variables may be represented as smooth functions in space and time[14]. Conservation of mass principle (continuity equation) is written as??? t +??u j?x j= 0 (2.1)where ? is the density of the fluid and u j is the velocity in the jth principal direction.Conservation of momentum is described by the following transport equation??ui? t +??uiu j?x j=? ? p?xi+ ??i j?x j+Fi. (2.2)where ?i j represents the stress tensor and p, the pressure at each point in the fluid,while Fi accounts for external forces. For the constant-property Newtonian fluids,this can be written as?i j =?p?i j +?(?ui?x j+ ?u j?xi)(2.3)wherein ?i j is the Kronecker delta and ? is the viscosity of considered fluid whereas?p is the normal static pressure observed at any point. In addition with an equa-tion for energy used in one of its various forms, the system is closed using anappropriate equation of state [15].Although this system of equations is readily solved for laminar flows, no knownanalytical solutions exist for turbulent flows even in physical simple domains suchas pipes. Features that characterise turbulent flows have been known for centuries,such as high rates of diffusivity, 3-dimensional irregularities and vorticity, but thetransition to turbulence has not been well described [16]. In the absence of suchcritical knowledge, turbulence research relies on physical observations and numer-ical analysis for valuable insight into this ubiquitous phenomenon. As established52. Background?? ?0 ?L?0? ?? ??? ?Universal equilibrium range Energy containing rangeFigure 2.1: Kolmogorov?s energy cascade hypothesis.by Reynolds [17], turbulent flows are characterized by a single non-dimensionalparameter based on the characteristic velocity (U ) and length scales (L ) of theflow, defined asReynolds number, Re?U ?L /? (2.4)where ? is the kinematic viscosity of the fluid. At high Reynolds numbers, typicalfor most turbulent flows, a disparity of scales is observed [14]. The energy cascadehypothesis illustrated in Fig. 2.1 describes the transfer of kinetic energy (k) fromthe largest scales (?0), which are determined by the scale of the domain geometry(L ), to the smallest scales where all the energy dissipates by viscous action. Kol-mogorov hypothesised that in a turbulent flow with high enough Reynolds number,the statistics of motion at the smallest scales are isotropic and are universal to allturbulent flows. The smallest scales can be determined solely by ? and the dis-sipation rate, ? , and thus all information of geometry is lost at these scales. Theunique scaling of various parameters at play results in a universal form for the ratioof largest scales, ?0, and the smallest, Kolmogorov scales defined by? ? (?3/?)1/4 (2.5a)for the length scales, and by?? ? (?/?)1/2 (2.5b)62. Backgroundfor the time scales [16, 18]. The relative magnitude of the smallest length and timescales, ? and ?? respectively, referred to as the Kolmogorov scales [14], comparedto the largest scales, l0 and ?0, is also approximated in a simple form by usingKolmogorov?s hypotheses as?/?0 ? Re?3/4 (2.6a)??/?0 ? Re?1/2 (2.6b)where ?0 and ?0 represent integral time and length scales respectively. An ex-act numerical solution of Navier-Stokes equations is necessary to obtain the flowvariables, such as velocity and pressure, as known functions of space and time.The computational cost of a numerical simulation are determined by the resolutionrequirements ? the solution domain must be large enough to contain the energy-containing motions, and the grid spacing ?x must be small enough to resolve thedissipative scales [14]. However, it is evident from Eq. 2.6 that the instantaneousrange of scales in a turbulent flow increases rapidly with the Reynolds number.Majority of engineering problems have too wide a range of scales to be directlycomputed ? given the necessary grid spacing, the number of grid points requiredvary as Ng ? Re9/4 [14], which when combined with the time resolution require-ments lends the alarming result given byComputational effort ? Re11/4. (2.7)It is evident that resolving all scales of motion using Direct Numerical Simu-lation (DNS) for practical flows is infeasible [19]. The vast information obtainedusing DNS of a practical combustion device often has little value since quantitiessuch as mean pollutant formation rates, average power, and average fuel consump-tion are of greater interest to the industry. Significant progress has been madetowards a faster estimate of these parameters by describing the flow as a chaoticprocess using random variables. Simpler descriptions seek to reduce computationalexpense: instead of solving for the instantaneous flow-field, a statistical time evo-lution of the flow is sought. Not only are the computational abilities inadequate,knowledge of initial conditions is seldom precise, accurate or sufficient. Be it ex-72. Backgroundperimental measurements or numerical estimates of initial and boundary values, thedifference between realisation and intention is often disparate. The mathematicaltreatment of turbulent flows using random variables is a useful exercise that helpsachieve reasonably accurate predictions. However, these methods suffer from thedrawback of introducing terms which cannot be obtained from the governing equa-tions and thus require modelling.2.1.1 Reynolds and Favre AveragingThe starting point of statistical analysis of turbulent flows is the decomposition ofeach quantity into a mean and a fluctuating component as follows:Q= Q+Q?, (2.8a)where the components follow the basic moment rules:Q= Q, Q? = 0. (2.8b)Reynolds Averaged Navier Stokes (RANS) equations describe the transport of av-eraged quantities through the Navier-Stokes equations using this decomposition.The immense computational expense involved in DNS is, therefore, avoided bycalculating only the mean flow field. In compressible flows such as turbulent reac-tive flows, the terms are instead Favre-averaged by considering a density weightedaverage to simplify the resulting equations while still including the effect of densityfluctuations. Decomposition into Favre-terms is written asQ= Q?+Q?? (2.9a)where the Favre-average is defined by the principlesQ?? ?Q? , and Q??? =?(Q? Q?)? = 0. (2.9b)No simple relation exists between Favre and Reynolds averages ? a mean-ingful relation depends on density fluctuation correlations ? ?Q? which are im-82. Backgroundplicit in the Favre-average. More importantly, since experimental techniques de-termine Reynolds averaged data, comparing experimental results with Favre av-eraged quantities obtained using numerical simulations might not be obvious [7].The main reason for Favre-averaging is compactness of resulting equations. TheFavre-averaged momentum equations are derived as?(? u?i? t + u? j? u?i?x j)=? ? p?xi+ ??x j(??i j?? u???i u??j ). (2.10)wherein u???i u??j is not closed within the governing equation system. The RANS orFANS decompositions thus result in the creation of non-linear terms which areunresolved in the equation system. In general, this term describing the averageexchange of turbulent momentum across the boundary of each finite volume isdescribed by a simple and inaccurate relationu???i u??j ?23k?i j = ?T(?ui?x j+ ?u j?xi)(2.11)Turbulence modelling then reduces to the art of estimating turbulent viscosity, ?T .An overwhelming majority of such models calculate this quantity using empiricaland scaling laws. Typically, turbulent energy spectrum parameters such as the dis-sipation rate ? and turbulent kinetic energy k are used. The most popular techniqueis to use the scaling?t =C??k?2?? (2.12)to obtain a rough estimate. The value of C? depends on the flow and hence, thismodelling technique is not complete in the sense that it requires to be tuned to suita particular problem. A similar approach solves a separate transport equation forthe turbulent frequency ? = k/? and is substantially more accurate in near-walllayers [20]. As a unique strategy, [21] developed a transport equation for turbulentviscosity which is a popular choice in external flow applications such as aircraftsimulations. In the recent past, several models have aimed to calculate Reynoldsstress using transport equations. Although this method can resolve the stress tensorwith greater accuracy it is not robust for industrial application.In general, the RANS concept must be used with discretion and care. The92. Backgroundapproach is constrained by time and spatial resolution requirements dictated bythe numerical solution procedure. A RANS simulation is, however, a fast meansto obtain basic insight into the flow field. The low computational expense andsimplicity makes RANS the most popular approach in industry.2.1.2 Filtering and Spatial AveragingThe objective of Large Eddy Simulation (LES) is to calculate the largest struc-tures of the flow exactly and model the small-scale behaviour. LES resolves thevelocity field obtained by applying a low-pass filter of characteristic width ? tothe underlying turbulent velocity field, ui(~x, t). Large-scale behaviour of turbulentflows is known to be governed by the domain geometry and motion at small-scalesis believed to have universal properties. In turbulent flows and turbulent reactiveflows, critical features including instabilities and unsteady mixing are governed bythe large structures which can be computed exactly by LES. The general filteringoperation is defined by?ui(~x, t)?=?VG(~r,~x)ui(~x?~r, t)d~r (2.13a)where the specified filter function G satisfies the normalization condition?G(~r,~x)d~r = 1 (2.13b)and integration is over the entire flow domain [14]. Standard filters include i) acut-off filter in spectral space which preserves length scales greater than 2?, ii)box filter in physical space which corresponds to a spatial averaging of Q over abox with size ?, and iii) a Gaussian filter. Unlike Reynolds averaging, the filteringoperations do not follow the moment relations, i.e.:??Q 6= Q?, Q?? 6= 0. (2.14)However, unresolved stresses defined as follows in LES have a similar expres-sion as in RANS and require sub-grid scale modelling:?Ri j ? u?iu j? u?iu? j. (2.15)102. BackgroundFiltered balance equations coupled with appropriate sub-grid models can beused to numerically solve the behaviour of filtered fields. Although LES lacks theaccuracy level of DNS, a drastically more detailed picture of the turbulent flowis obtained in contrast RANS. LES has advanced further as a powerful tool forthe simulation of turbulent flows in the past decade [22]. However, application toturbulent reactive flows of even canonical industrial problems are challenging par-ticularly due to the complex geometry that demands finer mesh resolution therebyincreasing computational expense [23, 24]. Pope [25] highlights the manifold chal-lenges including completeness and tractability of LES as a precise mathematicalmodelling tool; most industrial problems remain beyond the scope LES .2.1.3 Probability Density Function MethodsAn averaged quantity, computed using RANS, represents the first moment of thefield. However, given the one-point, one-time joint Cumulative Density Function(CDF) of velocity in a turbulent flow-fieldF(~u;~x, t)? P{ui(~x, t)<Vi, i = 1,2,3}, (2.16a)the joint probability density function (PDF) can be computed asf (~u;~x, t) = ?3F(~u,~x, t)?ui?u j?uk. (2.16b)At each point in in space, the PDF characterizes the random velocity vectorwithout holding any joint information for two or more points in space or time[14]. PDF methods solve a transport equation for these PDFs using stochastic La-grangian methods. At high Reynolds numbers, molecular diffusion has a negligiblecontribution to spatial transport and convection dominates the transport of momen-tum, chemical species and enthalpy. Despite the accurate treatment of convectiveterms in the Lagrangian PDF method, however, the computational cost involved inMonte Carlo simulations render these methods infeasible for industrial application.At one end of the spectrum of approaches for the simulation of turbulent flows,DNS provides complete detail at an extremely high computational cost and, hence,a limited range of applicability [19]. RANS-based models such as the k?? method112. Backgroundlie at the other end; these offer a diverse range of applicability but with poor accu-racy and limited insight [26]. LES and PDF methods form the middle range andare gaining popularity with the rapidly developing high performance computingtechnology.2.2 Turbulent CombustionTurbulent combustion has been a topic of intense scientific inquiry for severaldecades [27?30]. Several questions remain open due to the formidable complexityof problems involved in a turbulent reactive flow particularly due to heat releaseand differential diffusion [4, 31].There are two principal modes of combustion based on whether fuel and oxi-diser are mixed prior to combustion (premixed combustion) or not (non-premixedcombustion). In reality, most reactive flows exhibit a hybrid behaviour and arecategorised under partially-premixed combustion. Premixing provides for effec-tive control of stoichiometry of flames ? this allows lean-burning and reduction ofunburned hydrocarbons while providing an unsuitable environment for NOx forma-tion mechanisms as understood today. Due to these reasons premixed combustionis being employed in gas turbines, and is being theoretically investigated in or-der to develop efficient automobile engines and industrial furnaces. On the otherhand, premixed fuel storage is a known safety hazard. Moreover, premixed flamesare susceptible to convectively and acoustically coupled instabilities. From a the-oretical standpoint, the motion of turbulent premixed flames is a superposition offlame propagation and turbulent fluid convection. This renders their mathemat-ical modelling significantly more challenging than that of non-premixed flames.Greater understanding of premixed flames and development of accurate predictivetechniques is, therefore, absolutely necessary. Since the non-reactive aspects in aturbulent reactive flow have been already discussed, the following sub-sections fo-cus on the two main processes that chemical species undergo ? transformation andtransport.122. Background2.2.1 Combustion ChemistryIn hydrocarbon combustion, 50?7000 reacting species may be actually involved ina complete mechanism describing the transformation of the fuel-air mixture intoproducts. However, known chemical mechanisms for most fuels used with CFDanalysis comprise 150?250 [4]. Solving transport equations for even the decreasednumber of species is a computationally daunting task. Therefore, reduced mech-anisms that represent the chemistry in terms of few key species and scalars aresought as a chemistry model for most applications of practical concern.Species TransportThe generalised conservation equation of chemical species is written as?? t ?Yk+??x j(?u jYk) =???x jJ j,k + ??k, (2.17)where J j,k is the molecular diffusion term and ??k is the production rate of the kthspecies. Species molecular diffusion is generally described using Fick?s law:J j,k =??DkYkx j(2.18)where Dk is the coefficient of molecular diffusion for the kth chemical species.Premixed flames are adequately and inexpensively defined by a reaction progressvariable (c) in low-Mach number flows under adiabatic conditions and unity Lewisnumber [32]. Typically, this may be defined using a product species asc= YPYP,eq, (2.19a)or alternatively, using reactant concentrations asc= YR?YR,iYR,eq?YR,i. (2.19b)where YR,i and YR/P,eq are the initial and equilibrium concentrations of appropri-ately chosen reactant, R or product, P. In general, the reaction progress variable132. Backgroundmust transform monotonically across a flame, and preferably its gradient is evenlydistributed across the reaction. For such a progress variable, the Favre-averagedspecies conservation equation can be written as?? c?? t +?? u? j c??x j= ??x j(?Dc?c?x j)???u??j c???x j+ ??c (2.20)where ??c describes the averaged transformation rate of the progress variable. Thedensity weighted averaged chemical diffusivity remains roughly constant; thus, thefollowing simplification is achieved?Dc?c?x j? ?Dc? c??x j. (2.21)Chemical transport equations in turbulent combustion using RANS decompo-sition, therefore, results in yet another closure problem ? resolving averaged tur-bulent chemical transport and averaged chemical source-terms. Species turbulenttransport is typically modelled using a gradient diffusion assumption?u??jc?? =??DT? c??x j(2.22)where DT is an effective turbulent molecular diffusivity of the species associatedwith the progress variable, c. Although this model for species turbulent transport iswidely debated due to the known existence of counter-gradient diffusion regimes,the aforementioned model provides an inexpensive closure. It is, thus, widely used,particularly in the context of RANS.Chemical TransformationThe forward reaction rate of a generic chemical transformation of S species via Mreactions represented concisely asS?s=1?rm,s [A]sk f??S?s=1? pm,s [A]s with m= 1, . . . ,M, (2.23)142. Backgroundis a non-linear function of the concentration of the participating reacting species([A]s). The production rate of the kth reacting species is expressed as??k =M?m=1km(? pm,k??rm,k) S?s=1[Ys]?rm,s (2.24)wherein km is the reaction rate coefficient and Ys denotes the concentration of thesth species. To aggravate the non-linearity of the calculations, rate coefficientsof chemical reactions bear a strong non-linear dependence on temperature, whichgoverns activation energy of the reaction. The relation can be expressed in theArrhenius form [32] askm = A ?T? ? exp(? EaRT)(2.25)where the constant A is termed the pre-exponential factor and ? is the non-linearityconstant of temperature. The activation energy Ea, quantifies an energy barrier thatmust be overcome for the reaction to proceed in the forward direction.It is evident that averaged chemical reaction rates are non-linear functions ofscalar fields. These cannot be closed using other quantities solved for in the tur-bulent combustion system, i.e. the averaged source-term for kth species or thecorresponding progress variable cannot be described by the averaged flow scalarsfrom the chemistry model.??c(? ,Yk,T ) 6= ??c(?,Yk,T ), (2.26)2.2.2 Regimes of CombustionGiven the complex, non-linear law describing the burning rate ??c, a physical ap-proach is often sought to derive models for turbulent combustion. Physical analy-sis of chemical source-terms relies on a knowledge of the chemical and turbulentscales involved in a turbulent reactive flow. Representation of flame behaviour andits descriptive turbulence-chemistry interactions is the subject of the Borghi dia-gram [33, 34]. The objective is to analyse premixed turbulent combustion usingcharacteristic length and time scales corresponding to chemistry and turbulence.152. Background100100101101102102103103104104Integral length scale / flame thickness, ?0/?RMSvelocity/flamespeed,u?/SLReT = 1? c=? T?c=?ku? = SLWrinkled flameletsCorrugatedflameletsDistributedreactionzonesWell stirred reactorLaminarcombustionAircraft enginesPiston enginesFigure 2.2: Classical turbulent premixed combustion regime diagram.Fig. 2.2 shows the diagram wherein the oncoming turbulence intensity normalizedby the laminar burning velocity u?/SL is plotted against the largest turbulent eddylength scale normalized by the laminar flame thickness ?0/? . The turbulent flow ischaracterised by the turbulent Reynolds numberReT ?u??0? (2.27)where u? is the turbulent velocity fluctuation which is related to turbulent kineticenergy as u? ??k, ?0 is the turbulence integral length scale and ? is the kinematicviscosity of the fluid. When the turbulent Reynolds number is smaller than one,laminar combustion is observed. The remaining portion of the diagram representsturbulent combustion which is well scoped by two other dimensionless ratios.The regime diagram in Fig. 2.2, illustrates the approximate placement of tur-bulent premixed flames observed in two common combustion systems ? aircraft162. BackgroundSL?? ??? ?0Products Reactants?? ?Figure 2.3: Turbulence-chemistry interaction in premixed flames.and automobile engines. The Damk?hler number compares macroscopic turbulenttime scales (?t ) and chemical reaction time scales (?c) asDa? ?t?c. (2.28)For large values of the Damk?hler number (Da?1), the flame front is thin andits inner structure is not affected by turbulence motions at small scales. Large scaleturbulent structures wrinkle the flame surface imparting strain in different regionsof the flame-sheet. At each point on this sheet, the flame may be described as aone-dimensional laminar flame element called flamelet. Exceedingly low Damk?h-ler numbers (Da? 1), on the other hand, correspond to slow chemical reactionswherein mixing always occurs prior to reaction as regime termed perfectly stirredreactor. A transition away from the flamelet regime occurs when the smallestturbulence scales have a time scale (?k) smaller than (?c). The Karlovitz numberdefines the Klimov-Williams criterion corresponding the value of the ratioKarlovitz number, Ka? ?c?k. (2.29)being equal to unity. This criterion helps delineate the flamelet regime from the dis-tributed reaction zones regime where the inner flame structure is strongly affectedby turbulent motions.172. BackgroundFig. 2.3 illustrates the structure of turbulent premixed flames in a high turbulence-intensity flow wherein Kolmogorov scale turbulent motion may occur within theflame. It is noteworthy that in any given turbulent flow system, a wide range ofdissipation occurs ? therefore, a turbulent premixed flame is most generally rep-resented by a zone that may comprise smaller turbulent structures within. Thetrue structure of turbulence-chemistry interactions forms a topic of significant de-bate [4]. Qualitative analyses of turbulent premixed combustion provide a startingpoint for several quantitative chemical reaction rate modelling techniques. The di-mensionless parameters that characterise the location of a premixed flame on thediagram can guide the choice of the most appropriate modelling strategy.2.2.3 Chemical Closure MethodsClosure of the chemical source-term can be addressed using one of several combus-tion models; other than basic algebraic correlations, chemical closure techniquescomprise three classes: flamelet models, CMC methods and PDF methods.Algebraic ModelsThe chemical source-term can be described as a function of easily acquired quan-tities such as turbulent mixing. However, such models neglect local stochasticbehaviour and must be tuned to match the turbulent premixed flame problem con-sidered. Although these descriptions are commonly used in the context of RANS[35], more accurate models have been sought for LES [7]. An important class ofmodels, termed the eddy break-up model (EBU), for example, assumes that chem-istry is fast. It is then, the mixing of reactants that is the rate determining processof average reaction rate. If the time scale of mixing is taken to be the of the orderof turbulence time scale, the following model is obtained?? =C ? ?k?c??2 (2.30)wherein the progress variable is used as a measure of scalar fluctuation magnitude[32]. The EBU and other similar models are widely used because of their simplic-ity. Although such models can be implemented in LES contexts, the incorporationof chemistry effects is difficult besides the fact that these are entirely ad hoc and182. Backgroundmay result in non-physical solutions [32].Flamelet ModelsThe flamelet assumption forms the basis for geometrical analyses of turbulent pre-mixed flames. The flamelet method, conceptualized by [36], resolves turbulence-chemistry interactions by describing the turbulent flame locally with characteristicsof a steady, one-dimensional laminar flame. In the context of LES, two widely usedtechniques comprise flamelet-based chemical closure [23]: the flame surface den-sity approach (FSD) and the G-equation model. In the (FSD) approach, the flameis identified as a surface which can be convected, diffused, curved and strained bythe turbulent velocity field [37]. The mean chemical source-term is calculated fromthe total reaction rate and the flame surface density??k = ??k? (2.31)where ??k is the integrated reaction rate of the surface of the flame and ?, the meanflame surface density. A similar approach introduces a flame wrinkling descriptionusing the ratio of flame surface and its projection in the direction of projection[38]. However, the accurate closure of the turbulent flux and propagation velocityof the flame surface, and the effects of curvature and strain rate on the flame surfaceremains a challenge [15, 23, 32].In contrast, the G-equation model borrows from the description of a laminarflame surface using a level set [23, 39]. This concept eliminates the need for re-solving the flame-front structure, but consistency with LES filtering has only beenachieved recently by [40]. Further, the lack of resolution of the flame-front struc-ture challenges the accuracy of chemical closure when the flame structure is ade-quately affected by turbulence. The thickened flame model was proposed by [41]as a method to resolve the flame structure on the mesh. However, the empiricalparameters used to thicken the flame are cannot be applied to arbitrary flame spec-ifications.As a less restrictive alternative to the geometric analysis of flame-fronts, thestatistical behaviour of a turbulent reactive flow is represented using a PDF. Inpremixed laminar flames, the statistics of scalars are well characterized by the re-192. Background00.510 0.5 1Normalized Conditioning Variable, cNormalizedScalar,YkYk|cc?Yk|c?Yk|c??Y kY k? Illustrative experimental resultsFigure 2.4: Conditional averaging of a scalar in a turbulent reactive flow.action progress variable defined such that c = 0 in fresh gases and c = 1 in burntgas. One line of approach then is to presume the PDF of the progress variable andobtain chemical closure using conditional statistics (as shown in Fig. 2.4) basedon the solution of a laminar flame. The Bray-Moss-Libby model combines suchan approach, using a bimodal distribution, with a physical analysis which accountsfor differences in fresh and burnt gases. Variations of the flamelet/progress variablemethod have been applied to LES using a presumed PDF [42?44]. As mentionedearlier, and illustrated in Fig. 2.3, these flamelet models are only accurate in highDamk?hler regimes, but are insufficient for modelling high turbulence-intensityturbulent premixed flames [45].Conditional Moment ClosureFormulated by two researchers, Klimenko [46] and Bilger [47] independently, Con-ditional Moment Closure (CMC) is a well established chemical closure tool forthe simulation of non-premixed combustion [15, 32]. It has been used success-fully for industrial applications including spray combustion in diesel engines [48].CMC offers a promising closure tool for simulating premixed combustion since it202. Backgroundis devoid of the assumption of a prevailing combustion regime. Recently, Amzinet al. [49] simulated stoichiometric pilot-stabilized Bunsen flames using the RANS-CMC method. Further, an algorithm was recently proposed for the application ofCMC for LES of premixed flames [50].The CMC method is based on the concept that small fluctuations are experi-enced from the average within an ensemble of realizations which comply with agiven condition [47]. Fig. 2.4 illustrates that the deviation about an average acrossall possible conditions is seldom accurate. Whereas, for at each given condition,fluctuations about the average of the ensemble of realisations is minimal. Condi-tional averages of scalars are defined asQ(? ,~x, t)? ?Y (~x, t)|? (~x, t) = ? ? ? ?Y |? ? (2.32)where the angle brackets denote an ensemble average over an ensemble of realiza-tions of the flow and the vertical bar indicates that this average is conditional onthe condition, ? (~x, t) = ? that is, only those members of the whole ensemble thatmeet this condition are included in the average [47].Now, instead of presuming the reaction rate that is conditionally averaged onthe reaction progress variable, as in conventional flamelet based statistical ap-proaches, CMC calculates this term from the hypothesis??k|? ? ??k(T |? , Yk|? , ? |? ), (2.33)assuming small fluctuations of reaction rate about its conditional average. Trans-port equations are solved for conditionally averaged scalars required to estimatethe conditioned reaction rate. This adds to the dimensionality of the problem andmakes CMC computationally expensive. The conditionally averaged reaction rateobtained can be integrated by writing??k(~x, t) =? 10??k|c?(c?)P(~x, t;c?)dc? (2.34)where ??k|c?(c?) is estimated here through the CMC hypothesis and the presumedPDF, P(~x, t;c?) for any point in the domain is estimated using the first two mo-212. Backgroundments of the conditioning variable. As opposed to the flamelet method, CMC hasthe capability of accurately modelling the chemical source-term for flames beyondthe flamelet regime [33]. However, transport equations typically solved for con-ditionally averaged scalars consist of various unclosed terms which are not wellunderstood and are modelled inaccurately [31].Transported PDF modelsVarious other statistical methods have been developed such as the linear-eddymodel (LEM) and the one-dimensional turbulence (ODT) model by [51, 52] andthe broad class of transported PDF models developed for turbulent reactive flowsby Pope [53]. In this formulation, the chemical source-terms are closed exactly,whereas turbulent transport and molecular mixing must be modelled. The rela-tively recent multiple mapping conditioning (MMC) method combines features ofCMC, PDF methods and mapping closure models [54, 55]. Although these meth-ods provide excellent chemical closure, the computational expense of Monte Carlosolution discourages application to industry-scale problems.Conditional Source-term EstimationConditional Source-term Estimation (CSE) was formulated by Bushe and Steiner[6], Steiner and Bushe [56] as a chemical closure model for the LES of turbulentnon-premixed combustion. The general concept of CSE follows from the CMCmethod wherein the conditionally averaged source-term is approximated by condi-tionally averaged scalars relevant to the chemistry model from the CMC hypothesis(Eq. 2.33). Similar to CMC, the averaged chemical source-term is calculated by in-tegrating the corresponding conditionally averaged source-term with the presumedPDF at the given spatial coordinate with respect to the conditioning variable.However, while conventional CMC models solve transport equations for theconditionally averaged scalars, CSE employs an integral equation similar to Eq. 2.34.Each Favre-averaged scalar (? ) necessary for calculating the conditionally aver-aged reaction rate can be written in terms of the corresponding conditionally aver-aged scalar as??k(~x) =? 10?k|c?(~x,c?) P?(~x,c?)dc?, (2.35)222. BackgroundEnsemble, A~xTurbulent jet flameFigure 2.5: Chemical closure in conditional source-term estimation.for the corresponding spatial coordinate (~x) in the physical domain where theP?(~x,c?) represents the PDF at that point. This integral equation can be averagedover an ensemble of discrete points (A) in the domain (where the conditional scalaris constant) to yield an inverse problem for the calculation of the conditional scalar.Formally, clusters of points that preserve the statistical homogeneity of condition-ally averaged scalars are chosen as ensembles. The ensemble averaged conditionalscalar is then independent of the spatial coordinate and the integral equation iswritten as??k(~x j) =? 10?k|c?(c?) P?(~x j,c?)dc?, j ? A (2.36)for each point j within the ensemble. Conditionally averaged scalars have beenobserved to exhibit statistical homogeneity over spatially localized ensembles ofpoints [56]. Fig. 2.5 highlights a chosen ensemble of localized points in the com-putational domain of turbulent jet flame. With the conditionally averaged scalarnow invariant of the jth spatial coordinate (~x j) in these ensembles, the previousequation can be treated as the well-known Fredholm integral equation of the firstkind [57], written in the discrete form as~b= A~?, A ji =? c2c1P?(c?;~x j)dc?, (2.37)where P?(~x j,c?) is the kernel,~b is the known array of scalar averages in the entireensemble and ~? represents the conditionally averaged scalar. Integrations of the232. Backgroundpresumed PDF required for calculating A ji are done a priori and tabulated for mbins using lower (c1) and upper (c2) limits for the ith bin of conditioning variable.The Tikhonov method is implemented for regularising the solution of the integralequation using an L-curve approach for the choice of an optimal regularisationparameter as discussed by Salehi et al. [58].The application of CSE to the simulation of turbulent premixed flames wasaccomplished by Salehi [59]. As such, CSE has only ever been used for sim-ple canonical reference flames, both in the context of non-premixed and premixedflames. The following work focuses on advancing CSE for application to prob-lems of industrial relevance. Such problems generally involve transient phenom-ena which have not been modelled using CSE before. Moreover, industrial com-bustion chambers have complex geometries and bounding walls which make theaccurate implementation of CSE more challenging. In the next chapter, a robustmethod is developed to extend the applicability of localized CSE inversions to thecomplex geometries generally encountered in an industrial problem. The follow-ing chapter employs CSE for the simulation of a homogeneous-charge natural gasspark-ignition engine.243 | Spatial PartitioningAlgorithm3.1 IntroductionData classification is a frequently encountered solution approach in a wide rangeof scientific disciplines. Spatial partitioning is a classification technique that iscritical for obtaining accurate solutions of physical problems wherein constitutivelocal phenomena demand segregated treatment. In complicated domain geome-tries, solutions can be achieved conveniently by partitioning the problem spaceinto localized sub-domains. Moreover, with the advent of parallel processing tech-nology, mesh partitioning has become an essential tool to address large-scale prob-lems by dividing computational load across the available computing resources. Inthis context, size governs load balance whereas locality designates communicationbetween processors. The choice between load and communication is particularlydifficult to make for problems with multiple requirements.The partitioning technique presented in the current work is motivated by a strin-gent demand posed by a relatively recent model, termed Conditional Source-termEstimation (CSE), applied for the computational fluid dynamics (CFD) simula-tion of turbulent reacting flows [6, 58, 60, 61]. The key underlying assumption isthat conditionally averaged scalars requisite for the CFD solution are statisticallyhomogeneous in localized ensembles of points. In simple geometries, equi-sizedslices of the domain along its major axis have sufficed [56, 62]. However, auto-mated construction of ensembles is absolutely necessary for applying the modelPreprint to be submitted to Journal of Computational Physics.253. Spatial Partitioning Algorithmto the arbitrarily complex geometries that have industrial relevance. Also, spa-tial proximity must be maintained within each cluster to retain the validity of theunderlying assumption. In addition, parallel implementation of this large-scaleproblem demands that a measure of load balance be considered in the partitioningprocedure. Finally, frequent mesh changes enforced by industrial demands or byadaptive mesh refinement (AMR) impose a limit on the computational expense ofpartitioning.What makes this problem particularly challenging is that locality is not merelydesirable from a standpoint of computational efficiency, but is critical for calcu-lating a valid solution. Load balance has relatively less significance, and as such,were there to be a trade-off between load balance and strict locality, strict localitymust be preferred.3.1.1 Relevant BackgroundThe problem of clustering N objects into M groups of objects with similar prop-erties is a persistent challenge in several academic fields. Frequent applicationsare found in computer vision, pattern recognition, geographical information sys-tems, large-scale database management, cache performance and many other fields[63, 64]. Statistical classification, such as the cognitive mapping developed byJenks [65], seeks to reduce the variance within classes and maximize the vari-ance between classes. Hartigan [66] discusses the relevant statistical theory forclustering, whereas Dav? and Krishnapuram [67] present a unified view of robustclustering methods. These age old methods are efficient for applications with fewdimensions, but it is recognized that the curse of high-dimensionality pervades theproblem of clustering nearest point neighbours [68].In CFD problems, partitioning corresponds to clustering regions of mesh basedelements such as nodes, finite elements or finite control volumes (also called cells).For a multi-block finite volume mesh, popularly used in AMR implementations,blocks of cells can be divided into clusters assigned to different processors in thecomputing architecture. Such computational mesh partitioning has been the focusof extensive research [69?72]. The choice of adopted strategy depends on knowngeometric information and the trade-off between partition quality and computa-263. Spatial Partitioning Algorithmtional expense that best suits the problem. Dynamic problems are addressed usingfast methods, whereas relatively slower but highly customized procedures are for-mulated to match specific criteria such as those imposed by the turbulent reactiveflow model.Over the last few decades, a variety of mesh partitioning techniques have beendeveloped, majority of which fall into two broad classes. Geometric algorithmsincluding recursive bisection techniques and the SFC approach need only the lo-cal geometric information, such as block-centre coordinates. As an example, therecursive coordinate bisection (RCB) developed by Berger and Bokhari [73] iswidely accepted due to its simplicity and the wide range of its applicability. Theunbalanced recursive bisection strategy provides an advantageous modification ofthe RCB method [74] and Gilbert et al. [75] investigate yet another geometricmethod of dividing an irregular mesh into equal-sized pieces with few intercon-necting faces. In general, geometric partitioning is conceptually simple, so thecorresponding algorithms are fast and easy to implement. Additionally, the space-filling curve (SFC) approach, discussed at greater length in section 3.1.2, providesa unified and scalable data structure relevant for a wide range of CFD operations[76, 77].The other broad class of partitioning, graph partitioning, has been reviewed bySchloegel et al. [78]. In this method, nodes of a computational graph representtasks that can be executed concurrently and edges represent the communicationrequired between tasks from one iteration to the next. Domain partitioning usingsuch graphs is challenging and is recognized as an NP-hard problem [79]. Thegreedy algorithm developed by Farhat and Lesoinne [80], builds each partition bystarting with a vertex and adding adjacent vertices until the target size or expectedcomputational load has been reached. A relatively advanced technique, termedrecursive spectral bisection (RSB), was developed by Simon [69]. It provides ex-cellent quality partitions at exceedingly high computational costs. Many standardgraph-based partitioning packages have been developed to employ efficient globaland local graph-based methods for domain partitioning. Teresco et al. [81] providean extensive review of partitioning algorithms.As mentioned earlier, the choice of adopted strategy depends on known in-formation and the expense-quality trade-off that best suits the problem. The SFC273. Spatial Partitioning Algorithm(a) Hilbert-Peano curve(b) Morton order curveFigure 3.1: Contiguity of linear mapping within partitions using the (a)Hilbert curve, and the (b) Morton order curve. Leaps between consecutivepoints can be significantly larger in the Morton order curve.approach not only provides a cheap clustering technique, but also an efficient datastructure for performing CFD operations faster. In this work, we employ an SFCfor to partition a multi-block finite volume mesh.3.1.2 Space-filling CurvesAn SFC is a continuous mapping of a d-dimensional space onto a one-dimensionaldiscrete parametrized curve. The objective of this mapping is to retain spatial prox-imity of points in the linear space, i.e. neighbouring points in the physical domainremain close on the SFC. The central operation involved in the SFC partitioningapproach is the reordering of mesh elements based on a locality index assigneduniquely in each SFC. The ordered curve can then be partitioned into equal sec-tions wherein spatial locality relies on the accuracy of the original mapping func-283. Spatial Partitioning Algorithmtion. The clustering properties of SFCs inherently ensure a reasonable shape andlocality of clusters [82]. Also, in contrast to iterative or recursive partitioning al-gorithms, an SFC directly obtains a mapping for a given mesh. Additionally, thelinear mapping provides a method to access mesh elements for performing SFCoperations quickly. Due to these favourable properties, SFC have become an in-creasingly popular strategy for mesh partitioning.Among the various, well-documented SFCs, the Hilbert curve and the Mortonorder curve are standard choices in CFD applications [77, 83]. While the Hilbertcurve is known for its excellent proximity preserving behaviour [84], the Mortonorder curve is widely accepted for its simplicity and adaptability [85, 86]. Fig. 3.1illustrates the differences between the two mappings shown for grids with progres-sively increasing refinements. A partition of the entire domain is generally notaccessed completely by a contiguous part of the Morton order curve. The leapscharacteristic to this curve occasionally span the entire physical domain. On theother hand, the Hilbert curve retains measurable and bounded cluster locality [87],subjectively observable in the illustration. However, construction of Hilbert curvesis challenging, particularly for higher dimensional problems [88?90]. Moreover,Hilbert curves also contain leaps when used in arbitrary geometries. In the presentwork, we employ the Morton order curve although the partitioning algorithm de-veloped may be used in conjunction with any SFC.The Morton order curve is widely known for its compact construction rule [91].A locality index, called Z-value, is generated for every coordinate by interleavingthe bits in corresponding binary representations [91]. The Z-value of a coordinatein a 2-dimensional space is, hence, calculated asZ(xi,yi) = x1i y1i x2i y2i . . .xmi ymi (3.1)where x ji and y ji represents the jth significant bit of x and y coordinates of the ithpoint in its m bit representation. The interleaving concept can be extended to d-dimensions in general, at the expense of an increased frequency and size of thecharacteristic leaps between adjacent points [92]. In addition, unstructured gridscan be mapped, but bit-interleaving of the representative floating-point coordinatesis challenging. Connor and Kumar [93] introduced a Morton order curve construc-293. Spatial Partitioning Algorithmtion mechanism for floating-point coordinates which avoids assigning the Z-valueby sorting points based on their relative ordering. The algorithm is extendible tomultiple dimensions and easily implemented, but the ordering has been observedto be sensitive to coordinate values. We adopt the floating-point Morton orderingalgorithm after translating block-centre coordinates to avoid negative values. Thespatial partitioning algorithm presented uses properties of the generated curve toobtain highly localized partitions.(a) Mesh geometry. (b) Morton ordercurve.(c) Obtained clusters.Figure 3.2: Equi-sized partitioning of complex or refined mesh structuresresults in poor spatial locality within clusters. Regions of similar shading areassigned the same cluster.3.2 TheorySFC-based partitioning algorithms generally divide an SFC into equi-sized sectionswhich represent point-clusters in physical space. In Fig. 3.2, an irregularly refinedgrid (Fig. 3.2(a)) is used to illustrate partitions obtained for complex geometry ora refined mesh by dividing the Morton order curve (Fig. 3.2(b)) into sections ofequal size (Fig. 3.2(c)). Since two adjacent points on the Morton order curve maynot be spatially local, equi-sized partitioning may allocate these points to a singlepartition regardless of the size of the leap (distance in the spatial domain) betweenthe points. Such losses in partition quality may be acceptable for problems whereinspatial locality is not essential to the solution. However, in the context of our re-active flow problem, this would imply that statistical homogeneity of conditionallyaveraged scalars is not retained within a cluster containing a leap. Hence, the focus303. Spatial Partitioning AlgorithmAlgorithm 3.1: Spatial partitioning algorithm for multi-block meshes.Require: Morton ordered d-dimensional point set P of size nb.Ensure: C contains k sequential cluster sizes c j in point set P.1: function CLUSTER(set P, int k)2: R?{ri}|ri ? |~pi?~pi+1|? i ? {0,1, . . . ,nb?2}3: j? 1; c0? 1; c j? 0? j ? {1,2, . . . ,k?1}4: rth? THRESHOLD(R,k) ? using Algorithm 3.25: for all i = 0? nb?1 do6: if ri/rmax < rth then7: c j? c j +18: else9: j? j+1; nc? j10: end if11: end for12: C? BALANCE(C,R)13: C? PARTITION(C,k)14: return C15: end function16: function PARTITION(set C, int k)17: while nc < k do18: nc? nc+1; j0? j |A j = MAX(A)19: for all j = nc? j0+1 do20: c j+1? c j; j? j?121: end for22: c j, c j+1? c j/2 ? halve largest cluster23: end while24: return C25: end functionof our spatial partitioning algorithm is to avoid such locality losses by intelligentlypartitioning an SFC .The guiding principle of our partitioning heuristic is to partition the SFC be-tween points located across a leap. As presented in Algorithm 3.1, a normalizationscheme is employed to identify leaps between consecutive points on a Morton or-der curve. The distance (ri,i+1) between the ith point and its next adjacent pointis normalized by the maximum distance (rmax) between consecutive points on the313. Spatial Partitioning AlgorithmAlgorithm 3.2: Threshold consecutive point distance for partitioning.Require: Array R of consecutive point distances in Morton ordered set.Ensure: rth is maximum threshold distance that generates nc < k.1: function THRESHOLD(set R, int k)2: rmin? MIN(R), rmax? MAX(R); rstep? rmin/(2? rmax)3: nc? 1; rth? 14: while nc < k do5: for all i= 0? n?2 do6: if ri/rmax > rth then7: nc? nc+18: end if9: end for10: rth? rth? rstep11: end while12: return rth+ rstep ? generate nc < k13: end functioncurve asr?i =ri,i+1rmax. (3.2)The largest leaps are identified first and the curve is partitioned at those lo-cations. Using the step size defined in Eq. 3.3, progressively smaller leaps areidentified, creating more partitions. The Morton order curve has several identicallyspaced consecutive points throughout the mapping due to its repetitive Z-like pat-tern. Therefore, partitioning the curve at all leaps which are larger than a certainthreshold distance (rth) results in more clusters than the required number (k). Sincethis would add an additional challenging task of merging neighbouring clusters,Algorithm 3.2 provides a method to obtain the maximum possible threshold dis-tance which results in fewer clusters than required. The exact number of divisionscan then easily be produced by halving the largest clusters.rstep =rmin2 ? rmax(3.3)The partitioning algorithm has been illustrated in Fig. 3.3 for a Morton order curvethat has been constructed for a geometry identical to Fig. 3.2. The difference in323. Spatial Partitioning Algorithm(a) Morton ordercurve.(b) Partition at leaps. (c) Obtained clusters.Figure 3.3: Partitioning of the Morton order curve between points locatedacross leaps increases cluster locality at the expense of losing cluster sizebalance.spatial locality of obtained clusters is obvious. A quantitative analysis of a Mortonorder curve constructed for a simple geometry is performed in Section 3.4. Arange of point-clouds representing relevant geometries (described in the followingsection) were partitioned using the algorithm.Normalization of distances between adjacent points serves to maintain gener-ality of the algorithm. However, a strict adherence to this metric for the identi-fication of leaps ignores the size of clusters obtained. Often such a partitioningalgorithm produces exceedingly wide range of cluster sizes, particularly for irreg-ularly spaced points, resulting in a severe load imbalance. As a method to limitcluster size imbalance, Algorithm 3.3 merges clusters smaller than a threshold size(cth) with their neighbouring clusters. A cluster is merged with its closest adjacentcluster only if the inter-cluster distance is smaller than a pre-defined value.Discrete SFCs are not uniquely defined. In general, additional criteria are re-quired to specify the curve that best preserves locality. Appropriately defined mea-sures of spatial locality allow the construction of curves that minimize the metricin each cluster. Additionally, this provides a tool for meaningful comparison be-tween two SFCs. The measures defined by Perez et al. [94] and by Gotsman andLindenbaum [87] calculate the weighted sum of distances within a cluster, wherethe weights are inversely proportional to the spatial distance between the points.Another formal analytical measure for parallel domain decomposition using SFCshas been established for Hilbert curves by Tirthapura et al. [95]. A universal metric333. Spatial Partitioning AlgorithmAlgorithm 3.3: Cluster coarsening for reviving computational load balance.Require: Cluster size set C with nc elements. Point-to-point distance set.Ensure: C is devoid of small clusters that skew load balance.1: function BALANCE(set C, set R, int k)2: i? c0; cth? np/(2? k)3: rmax? MAX(R) rco? rmax/24: for all j = 1? nc?1 do5: i? i+C j6: if c j < cth then7: if ri/rmax < rco or ri?C j/rmax < rco then8: jco?{ri < ri?C j ? j : j?1}9: C jco?C jco+C jco+110: for all k = jco+1? nc do11: Ck?Ck+112: end for13: nc? nc?1; j? j+114: end if15: end if16: end for17: return C18: end functionfor the Morton order curve is more challenging because of the irregular frequencyand length of leaps, particularly in unstructured meshes.We define a simple metric to quantify the locality preserved within clustersobtained using our spatial partitioning algorithm. Within a cluster of points, theroot-mean-square distance (rrms) from the centroid (~rc) is a common method ofrepresent the locality preserved for a given number of points. The maximum dis-tance (rmax) between any points within the cluster is an estimate of the skewness.We define an index for each cluster as the ratioLR =rmaxrrms, (3.4a)343. Spatial Partitioning Algorithmwith the RMS distance written asrrms =?????nb?i=1|~ri?~rc|nb, (3.4b)where vectors~ri represent block-centre coordinates and nb is the number of blocksin the cluster. The ratio LR quantifies how skewed a given cluster is from the idealspherical distribution of points considering the number of points in the cluster.In our locality analysis, we use LR to compare the performance of partitioningmethods. A similar ratio (LR,c) is applied to quantify inter-cluster locality using thecoordinates of cluster centroids instead of block-centres.3.3 Test GeometriesThe partitioning algorithm has been tested on three different geometries pertinentto turbulent reactive flow simulations. The experimentally studied cases comprise aBunsen burner [96] and a bluff-body burner [97]. The third case is the combustionchamber in a Ricardo Hydra internal combustion engine, including a relativelysimple piston [98]. Details of each geometry including prominent geometricalfeatures have been listed in Table 3.1.Specification Dimensions (m)Case Type ?x ?y ?z Geometry FeatureGulder Burner 5 ?10?2 5 ?10?2 0.1 -Cambridge Burner 8 ?10?2 0.15 0.25 bluff-bodyRicardo Hydra SI-Engine 8 ?10?2 4 ?10?3 0.21 piston-bowlTable 3.1: Specification and characteristic dimensions of test cases.Computational meshes for the two burner cases were constructed using an in-house multi-block meshing code. An axi-symmetric mesh representing the RicardoHydra engine geometry (at 330? crank angle) was constructed using blockMeshutility of OpenFOAM [99]. Details of all three meshes have been summarised inTable 3.2. These meshes form a broad range of point distributions suitable fora rigorous testing of the partitioning algorithm. All tests were performed on a353. Spatial Partitioning Algorithmdesktop computer with 64-bit Intel R?Xeon R?CPUE31225 with 3.10GHz? 4 coresand 8GB memory.3.4 ResultsWe analyse the performance of our spatial partitioning algorithm primarily in thecontext of load and locality. Qualitative illustrations have been used to highlightthe inaccuracy involved in equal sized partitions along with an assessment of com-putational expense. Finally, validation using an a priori case has been reported.As an illustration of the partitioning mechanism employed, properties of theMorton order mapping have been plotted in Fig. 3.4 and Fig. 3.5 for various stagesof Algorithm 3.1. A simple cylindrical mesh geometry with 3200 blocks repre-senting Bunsen burner is being partitioned here to obtain 128 localized clusters.Fig. 3.4(a) graphs the normalized distance between consecutive points (which rep-resent blocks in this geometry) on the Morton order curve. The conspicuous peaksdenote leaps between adjacent points ? the largest leaps have unity normalizedlength by definition. 125 point-clusters are obtained on partitioning the curve at alllocations where the consecutive block distance is greater than the threshold (rth).In Fig. 3.4(b) the inter-cluster separations (which are the lengths of correspondingleaps) have been plotted or each cluster. Whereas equi-sized partitions would con-tain 25 blocks each, partitioning at leap positions produces some clusters (groupsof closely placed points in Fig. 3.4(b)) with fewer than 10 blocks. Fig. 3.5(a)shows the result of the coarsening procedure implemented using Algorithm 3.3.Closely placed clusters below a threshold distance are merged together, decreasingthe number of clusters to 101. Finally, k = 128 clusters are obtained by halvingthe largest clusters. In Fig. 3.5(b), the troughs represent regions where halving wasMesh Type Mesh SpecificationsCase Geometry Case Blocks CellsBunsen cylindrical block-based 3,200 1.63 ?106Bluff-body cylindrical block-based 588 4.7 ?105Ricardo Hydra axi-symmetric cell-based 1 2.1 ?103Table 3.2: Mesh type and features for different test geometries.363. Spatial Partitioning Algorithm0 800 1600 2400 320000. indexConsecutiveblockdistance,r?i(a) Separations between adjacent blocks on the Mortonorder curve.0 32 64 96 1280. indexConsecutiveclusterseparation,r? i(b) Separations between adjacent clusters of blocks onthe Morton order curve.Figure 3.4: Distances between consecutive points (block-centres) on the Mor-ton order curve displayed for the first two steps of the partitioning algorithm.The normalized distance between points (r?i ? ri/rmax) has been employed forthe case with nb = 3200, and k = 128.373. Spatial Partitioning Algorithm0 32 64 96 1280. indexConsecutiveclusterseparation,r? i(a) Separations between adjacent clusters after themerging procedure.0 32 64 96 12800. indexConsecutiveclusterseparation,r? i(b) Resultant adjacent cluster separations after halvingthe largest clusters.Figure 3.5: Distances between consecutive points (block-centres) on the Mor-ton order curve displayed for final steps of the partitioning algorithm. Thenormalized distance between points (r?i ? ri/rmax) has been employed for thecase with nb = 3200, and k = 128.383. Spatial Partitioning Algorithmdone ? clusters that are halved generally have high spatial locality within them re-sulting in very small inter-cluster distance between the halved partitions. Based onthe strictness of spatial-locality desired, the threshold size of cluster can be changedwhich is expected to affect coarsening results and hence, the cluster-size balance.Load balance is assessed by comparing resulting cluster sizes (ci) with sizesobtained from equal partitions (ceq,i). The ratioc?i ?ciceq,i(3.5)then provides a relative workload comparison for problems where each point con-tributes equally to the computational load. The locality metric (LR) defined inSection 3.2 has been used to quantify preserved spatial locality within each cluster.A similar locality metric (LR,c) estimates the packing of clusters within a geometry.Through Fig. 3.6 and Fig. 3.7, a comparison of the load balance and spatiallocality within clusters has been made between equi-sized partitioning and ourpartitioning algorithm for the simple Bunsen burner geometry. In addition, thecomparison is extended to manually partitioned axial slices obtain for the cylindri-cal mesh in Fig. 3.8(a). Relative cluster sizes obtained by the partitioning algorithmare roughly bounded between 0.5 and 2 with the exception of a few partitions con-structed for the sparse case, k = 32. Alongside, the packing of clusters has beengraphed in Fig. 3.8(b) using the inter-cluster locality metric. Equi-sized partition-ing results in a number of poor quality partitions particularly for k = 128.To visualize the poor quality of partitions often obtained from equi-sized SFCpartitioning, a section of a such partitioned Bunsen burner mesh has been shownin Fig. 3.9. Although several partitions may be mapped by a single contiguous partof the curve, Fig. 3.9(c) and Fig. 3.9(d) clearly illustrate that several partitions con-tain leaps. Such clusters are unfavourable in terms of minimizing communicationsince they comprise of two (or sometimes more) sub-clusters of blocks. A parallelcomputation (with each processor assigned a cluster) would, therefore, have highercommunication overhead for these clusters when compared to the cluster such asFig. 3.9(b). More importantly, for the model implemented in our turbulent reactiveflow simulations, these clusters would be entirely unacceptable since conditionalaverages of reactive scalars (such as temperature, for example) are expected to be393. Spatial Partitioning Algorithm0 32 64 96 12800.511.52Cluster indexClusterSize,c?i32 64 128(a) Cluster sizes from strictly equal partitioning of theMorton order curve.0 32 64 96 1282.533.544.5Cluster indexLocalitymeasure,LR32 64 128(b) Spatial locality of equi-sized clusters on the Mor-ton order curve.Figure 3.6: Degree of computational load balance for the Bunsen burner ge-ometry gauged using relative cluster (c?i) for strictly equal partitioning. Inter-processor communication is indicated by the locality metric (LR,c).403. Spatial Partitioning Algorithm0 32 64 96 12800.511.52Cluster indexClusterSize,c?i32 64 128(a) Cluster sizes obtained from the spatial partitioningalgorithm.0 32 64 96 1282.533.54Cluster indexLocalitymeasure,LR32 64 128(b) Spatial locality of corresponding clusters on theMorton order curve.Figure 3.7: Degree of computational load balance for the Bunsen burner ge-ometry gauged using relative cluster (c?i) for leap-based partitioning. Inter-processor communication is indicated by the locality metric (LR,c).413. Spatial Partitioning Algorithm0 2 4 6 8 10?10?22.533.544.5Axial distance, zLocalitymeasure,LRMP EP SPA(a) Spatial locality within each cluster.0 2 4 6 8 10?10? distance, zLocalitymeasure,LR,cMP EP SPA(b) Inter-cluster spatial distribution.Figure 3.8: Comparison of intra-cluster (LR) and inter-cluster (LR,c) localitymeasures for manual partitioning (MP), equi-sized partitioning (EP) and us-ing the spatial partitioning algorithm (SPA) on the Morton order curve in theBunsen burner geometry for k = 128 divisions.423. Spatial Partitioning Algorithm(a) Clusters obtained from equal par-titions of Morton order curve.(b) Cluster with no leaps.(c) Pathological cluster with leap. (d) Pathological cluster with leap.Figure 3.9: Clusters produced from the equi-sized partitioning of Morton or-der curve often contain leaps that reduce spatial locality within. Even for asimple cylindrical Bunsen burner geometry, leaps occur at frequent and irreg-ular intervals.completely different in the individually localized, but separated, sub-clusters. Es-sentially, for these two pathological cases, the underlying assumption made in theformulation of the combustion model would be violated.Results of leap-based partitioning have been illustrated using sections of par-titioned geometries in Fig. 3.10 for the Bunsen burner, and in Fig. 3.13 for therelatively complex geometries accompanied with respective sectional views. Spa-433. Spatial Partitioning Algorithm(a) Cross-section view. (b) Perspective view.(c) Section of clusters obtain fromequi-sized partitioning.(d) Clusters as obtained from the spa-tial partitioning algorithm.Figure 3.10: A qualitative picture of the difference between regular equi-sized partitioning and our spatial partitioning algorithm applied to Bunsenburner geometry shown above.tial locality is well preserved in the axi-symmetric mesh despite a low number ofcells. This emphasises the weakness of space-filling curves, especially the Mor-ton order curve, in accurately preserving locality in 3-dimensional domains. In thegeometries used in our study, points on the Morton order curve designate centresof blocks. Due to the varying shape and size of each block, an SFC mapping isexpected to have lower quality for these cases than cell-based construction.Fig. 3.11 and Fig. 3.12 assess the performance of the partitioning algorithm for443. Spatial Partitioning Algorithm0 32 64 96 1280123Cluster indexClusterSize,c?i32 64 128(a) Cluster sizes obtained for the bluff-body burner ge-ometry.0 32 64 96 12822.533.54Cluster indexLocalitymeasure,Lrms,c32 64 128(b) Spatial locality of corresponding clusters on theMorton order curve.Figure 3.11: Computational load balance measures (c?i and Lrms,c) estimatedfor the bluff-body burner for different cluster divisions, k = 32,64 and 128.453. Spatial Partitioning Algorithm0 32 64 96 1280.511.52Cluster indexClusterSize,c?i32 64 128(a) Cluster sizes obtained for the Ricard Hydra enginegeometry.0 32 64 96 1282.533.54Cluster indexLocalitymeasure,Lrms,c32 64 128(b) Spatial locality of corresponding clusters on theMorton order curve.Figure 3.12: Computational load balance measures (c?i and Lrms,c) estimatedfor the Ricardo Hydra engine for cluster divisions, k = 32,64 and 128.463. Spatial Partitioning Algorithm(a) Cambridge bluff-body burner. (b) Section of burner showing clus-ters obtained from algorithm.(c) Ricardo Hydra engine. (d) Perspective of engine showingclusters obtained from algorithm.Figure 3.13: Sections of complex geometries clustered using the spatial par-titioning algorithm along with appropriate sectional views.the two relatively complex geometries, the bluff-body burner and the Ricardo Hy-dra engine combustion chamber. It is observed that on average, clustering qualitydecreases with progressively more challenging geometrical features in the domain.Computational load is, however, maintained within the [0.5,2] bounds barring afew exceptions. This is expected since, the floating-point Morton order algorithmis sensitive to domain irregularities. Moreover, it is observed that, an optimumnumber of partitions, k = kopt provide the best load balance for a given geome-try. Often in large-scale parallel computing, one is constrained to use a numberof processors that is a power of 2 to take maximal advantage of the available re-473. Spatial Partitioning AlgorithmSensitivity Parameters (m) Computational Expense (ms)Case rmax rmin k = 32 k = 64 k = 128Bunsen 0.03 0.0008 72.54 41.6 35.8Bluff-body 0.25 0.0019 7.1 7.2 7.61Ricardo Hydra 0.04 0.0002 33.9 28.9 27.9Table 3.3: Performance of spatial partitioning algorithm in different geome-tries and inter-point distance parameters defining the corresponding Mortonorder curve.sources. As such, the distribution of cluster-sizes for each case considered here hasbeen graphed for 32, 64 and 128 clusters. For serial processing architectures, anydesirable number of partitions (preferably, k ? kopt ) may be constructed for serialcalculations.The space-filling curve partitioning approach is sensitive to the coordinate val-ues provided since inter-point distances can significantly affect the algorithm. Sev-eral parameters such as the maximum and minimum inter-point distances have beentabulated in Table 3.3. In addition the run-times have been tabulated for differentcluster numbers. A consistent correlation of execution time with mesh parameterssuch as cell numbers and inter-point distances is not be observed. However, theoverall execution time is negligible, particularly for an operation that is expectedto be done only once on a simulation with no moving mesh and without any meshadaptation. On simulations with moving mesh or mesh adaptation, while this com-putational cost would also be negligible, it is likely that there would be considerableredistribution of data among nodes on a parallel machine; nevertheless, the com-putational overhead of recalculating the distribution would be very small and theredistribution is almost certain to be needed regardless of locality considerations,simply to maintain load balance. In comparison to graph-partitioning methods, thecomputational expense is low and encourages the use of SFCs for low to mediumend parallel processing architectures.483. Spatial Partitioning Algorithm3.5 ConclusionsA new algorithm for partitioning mesh for arbitrarily complex geometries has beenproposed that favours maximizing locality at the possible expense of increasedload imbalance. The method is based on the floating-point Morton order algo-rithm, which provides a computationally inexpensive alternative for the dynamicmesh partitioning in arbitrary geometries. Careful translation of negative coordi-nates is required prior to the mapping in order to avoid inconsistencies. Also, thelack of a unique Morton order curve for a given domain lends unwanted geometricsensitivity to the partitioning algorithm. Our spatial partitioning algorithm can beused in conjunction with the floating-point Morton ordering algorithm or with anyexisting SFC-based domain decomposition. Cluster locality is improved againststrictly equal partitioning based on the Morton order curve. The technique canbe applied to geometries relevant to industrial problems. An acceptable degree ofcomputational load balance is achieved for large-scale parallel processing environ-ments. The LR locality metric quantifies spatial locality within clusters and roughlyindicates inter-processor communication time.494 | Engine Simulations4.1 IntroductionLean premixed combustion is increasingly employed by modern engines in orderto reduce fuel consumption and NOx emissions. Design and development of ro-bust premixed combustion systems hinge on predictive computational tools that notonly reduce dependence on experiments but may offer additional information. Thechoice of computational approach is a product of the trade-off between solutionaccuracy and the associated computational time. In this regard, LES is graduallybecoming a computationally feasible approach for understanding cyclic variationsin spark-ignition engines. In the meantime, RANS continues to offer an inexpen-sive alternative and provide acceptable design solutions in a time-frame suitable tothe industry.Closure issues arise inevitably from both these solution strategies; in particu-lar, the closure of species turbulent fluxes (?u??i u??j ) and chemical source-terms (??k)is the objective of combustion modelling. Few models account for the complexi-ties involved in turbulent chemical transport such as counter-gradient transport inreactive flows [100]. Although models based on the gradient transport hypothesisare relevant only for specific premixed combustion regimes, these are widely useddue to their simplicity. Estimation of chemical source-terms has received greaterinterest due to its contribution to the overall accuracy of the solution [15, 32].Preprint to be submitted to Combustion Science and Technology.504. Engine Simulations4.1.1 Relevant BackgroundChemical closure for RANS can be broadly divided into three different approachesthat respectively address the non-linear chemical source-term using an algebraic,a geometrical or a statistical analysis. Algebraic models represent the source-termas a function of easily acquired quantities such as turbulent mixing. Such modelsinvolve inaccuracies and must be tuned for the behaviour of turbulent premixedflame in context. These models have been hitherto commonly preferred in thecontext of RANS [35]. However, more accurate models have been developed andgaining popularity in recent years.Geometrical analyses of turbulent premixed flames are based on the flameletassumption conceptualized by [36]. Turbulence-chemistry interactions are resolvedby describing the turbulent flame as an ensemble of locally steady, one-dimensionallaminar flames. Two widely used techniques comprise flamelet-based chemicalclosure: the Flame Surface Density (FSD) approach and the G-equation model. Inthe FSD approach, the flame is identified as a surface which can be convected, dif-fused, curved and strained by the turbulent velocity field [37]. A similar approachintroduces a flame wrinkling description using the ratio of flame surface and itsprojection in the direction of projection [38]. In an experimental investigation byVeynante et al. [101], the closures generally used for the flame surface densityequation have been debated due to their inaccurate representation of propagationand curvature terms. The accurate closure of the turbulent flux and propagation ve-locity of the flame surface and the effects of curvature and strain rate on the flamesurface continues to be a persistent challenge [15, 32].In contrast, the G-equation model borrows from the description of a laminarflame surface using a level set [39]. This concept eliminates the need for resolv-ing the flame-front structure, but the lack of resolution of the flame-front structuredisrupts the accuracy of chemical closure when the flame structure is adequately af-fected by turbulence. The thickened flamemodel was proposed by [41] as a methodto resolve the flame structure on the mesh. However, the empirical parameters usedto thicken the flame are cannot be applied to arbitrary flame specifications. As anadvancement, a three-zone coherent flame model based on a flame surface densityequation and a conditional averaging technique that allows precise reconstruction514. Engine Simulationsof local properties has been used with success in gasoline engines [102].As a less restrictive alternative to the geometric analysis of flame-fronts, thestatistical behaviour of a turbulent reactive flow is represented using a PDF of asingle variable that characterizes the one-point, one-time state of a turbulent re-active flow. In premixed flames, the reaction progress variable defined as c = 0in fresh gases and c= 1 in burnt gas is employed. One line of approach then is topresume the PDF of the progress variable and obtain the averaged chemical source-term by integrating the source-term conditionally averaged on the progress variablewith its local PDF as??k(~x, t) =? 10??k,lam(c?)P(~x, t;c?)dc? (4.1)where P(~x, t;c?) is the PDF and ??k,lam(c?) is the presumed conditional moment(PCM) written as a function of c?, the statistical random variable associated withthe reaction progress variable. The Bray-Moss-Libby (BML) model combines suchan approach using a bimodal distribution with a physical analysis that accounts fordifferences in fresh and burnt gases. Veynante et al. [101] observe that the BMLmodel has good overall trends but lacks precise closure of the flame wrinklinglength scale.Although variations of the progress variable approach exist, most such mod-els describe constant pressure premixed combustion; advancements to accuratelyrepresent varying pressure conditions have been relatively recent [103]. Recently,a model has been proposed for simulating turbulent premixed flames in the cor-rugated flamelet regime by combining the BML approach in a PDF frameworkwith elements of the FSD equation using a reaction progress variable [104]. Also,the Eulerian particle flamelet model has been applied to the simulation of non-premixed combustion in diesel engines [105] and bluff-body flames [106]. Thisapproach includes unsteady effects in stretched laminar flamelet models by intro-ducing marker particles that correspond to flamelet histories associated with thepath of the particle in the turbulent flow field. The range of flamelet PDFs requiredfor simulating a particular combustion process may differ from the ensemble cho-sen manually, thereby demanding an a priori knowledge of the solution. Experi-mental observations of [45] highlight that the flamelet assumption is unwarranted524. Engine Simulationsfor high turbulence intensity turbulent premixed flames.Various other statistical methods have been developed such as the broad classof transported PDF models developed for turbulent reactive flows by Pope [53].In these models, the complete statistical description of the state of the flow is pro-vided through a velocity-composition joint-PDF. A transport equation is writtenfor the joint-PDF and solved using Monte Carlo methods. The relatively recent,Multiple-mapping Conditioning (MMC) method combines features of CMC meth-ods, PDF methods and mapping closure models [54, 55]. These methods provideexcellent chemical closure; but, the associated computational expense discouragestheir application to industry-scale problems.Originally formulated by Klimenko [46] and Bilger [47], CMC has been wellestablished as a chemical closure tool for the simulation of non-premixed combus-tion [15, 32]. It has been used successfully for industrial applications includingspray combustion in diesel engines [48]. CMC offers a promising closure tool forsimulating premixed combustion since it is devoid of the assumption of a prevail-ing combustion regime. Instead of presuming the reaction rate that is conditionallyaveraged on the reaction progress variable, as in conventional flamelet based sta-tistical approaches, CMC calculates this term for the kth filtered species transportequation from the hypothesis??k|? ? ??k(T |? , Yk|? , ? |? ), (4.2)assuming small fluctuations of scalars about their conditional averages. The con-ditionally averaged reaction rate obtained is integrated in a similar fashion to thePCM approach by writing??k(~x, t) =? 10??k|c?(c?)P(~x, t;c?)dc? (4.3)where ??k|c?(c?) is estimated here through the CMC hypothesis. As opposed tothe flamelet method, CMC has the capability of accurately modelling the chemicalsource-term for flames beyond the flamelet regime [33]. However, transport equa-tions typically solved for conditionally averaged scalars consist of various unclosedterms which are not well understood and are modelled inaccurately [31].534. Engine SimulationsCSE is a novel CMC-based chemical closure model. In CSE, an inverse prob-lem is solved over spatially localized ensembles of finite volumes to obtain con-ditionally averaged scalars. Therefore, CSE is devoid of the closure issues asso-ciated with conventional CMC methods. Recent advances have allowed CSE tobe applied to RANS simulations of premixed combustion in a Bunsen burner [58].However, CSE has never been applied to practical problems which are governedby transient phenomena and characterised by complex geometries. The followingsections present the theory and results for application of CSE to natural gas fuelledspark-ignition engines.4.2 Theory4.2.1 Conditional Source-term EstimationThe general concept of CSE follows from the CMC method wherein the condi-tionally averaged source-term is approximated by conditionally averaged scalarsrelevant to the chemistry model from the CMC hypothesis (Eq. 4.2). The averagedchemical source-term is calculated by integrating the corresponding conditionallyaveraged source-term with the presumed PDF at the given spatial coordinate withrespect to the conditioning variable as follows??k =? 10??k|c?P(c?)dc?, (4.4)where c? is the statistical random variable associated with the conditioning vari-able. The presumed PDF for any point in the domain can be estimated with rea-sonable accuracy using the first two moments of the conditioning variable. Asdescribed in the following section, a normalized oxygen based reaction progressvariable has been used in this work with the presumed PDF shape chosen as the ?function. While conventional CMC models solve transport equations for the condi-tionally averaged scalars that are necessary for obtaining the conditionally averagedreaction rates using the CMC hypothesis [31], CSE employs an integral equationsimilar to Eq. 4.3. Each Favre-averaged scalar (? ) necessary for calculating theconditionally averaged reaction rate can be written in terms of the corresponding544. Engine Simulationsconditionally averaged scalar as??k(~x) =? 10?k|c?(~x,c?) P?(~x,c?)dc?, (4.5)for the corresponding spatial coordinate (~x) in the physical domain. The condition-ally averaged scalar can then be averaged over an ensemble of discrete points (A)in the domain. If the ensemble of discrete points is chosen in a way to preserve sta-tistical homogeneity of conditionally averaged scalars, then the ensemble averageis independent of the spatial coordinate and the integral equation is written as??k(~x j) =? 10?k|c?(c?) P?(~x j,c?)dc?, j ? A (4.6)for each point j within the ensemble. Typically, conditionally averaged scalarsare statistically homogeneous over spatially localized ensembles of points [56].In this work, however, the entire computational domain is considered as a singleensemble of points where the above equation is assumed to be accurate. With theconditionally averaged scalar now invariant of the jth spatial coordinate (~x j) in theseensembles, the previous equation can be treated as a Fredholm integral equation ofthe first kind with P?(~x j,c?) as the kernel. The conditionally averaged scalar can beobtained for the ensemble from the deconvolution of this integral equation. Usingm bins for the conditioning variable, the integral equation is written in the discreteform as~b= A~?, A ji =? c2c1P?(c?;~x j)dc?, (4.7)where b j = ??k(~x j) and ?i = ?k|c?i is the conditional average in the ith bin. Inte-grations of the presumed PDF required for calculating A ji are done a priori andtabulated for m bins using lower (c1) and upper (c2) limits for the ith bin of condi-tioning variable. As in the work of Salehi et al. [58], the Tikhonov method is im-plemented for regularising the solution of the integral equation using an unstrained1-dimensional laminar flame solution appropriate for the mixture properties andinitial conditions. Fig. 4.1 shows different conditionally averaged temperaturesobtained from laminar solutions of mixtures with varying relative air-fuel ratio (? )used as regularisation solutions. An L-curve approach is employed for the choice554. Engine Simulations0 0.2 0.4 0.6 0.8 1600800100012001400160018002000220024002600Smaller ?Conditioning variable, c?ConditionedTemperature,T|c?(K)Figure 4.1: Laminar flame solutions of mixtures with varying relative air-fuelratio (? = 1,1.1, . . . ,1.5) used for the regularisation of CSE solutions.of an optimal regularisation parameter as discussed by Salehi et al. [58]. In thiswork, the entire physical domain is treated as a single ensemble of points for theinversion procedure.4.2.2 Presumed Probability Density FunctionThe PDF is central to the accuracy of CSE related integral inversions and the choiceof conditioning variable is important for the overall description of the solution.Swaminathan and Bilger [107] have discussed various choices appropriate for tur-bulent premixed flames. It is critical for the progress variable to vary linearly acrossthe turbulent flame so that the computational description is accurate. In this work, areaction progress variable defined using O2 mass fraction is used due to the roughlylinear transformation of oxygen across an ignition process. The reaction progressvariable is defined asc= Y0R ?YRY 0R ?Y?R(4.8)where O2 is the chosen reactant (R), Y?R is its equilibrium mass fraction, and Y0Ris the initial mass fraction of oxygen in the premixed mixture. The ? function564. Engine Simulations0 0.2 0.4 0.6 0.8, c?Variance,c???2(a) Map of means and variances of stored PDFs.0 0.2 0.4 0.6 0.8 variable, c?PDF,P(c? ,c??,c???2c?= 0.1, c???2 = 0.1c?= 0.9, c???2 = 0.1c?= 0.5, c???2 = 0.9c?= 0.5, c???2 = 0.25(b) Sample distributions of the ? -PDF.Figure 4.2: The ? -PDF has been tabulated a priori for different means andvariances. A sample distribution of the two presumed PDFs employed hasalso been illustrated.574. Engine Simulationshas been used as the presumed shape of the PDF despite its failure in reproducingthe true shape of the PDF for premixed flames as discussed by Jin et al. [61]. Asa first implementation of the CSE framework in an industrial context, the ? -PDFoffers an acceptable ad hoc solution and recovers the extreme properties expectedof the true PDF such as having ? functions appear at 0 and 1 for maximal varianceor a single ? appear at the mean for zero variance. In practice, the cumulativedistribution function is used instead of a PDF because of its monotonic propertieswhich allows better numerical accuracy.Given the reaction progress variable, its mean and variance characterise thepresumed PDF in an approximate sense. Therefore, the one-time, one-point PDFis calculated asP(c?;~x, t)? P(c?; c?, c???2) (4.9)where c? is the statistical random variable associated with the progress variable c?and c???2 are the mean and variance of the progress variable. In the a priori tabulationprocess, the binning is performed as(Pdc)1 =? ?c/20Pdc ; (Pdc)m =? 11??c/2Pdc (4.10a)for the first and last bins, and as(Pdc)i =? ci+?c/2ci??c/2Pdc (4.10b)for the general ith bin. Careful tabulation of the first and last bins is done usinghalf-bins in order to preserve the precise placement of the ? functions at 0 and/or1. Fig. 4.2 shows the mean-variance map for which the ? -PDF is tabulated prior tosimulations; several sample distributions of the ? -PDF have been shown alongsideto demonstrate the range of behaviour described by ? -PDF.4.2.3 Chemistry ReductionThe conditionally averaged reaction rate is calculated using Eq. 4.2 from the knowl-edge of conditionally averaged scalars obtained via the inversions discussed previ-ously. In general, it is possible to consider the entire state space of the turbulent584. Engine Simulations0 0.2 0.4 0.6 0.8 101020304050Progress variable, c?Conditionalreactionrate,?? c|c?(s?1 )10 bar30 bar25 barInterpolatedFigure 4.3: Conditional reaction rates at 25 bar obtained directly comparedwith an interpolation between tabulated values at 10 bar and 30 bar.reactive flow involving detailed chemistry into the chemistry model necessary forthe CMC hypothesis. In such a case, CSE can be used to calculate all requiredconditionally averaged scalars which can then be used to obtain the conditionallyaveraged reaction rate. However, in practice the chemistry model is simplified byreducing the chemistry space in order to decrease computational expense. Further,based on the gradients involved in the reaction rate dependence with a scalar, theexact number of bins may be varied. In this work, the Trajectory Generated Low-dimensional Manifold (TGLDM) proposed by Wang et al. [62] is implemented totabulate a modified GRI-Mech chemistry developed by Huang et al. [108] usingtemperature, reaction progress variable and pressure. Due to the near-linear depen-dence of reaction rate on average pressure at any given point, only 3 manifolds ofpressure (10, 30 and 60 bar) were tabulated; any conditional reaction rate is wellapproximated from an interpolation between the available manifolds as illustratedin Fig. 4.3. Therefore, conditionally averaged reaction rate is expressed as??k|c?(c?)? ??k(P, T |c?(c?)) (4.11)594. Engine Simulationswhere P is the average pressure and T |c? is the conditionally averaged temperatureobtained via the CSE model. While 51 bins of the reaction progress variable areused (including the first and last half-bins), 100 bins of temperature have been useddue to the sensitivity of reaction rate to temperature. Since the tabulation is doneusing the reaction progress variable, the conditioning variable is known for eachbin. As a result, the only unknown parameter required for closure of chemicalreaction source is the conditional average of temperature. This is in contrast to thework of Salehi et al. [58] wherein a Bunsen burner flame was simulated using CO2and H2O mass fractions to tabulate the chemistry with normalized CO2 fractionused as the conditioning variable.4.2.4 Governing EquationsAs discussed earlier, the CSEmethod employs a presumed shape for the PDF of theconditioning variable. Construction of the simple ? -PDF relies on the knowledgeof mean and variance of the conditioning variable at the given point and time.These quantities can be obtained using transport equations: since the conditioningvariable is normalized oxygen mass fraction, the equation is identical to that writtenfor any species mass fraction:? ?? c?? t +? ?? u?ic??xi=?? ?? u???i c???xi+ ??xi(??Dc? c??xi)+ ??c, (4.12)where Dc is the Fick?s coefficient of molecular diffusivity for the species related tothe progress variable. Two unclosed terms exist in the equation mentioned above:i) the turbulent scalar flux, u???i c?? and ii) the chemical source-term ??c which is es-timated using CSE. The widely applied method of modelling the turbulent scalarflux using gradient transport hypothesis has been employed in this work. This isapplicable for flows where counter-gradient diffusion is not dominant [100]. Atransport equation for the variance of progress variable is also solved:? ?? c???2? t +? ?? u?ic???2?xi= ??xi(?? ?TSc1? c???2?xi)+2?? ?TSc2? c??xi? c??xi?2?Dc ?c???xi?c???xi+2c????c604. Engine Simulationswhere ?T is turbulent viscosity. Similar to the work of Salehi et al. [58], all Schmidtnumbers are set to 0.7 based on the work of Yimer et al. [109]. The two unclosedterms in the equation namely, the correlation between source term and fluctuationsof the conditioning variable (c????c) and the Favre-averaged scalar dissipation rate(?D ?c???xi?c???xi)have been modelled in an identical fashion as used for the RANSsimulation of Bunsen burners by Salehi et al. [58]. Vervisch et al. [110] explainthese models at great length; the discussion is beyond the scope of this work.In comparison to a steady state Bunsen burner flame, an SI engine consists ofwalls on the domain boundaries. Heat transfer from these walls contributes signif-icantly to the governing enthalpy equation and dictates the decay of flame speednear the wall. In this work, walls have been maintained at constant temperature(T = 500K). The overall wall heat transfer is calculated using an effective turbu-lent thermal diffusivity and the surface normal gradient of enthalpy. However, theisothermal assumption disrupts the validity of conditional averages calculated overa single domain-wide ensemble of points. Even though the chemical composi-tion may change and result in temperature fluctuations, the assumption suppressesany such deviations close to the wall. A more accurate method of calculation ofthe scalar conditional averages would be to construct ensembles of points that areequi-distant from the wall. This method is avoided here since in an axi-symmetriccase, with few finite volumes, this method inevitably results in a rank-deficient ker-nel for the inverse problem; a fine mesh simulation in a 3D domain is appropriatefor such ensemble construction methods.A summary of the structure of CSE-TGLDM model has been illustrated inFigure 4.4. The structure shows the flow of control in the CFD solver which origi-nally provides the mean and variance of conditioning variable through the solutionof transport equations. Once the ? -PDF is constructed for all points in the meshfrom the corresponding mean and variance values, the kernel containing all PDFsis passed along with necessary scalars (Favre-averaged temperature and pressure)calculated in the previous time step to the CSE routine. The desired regularisationprocedure is employed to perform inversions in the CSE routine to calculate theconditional averages of scalars. Since the conditioning variable is an oxygen basedreaction progress variable, also used in the TGLDM tabulation, the conditional av-erages of this species are already known. All requisite conditional averages are then614. Engine SimulationsCFD??k, T?CSE??k??k|c?k|c,T |cc?, c???2PDFTGLDMTGLDMPDFFigure 4.4: Iterative algorithm of the CSE-TGLDM chemical closure model.used to retrieve conditionally averaged reaction rates from the TGLDM chemistrytabulation; this vector of values is integrated with the PDF at each point to obtainthe average reaction rate in the corresponding cell.4.2.5 Spark-Ignition ModelFor a mixture with given relative air-fuel ratio (? ), the experimentally observedinitial conditions are matched by using the reported pressure at intake valve clos-ing IVC. Another variable important for matching combustion relevant conditionsis the timing of spark-ignition. A realistic high-temperature charge plasma is chal-lenging to implement in the context of numerical simulations because placing atiny zone of extremely high temperature is constrained by mesh resolution and tab-ulation of chemical kinetic constants and species chemistry at the correspondingscale of temperature. Thiele et al. [111] have noted that due to the fast expansionof the plasma channel, a complicated flow-field is developed after the emission ofa shock wave. Development of a propagating flame after the initiation of the flamekernel is dominated by reactive and diffusive processes [111, 112].Pischinger and Heywood [113] derive a first-law-based model for the growthof spark-generated flame kernel. They focus on the effect of heat losses on the rateof flame kernel development and observe significant cycle-by-cycle variations. Tanand Reitz [114] explain that the commonly used method of increasing internal en-624. Engine Simulationsergy in specified ignition cells during ignition by a factor each time step is verysensitive to the computational mesh size. Instead, they derive an equation to calcu-late spark-ignited kernel growth rate by considering the effects of the spark ignitiondischarge energy and flow turbulence on the ignition kernel growth. In a more re-cent study, Enaux et al. [115] advance an Arc and Kernel Tracking Ignition Model(AKTIM) using a three phase progress variable tracking procedure for simulat-ing the kernel development. Another detailed physical model developed by Yas?ar[116] accurately describes the momentum and energy exchange between the gasand plasma.A simpler approach referred by Mastorakos [117] and explained in detail inmany texts including that by Law [118] is adopted for this first implementation ofCSE in an industrial context. Typically, ignition is forced in the uniform mixtureby a sudden insertion of a layer of burnt gas. In this work, instead of initiating aspark at the ignition timing, we place a developed progress variable profile (??c) in asmall volume (typically, ? 2.5% of the chamber volume) around the region of sparkplug. This developing kernel is made highly reactive by setting the mean (c?= 0.5),variance (c???2 = 0.25) and temperature (T ? 0.5Tad) to appropriate conditions. Theradius of the flame kernel is then calculated from2piR3kernel3 = ??V? (4.13)where ? is the fraction of chamber volume occupied. Special care must be taken inchoosing ? since a small kernel may not be able to withstand the turbulent straininduced in the engine at the corresponding timing and may, therefore, lead to ex-tinction of the kernel. The timing of kernel placement is decided by the displace-ment required to achieve 5% fuel mass fraction burned (MFB) at the same timingas reported in the experimental study of Reynolds [98].634. Engine Simulations81.1257.920.810.33R6.35Dimensions in mmFigure 4.5: Bowl-in piston geometry of the Ricardo Hydra engine.4.3 Problem Specifications4.3.1 Experimental Validation CaseThe engine used for experimental validation is a Ricardo Hydra single-cylinder re-search engine rebuilt for natural gas fuelling. The piston geometry is a modifiedform of the standard Ford Festiva bowl-in piston with dimensions as specified inFig. 4.5. Other important specifications of the engine geometry have been listed inTable 4.1. Measurements have been provided with varying relative air-fuel ratios(? = 1,1.1, . . . ,1.5) of the premixed homogeneous natural gas-air mixture. Thetests were reported for different engine speeds, but we have made comparisons at2500 rpm, in particular for a lean mixture with ? = 1.5. The natural gas compo-sition reported in the experimental study was scaled to provide a fuel compositionthat matched the chemistry tabulation employed.In the experimental work, heat release measurements calculated from the in-cylinder pressure trace measured using a pizeo-electric pressure transducer were644. Engine SimulationsSerial Geometrical Feature Value Unit1 Number of cylinders 1 ?2 Bore 81.12 mm3 Stroke 88.90 mm4 Connecting rod length 158.01 mm5 Swept volume 459.46 cm36 Clearance volume 41.988 cm37 Compression ratio 11.94:1 ?12 Inlet closes 56 ?ABDC13 Exhaust opens 56 ?BBDCTable 4.1: Specifications of the Ricardo Hydra research engine.used to estimate combustion duration and the fuel Mass Fraction Burned (MFB) asa function of crank angle. No turbo-charger is employed in the experiments; how-ever, reported pressures have large values in the low-pressure zones away from theTop Dead Center (TDC). These have been attributed to the inaccuracies involved inmeasuring low pressures using a piezoelectric sensor. A corrected intake manifoldpressure (PIVC = 2.32 bar) has been used in our calculations in order to match theearly pre-combustion pressure trace.The primary emissions examined in the experimental study were the brake spe-cific exhaust emissions of NOx, tHC and CO; in our work, NOx emissions havebeen reported. Experimentally observed emissions were reported as the mass flowrate of exhaust emissions normalized with the brake power produced by the engine:NOx(g/kW ?hr) =NOx(g/h)Pb(kW )(4.14)where Pb is the brake specific power generated. In our study, the emission massoutput has been normalized with the Gross Indicated Work (GIW) generated.4.3.2 Numerical SolverSimulations have been performed on a 2-dimensional axi-symmetric mesh shownin Fig. 4.6 generated using the blockMesh utility of the OpenFOAM CFD software[99]. The cell sizes have a gradient to provide a fine mesh structure near the spark654. Engine Simulationsignition region and the entire mesh was refined at various pre-specified intervals(crank angle, CA = 330,390,420) to avoid cell skewness as much as possible. Amaximum Courant number,Cmax = 0.25 was used with adjustable time step startingwith an initial value of 0.05 s. Solution was written at unity intervals of crank angle.Limited Gauss linear schemes were used for the solution of Laplacian and di-vergence operators in transport equations of species, enthalpy and momentum. TheGauss upwind scheme was used for the solution of transport equations relevant tothe k? ? turbulence model used from the OpenFOAM library.Figure 4.6: Axi-symmetric mesh constructed at 270? crank angle.4.4 ResultsThe following section describes our results; flame kernel development has beenhighlighted and the pressure map obtained has been used to calculate heat releaserates and fuel mass fraction burned. Further, trends of major pollutants, NOx andCO have been discussed for mixtures of varying ? .664. Engine Simulations345 345.5 346 346.5 347051015Crank angleMassFractionBurned,%?? 1.5%?? 2.5%?? 3.5%(a) Kernel growth from initiation to t5% m fb.345.5346346.5 0 0.250.5 0.75170010001300160019002200Crank Angle c?T|c?(b) Conditional averages post kernel initiation.Figure 4.7: Development of flame kernel shown as the degree of MFB until5% MFB and corresponding conditionally averaged temperatures obtained fora lean mixture with ? = 1.5 at 2500 rpm.674. Engine Simulations(a) Mean of progress variable, c?(b) Variance of progress variable, c???2Figure 4.8: Mean and variance contours of progress variable at 5% MFB fora mixture with ? = 1.5 at 2500 rpm.4.4.1 Flame Kernel DevelopmentUpon the imposition of the impulsive high reactivity field on the turbulent reactiveflow, quantities averaged over the geometry, such as the MFB, exhibit sensitivitytowards the size of the initial kernel. Fig. 4.7(a) shows the increase in MFB frominitiation of the flame kernel until 5% fuel MFB (t5%MFB) is attained in the com-bustion chamber for different initial volumes (as a fraction of chamber volume) ofthe kernel. This characterises the growth of the flame kernel and has different pro-684. Engine Simulationsfiles for varying mixture fractions. The profile of MFB growth has been used as amarker for judging whether the kernel will propagate as a flame. A tiny kernel doesnot provide enough energy sufficient for successful propagation [117] and is oftenextinguished by turbulent strain induced by the flow as outlined by Eichenbergerand Roberts [119]. In contrast, a large kernel may be characterised by rapid ker-nel growth and may require that the kernel be placed too close to t5%MFB therebyallowing little time for the evolution of kernel properties and conditional averages.For the case studied in Fig. 4.7(a), the experimentally observed displacement be-tween the spark and 5% MFB was approximately 31?. The effects of the spark onthe turbulent flow field from the initial discharge to the timing of kernel placementcannot be adequately modelled within the 2? displacement permitted by our empir-ical progress variable approach. Therefore, several measurable parameters in thestudy such as the timings of peak pressure and the associated crank angle cannotbe accurately estimated.Fig. 4.7(b) shows the variation of conditionally averaged temperature obtainedin the flame kernel development phase between initiation and t5%MFB. Significantvariations are noted from the steady laminar flamelet behaviour highlighting theeffects of turbulence on the flame. It is evident that the CSE model captures suchdeviations from flamelet behaviour of the flame which might be characterised bythe distributed reaction zone regime. In Fig. 4.8, contours of reaction progressvariable variance highlight the regions that roughly characterise the flame surface.Flame progress is observed to be slower close the engine walls as expected due toheat transfer away from the flame and possibly due to flame quenching.Kernel placement timings have been compared for a range of natural gas-airmixtures (varying ? ). These have been graphed along with the respective timingsof 5% MFB and of spark ignition. Fig. 4.9(b) compares the combustion duration(defined as the displacement between 5% and 95% fuel (MFB) trends observed ex-perimentally and using simulations. In addition to decreasing fuel content, the sizeof reactive flame kernel zones is also expected to affect the duration of combustion.Although a slight increase in the combustion duration is noted for leaner mixtures,the expected values are generally overestimated. To some extent, the errors may beattributed to the initial uniform progress variable profile of the flame kernel and itstiming which collectively result in an incorrect estimate of the initial kernel growth694. Engine Simulations1 1.1 1.2 1.3 1.4 1.5?40?30?20?10Relative Air-Fuel Ratio, ?Kernelplacementtiming,?BTDCIgnition2.5% Kernel5% MFB(a) 2.5% volume kernel initiation timing.1 1.1 1.2 1.3 1.4 1.55960616263Relative Air-Fuel Ratio, ?CrankdisplacementangleExperimentCSE(b) Combustion duration (5%?95% MFB).Figure 4.9: Initiation timing of high-reactivity flame kernel with ? 2.5%chamber volume for varying ? and their combustion durations.704. Engine Simulations270 300 330 360 390 420 450024?106Crank AnglePressure,(Pa)ExperimentCSEFigure 4.10: Pressure map for a lean mixture (? = 1.5) at 2500 rpm.rate. Moreover, the profile of progress variable and temperature in the initiatedkernel does not model the flame-generated turbulence within the kernel leading toan inaccurate growth rate estimate.4.4.2 Pressure TraceThe pressure trace has been obtained from intake valve opening (IVC) at 270? toexhaust valve opening (EVO) at 484? crank angle (CA). Fig. 4.10 shows the com-parison of pressure profile with values obtained experimentally. The piezoelectricpressure transducer used in the experiment by Reynolds [98] is susceptible to er-ror at low pressures. Therefore, the intake manifold pressure observed by [98] hasbeen used as the initial pressure at 270? CA instead in order to achieve accuratematching of the pressure profile prior to combustion. Despite kernel placementat the appropriate timing as discussed previously, calculated pressures are slightlyhigher until close to TDC. Severe fluctuations in the pressure map are observedpost TDC; these fluctuations are attributed to the skewness of moving mesh cellsat TDC.714. Engine Simulations270 300 330 360 390 420 450?0.500.511.52Crank AngleHeatReleaseRate,(J/?CA) ExperimentCSEFigure 4.11: Heat release rate for a lean mixture (? = 1.5) at 2500 rpm.Heat release rates (HRR) have been calculated using the pressure profile byemploying the following relation obtained using the first law of thermodynamicsand assuming uniform and identical properties of reactants and products:dQd? =1??1(VdPd? + ?PdVd?)(4.15)where ? is the ratio of specific heats for natural gas. The HRR profile obtained fromboth the experimental and simulation pressure trace were filtered using a low-passfilter in MATLAB followed by a smoothing function. Significant negative heatrelease is observed for the experimental profile in Fig. 4.11; this may be attributedto incorrect pressure values read by the transducer at earlier crank angles. Due tothe fluctuations of pressure near the TDC, calculation of heat release is prone toextreme error in this region.The technique used here for computing the mass fraction burned (MFB) isbased on the method developed by Rassweiler and Withrow [120]: the pressurerise due to combustion is estimated from the total pressure rise for small crankangle intervals as?P= ?Pv+?Pc (4.16)724. Engine Simulations270 300 330 360 390 420 45000. AngleMassfractionburned,%ExperimentCSEFigure 4.12: Variation of MFB for a lean mixture (? = 1.5) at 2500 rpm.where ?Pv and ?Pc are the change in pressure due to volume change and combus-tion respectively. MFB can then be calculated as a fraction of the pressure changedue to combustion up to a given crank angle with respect to the entire range fromthe Intake Valve Closes (IVC) to Exhaust Valve Opens (EVO):MFB= ???i ?Pc??EOC?i ?Pc(4.17)where ?i is the initial crank angle (270?) chosen to be close to IVC. Assumingpolytropic compression in each small interval of crank angle change, the pressurerise attributed to combustion can be calculated for the ith interval as?Pc,i = Pi?Pi?1(Vi?1Vi)n(4.18)where n is the polytropic index. The mass fraction burned for a lean mixture with? = 1.5 is shown in Fig. 4.12. It is to be expected that the profile is non-smoothbecause the pressure oscillates severely in the region of TDC.In addition to the pressure profile, the peak pressure and corresponding crank734. Engine Simulations1 1.1 1.2 1.3 1.4 1.555.566.577.5?106Relative Air-Fuel Ratio, ?Peakpressure,P max(Pa)ExperimentCSE(a) Maximum average cylinder pressure.1 1.1 1.2 1.3 1.4 1.581012141618Relative Air-Fuel Ratio, ?Crankangle,? ATDCExperimentCSE(b) Crank angle of maximum average pressure.Figure 4.13: Peak pressure and associated crank angle for varying ? .744. Engine Simulationsangles are calculated and have been compared with the experimental estimatesin Fig. 4.13. Although the pressure maxima show a decreasing trend similar toexperiments, the exact values are overestimated by the CFD calculation. On theother hand, the crank angle of peak pressure does not show a consistent trend as inexperiments perhaps due to small mismatches in the 5% MFB timing.A sensitivity test was performed to gauge the effect of intake manifold pressureon the simulated peak pressures and corresponding crank angles. Fig. 4.14 showsthe different maximum pressure values obtained for slightly different intake pres-sures. Matching the exact values from experimentally obtained pressure profilesmeans using erroneous results as opposed to noting the corrected intake manifoldpressure that is measured with higher precision in experiments. Similarly, crankangle of maximum pressure shows reasonably large variations with the choice ofintake manifold pressure. The mismatch between experimentally observed pres-sures at IVC and the pressure used for simulations (P= 2.32 bar) causes significantdeviations in estimated combustion duration and MFB among other parameters.4.4.3 EmissionsNOx pollutant emissions have been calculated using the chemistry mechanism pro-posed by Hanson and Salimian [121]. A transport equation for NO is solved and thechemical source-term is closed using the CSE method, an approach similar to thatused by Wang et al. [62]. Fig. 4.15 shows that CSE predicts the general decreas-ing trend of NOx calculated at EVO with decreasing fuel content in the homoge-neous charge. For an appropriate comparison with the experimental measurementof NO2, a scaling factor (44/30, based on the molecular masses of the two com-pounds) is applied to the NO estimates calculated by the numerical solver. Thesecalculations have been normalized with the gross indicated work (GIW) estimatedfrom the pressure obtained via simulations, whereas the experimental data consistsof calculations made in the exhaust flow normalized with the brake power. Thevalues of normalized emissions are only expected to increase (thereby, matchingwith experimental results better) if we were to account for friction in the GIW. Ex-perimental error in NOx mass measurement and exhaust flow rate are considered? error in brake power (Pb) measurement is not considered in the graphed error754. Engine Simulations2.1 2.2 2.3 2.4 2.5?10666.577.5?106Intake manifold pressure, (Pa)Peakpressure,P max(Pa)(a) Maximum average cylinder pressure.2.1 2.2 2.3 2.4 2.5?106368370372374376378Intake manifold pressure, (Pa)Crankangle,? ATDC(b) Crank angle of maximum pressure.Figure 4.14: Sensitivity of maximum cylinder pressure and the correspondingcrank angle to the pressure at IVC in the simulations.764. Engine Simulations1 1.1 1.2 1.3 1.4 1.50510152025Relative Air-Fuel Ratio, ?DryNOxemission,g/kW?hr ExperimentCSEFigure 4.15: Pollutant emission trends with varying ? . Numerical estimateshave been scaled appropriately.bars. Given that a simplified chemistry mechanism is used in a RANS cotnext,the CSE trends show promising agreement with experiments despite considerableunderestimation.4.4.4 DiscussionAccurate calculation of conditionally averaged temperature plays a significant rolein chemical closure using CSE. In a moving-mesh scenario in complex geometries,the statistical homogeneity of conditional averages might be sacrificed near TDCdue to the skewness of finite volume cells. In addition, the application of ? -PDF todescribe the stochastic behaviour of reaction progress variable in premixed flamesis known to be erroneous and is expected to affect the CSE solution.4.5 ConclusionsCSE has been established as a combustion model for simulating industry relevantturbulent reactive flow problems. In view of its first application to simulating com-bustion in a spark-ignition engine, the obtained results are promising and accept-774. Engine Simulationsable in the context of RANS. The Pressure profile, through which heat release andcombustion duration is calculated, is an important quantity of interest along withemissions. The pressure profile is matched to a fair extent, particularly for leanermixtures where more accurate matching is possible for the initial flame kernel evo-lution. CSE has predicted identical trends for NOx emissions despite the use of asimplified chemistry. In conclusion, this work validates the consistency of CSE fortransient turbulent combustion phenomena in an industrial setting. Further workon precise ignition modelling and the implementation of a robust regularisationparameter choice method will be needed to improve the predictive accuracy.785 | Conclusions5.1 Inverse Problem ImplementationA novel methodology has been developed to partition finite volume meshes of ar-bitrarily complex geometries into spatially localized clusters of points.? In the context of parallel processing architectures, spatial locality of parti-tions is chosen over computational load balance in order to maintain the va-lidity of CSE inversion within the partitions treated as ensembles. Althoughthe distribution of mesh elements is not equal, appreciable gains are made ininter-processor communication by preserving locality.? The spatial partitioning method relies on theMorton order space-filling curvegenerated based on a floating-point algorithm which is applicable to unstruc-tured finite volume meshes. In conjunction with the mapping algorithm, thedeveloped algorithm automates the process of localized ensemble construc-tion by finding locality destructive leaps.? While cluster locality is improved against equi-sized partitioning of the Mor-ton order curve, computational load balance is sacrificed to some degreedespite algorithms in place to maintain equal distribution of finite volumesacross partitions. The degree of computational load balance achieved is ac-ceptable for large-scale parallel processing environments and is expected toimprove with mesh refinement.Therefore, a robust technique has been developed that can be applied for inexpen-sive and unsupervised partitioning of industrial geometries.795. Conclusions5.2 Industrial Modelling ApplicationCSE has been established as a combustion model for simulating industry relevantturbulent premixed combustion problems.? As a first application of CSE to combustion simulations in spark-ignition en-gines, the pressure trace has been reasonably well predicted. Other quantitiesof interest such as heat release rate are derived from the pressure trace andtheir inaccuracies are attributed to the deviations from experimentally mea-sured pressure traces. In particular, the severe oscillations of the pressureprofile close to TDC are expected to result in extreme errors in the calcula-tion of heat release rates.? The NOx pollutant emissions calculated have shown a good agreement withthe experimentally observed trends for varying mixture fractions. Consider-ing the simplified chemical mechanism employed, emission predictions areof an acceptable standard.? Accurate calculation of conditionally averaged temperature plays a signifi-cant role in chemical closure using CSE. The inversions are affected by thewall; in the case of using a single domain-wide ensemble, the underlyingassumptions are, in fact, invalidated. Moreover, the application of ? PDF todescribe the stochastic behaviour of reaction progress variable in turbulentpremixed flames is known to be erroneous and is expected to affect the CSEsolution.? Ignition modelling is critical to the simulation of combustion in a homogeneous-charge SI engine. While matching the initial kernel growth rate is challeng-ing with the simplified empirical model in place, it nevertheless gives rea-sonable estimates of key variables.This investigation validates the applicability of CSE for transient turbulent com-bustion phenomena in an industrial setting.805. Conclusions5.3 Future Directions? Simulation of turbulent combustion in complex industrial geometries withswirl motion may be investigated as the next step. With the development ofthe spatial partitioning algorithm, an appropriate choice of ensembles can bemade to retain the validity of CSE in industrial burner domains. In particular,specific attention could be given to wall effects such as heat transfer that oth-erwise sacrifice the accuracy of CSE with a single domain-wide ensemble.? The stochastic behaviour of the conditioning variable significantly affects theaccuracy of CSE and, hence, the precision with which transient phenomenasuch as spark-ignition can be modelled. For premixed flames, better descrip-tions of the presumed PDF shape are available. While the ? -PDF capturesthe placement of delta functions, the probability distribution is matched bet-ter by models that account for the flame reactivity and turbulence. The im-plementation of such improved PDFs could lead to sizeable improvementsin the predictive accuracy in SI engines.? The regularisation of CSE solutions done using the Tikhonov method isprone to error due to the L-curve method of optimum regularisation param-eter choice. A variety of parameter choice methods are available and can beexplored to lend robustness to the computational implementation of CSE.? In the context of SI engines, a more detailed physical representation of theplasma charge could replace the employed simplistic model. If the evolutionof the flame kernel is matched better, critical parameters such as maximumpressure and the corresponding crank angle timing can be predicted withgreater accuracy. Moreover, to improve predictions of emissions, a transportequation is necessary for CO, whereas an improved chemical mechanismmust be implemented for NOx chemistry.? The CSE approach could be implemented for non-premixed combustion sim-ulations in diesel engines. Modelling of sprays and their ignition would bethe next challenging leap of this promising chemical closure technique.81References[1] N. Stern. Review on the economics of climate change. London HMTreasury, 2006.[2] J. Tollefson and R. Monastersky. The global energy challenge: awash withcarbon. Nature, 491:654?655, 2012.[3] P. Moin and J. Kim. Tackling turbulence with supercomputers. ScientificAmerican, 1997.[4] S. B. Pope. Small scales, many species and the manifold challenges ofturbulent combustion. Proceedings of the Combustion Institute, 2012.[5] C. L. Fefferman. Existence and smoothness of the navier-stokes equation.The millennium prize problems, pages 57?67, 2000.[6] W. Kendal Bushe and Helfried Steiner. Conditional moment closure forlarge eddy simulation of nonpremixed turbulent reacting flows. Physics ofFluids, 11:1896, 1999.[7] D. Veynante. Turbulent combustion modeling. Progress in Energy andCombustion Science, 28, 2002.[8] J. A. Miller and C. T. Bowman. Mechanism and modeling of nitrogenchemistry in combustion. Progress in Energy and Combustion Science, 15(4):287?338, 1989.[9] D. B. Olson and W. C. Gardiner Jr. An evaluation of methane combustionmechanisms. The Journal of Physical Chemistry, 81(25):2514?2519, 1977.[10] U. Warnatz, J.and Maas and Robert W. Dibble. Combustion: physical andchemical fundamentals, modeling and simulation, experiments, pollutantformation. Springer, 2006.82References[11] C. K. Westbrook and F. L. Dryer. Chemical kinetic modeling ofhydrocarbon combustion. Progress in Energy and Combustion Science, 10(1):1?57, 1984.[12] W. C. Gardiner. Gas-phase Combustion Chemistry. Springer, 2000.[13] F. M. White. Fluid mechanics. 5th. Boston: McGraw-Hill Book Company,2003.[14] S. B. Pope. Turbulent flows. Cambridge University Press, 2000.[15] T. Poinsot and D. Veynante. Theoretical and Numerical Combustion. RTEdwards Inc., 2011.[16] G. K. Batchelor. The theory of homogeneous turbulence. Cambridgeuniversity press, 1982.[17] O. Reynolds. On the dynamical theory of incompressible viscous fluids andthe determination of the criterion. Philosophical Transactions of the RoyalSociety of London. A, 186:123?164, 1895.[18] A. A. Townsend. On the fine-scale structure of turbulence. Proceedings ofthe Royal Society of London. Series A. Mathematical and PhysicalSciences, 208(1095):534?542, 1951.[19] P. Moin and K. Mahesh. Direct numerical simulation: A tool in turbulenceresearch. Annual Review of Fluid Mechanics, 30(1):539?578, 1998.[20] F. R. Menter, M. Kuntz, and R. Langtry. Ten years of industrial experiencewith the sst turbulence model. Turbulence, heat and mass transfer, 4:625?632, 2003.[21] P. R. Spalart. Strategies for turbulence modelling and simulations.International Journal of Heat and Fluid Flow, 21(3):252 ? 263, 2000.[22] R. Bouffanais. Advances and challenges of applied large-eddy simulation.Computers & Fluids, 39(5):735?738, 2010.[23] H. Pitsch. Large-Eddy Simulation of Turbulent Combustion. AnnualReview of Fluid Mechanics, 38(1):453?482, January 2006.[24] T. Poinsot and D. Veynante. Theoretical and numerical combustion. RTEdwards Incorporated, 2005.83References[25] S. B. Pope. Ten questions concerning the large-eddy simulation ofturbulent flows. New Journal of Physics, 6(1):35, 2004.[26] B. E. Launder and D. B. Spalding. The numerical computation of turbulentflows. Computer Methods in Applied Mechanics and Engineering, 3(2):269? 289, 1974.[27] F. J. Weinberg. The first half-million years of combustion research andtoday?s burning problems. Symposium (International) on Combustion, 15(1):1 ? 17, 1975.[28] F. A. Williams. The role of theory in combustion science. Symposium(International) on Combustion, 24(1):1 ? 17, 1992.[29] K. N. C. Bray. The challenge of turbulent combustion. Symposium(International) on Combustion, 26(1):1 ? 26, 1996.[30] N. Peters. Multiscale combustion and turbulence. Proceedings of theCombustion Institute, 32(1):1 ? 25, 2009.[31] R. W. Bilger, S. B. Pope, K. N. C. Bray, and J. F. Driscoll. Paradigms inturbulent combustion research. Proceedings of the Combustion Institute, 30(1):21?42, January 2005.[32] T. Echekki and E. Mastorakos. Turbulent combustion modeling: Advances,new trends and perspectives, volume 95. Springer Science, 2011.[33] R. Borghi. On the structure and morphology of turbulent premixed flames.In C. Casci and C. Bruno, editors, Recent Advances in the AerospaceSciences, pages 117?138. Springer US, 1985.[34] T. Poinsot and D. Veynante. Theoretical and Numerical Combustion. RTEdwards Inc., 2001.[35] D. B. Spalding. Mixing and chemical reaction in steady confined turbulentflames. Symposium (International) on Combustion, 13(1):649 ? 657, 1971.[36] N. Peters. Laminar diffusion flamelet models in non-premixed turbulentcombustion. Progress in Energy and Combustion Science, 10(3):319?339,1984.[37] E. R. Hawkes and R. S. Cant. A flame surface density approach tolarge-eddy simulation of premixed turbulent combustion. Proceedings ofthe Combustion Institute, 28(1):51?58, 2000.84References[38] H. G. Weller, G. Tabor, A. D. Gosman, and C. Fureby. Application of aflame-wrinkling les combustion model to a turbulent mixing layer. InSymposium (International) on Combustion, volume 27, pages 899?907.Elsevier, 1998.[39] W. Kim, S. Menon, and H. C. Mongia. Large-eddy simulation of a gasturbine combustor flow. Combustion Science and Technology, 143(1-6):25?62, 1999.[40] H. Pitsch. A consistent level set formulation for large-eddy simulation ofpremixed turbulent combustion. Combustion and Flame, 143(4):587?598,2005.[41] O. Colin, F. Ducros, D. Veynante, and T. Poinsot. A thickened flame modelfor large eddy simulations of turbulent premixed combustion. Physics ofFluids, 12:1843, 2000.[42] J. A. van Oijen and L. P. H. de Goey. Modelling of premixed laminarflames using flamelet-generated manifolds. Combustion science andtechnology, 161(1):113?137, 2000.[43] M. Ihme and H. Pitsch. Prediction of extinction and reignition innonpremixed turbulent flames using a flamelet/progress variable model:Part 2. application in LES of sandia flames d and e. Combustion andFlame, 155(1?2):90 ? 107, 2008.[44] M. Ihme and H. Pitsch. Prediction of extinction and reignition innon-premixed turbulent flames using a flamelet/progress variable model.part 1: A priori study and presumed PDF closure. Combustion and Flame,155(1?2):70 ? 89, 2008.[45] F. T. C. Yuen and ?. L. G?lder. Turbulent premixed flame front dynamicsand implications for limits of flamelet hypothesis. Proceedings of theCombustion Institute, 34(1):1393 ? 1400, 2013.[46] A. Y. Klimenko. Multicomponent diffusion of various admixtures inturbulent flow. Fluid Dynamics, 25(3):327?334, 1990.[47] R. W. Bilger. Conditional moment closure for turbulent reacting flow.Physics of Fluids A: Fluid Dynamics, 1993.[48] G. De Paola, E. Mastorakos, Y. M. Wright, and K. Boulouchos. Dieselengine simulations with multi-dimensional conditional moment closure.Combustion Science and Technology, 180(5):883?899, 2008.85References[49] S. Amzin, N. Swaminathan, J. W. Rogerson, and J. H. Kent. Conditionalmoment closure for turbulent premixed flames. Combustion Science andTechnology, 184(10-11):1743?1767, 2012.[50] B. Thornber, R. W. Bilger, A. R. Masri, and E. R. Hawkes. An algorithmfor LES of premixed compressible flows using the conditional momentclosure model. Journal of Computational Physics, 230(20):7687 ? 7705,2011.[51] A. R. Kerstein. One-dimensional turbulence: model formulation andapplication to homogeneous turbulence, shear flows, and buoyant stratifiedflows. Journal of Fluid Mechanics, 392(1):277?334, 1999.[52] A. R. Kerstein. A linear- eddy model of turbulent scalar transport andmixing. Combustion Science and Technology, 60(4?6):391?421, 1988.[53] S. B. Pope. PDF methods for turbulent reactive flows. Progress in Energyand Combustion Science, 11(2):119 ? 192, 1985.[54] A. Y. Klimenko and S. B. Pope. The modeling of turbulent reactive flowsbased on multiple mapping conditioning. Physics of Fluids, 15:1907, 2003.[55] M. J. Cleary and A. Y. Klimenko. A generalised multiple mappingconditioning approach for turbulent combustion. Flow, turbulence andcombustion, 82(4):477?491, 2009.[56] H. Steiner and W. K. Bushe. Large eddy simulation of a turbulent reactingjet with conditional source-term estimation. Physics of Fluids, 13(3):754,2001.[57] P. C. Hansen. Numerical tools for analysis and solution of Fredholmintegral equations of the first kind. Inverse problems, 849, 1992.[58] M. M. Salehi, W. K. Bushe, and K. J. Daun. Application of the conditionalsource-term estimation model for turbulence-chemistry interactions in apremixed flame. Combustion Theory and Modelling, 16(2):301?320, April2012.[59] M. M. Salehi. Numerical Simulation of Turbulent Premixed Flames withConditional Source-Term Estimation. PhD thesis, 2012.[60] W. K. Bushe and H. Steiner. Laminar flamelet decomposition forconditional source-term estimation. Physics of Fluids, 15(6):1564, 2003.86References[61] B. Jin, R. Grout, and W. K. Bushe. Conditional source-term estimation as amethod for chemical closure in premixed turbulent reacting flow. Flow,turbulence and combustion, 81(4):563?582, 2008.[62] M. Wang, J. Huang, and W. K. Bushe. Simulation of a turbulentnon-premixed flame using conditional source-term estimation withtrajectory generated low-dimensional manifold. Proceedings of theCombustion Institute, 31(2):1701?1709, 2007.[63] A. K. Jain, M. N. Murty, and P. J. Flynn. Data clustering: a review. ACMcomputing surveys (CSUR), 31(3):264?323, 1999.[64] A. K. Jain. Data clustering: 50 years beyond k-means. Pattern RecognitionLetters, 31(8):651?666, 2010.[65] G. F. Jenks. The data model concept in statistical mapping. In KonradFrenzel, editor, International Yearbook of Cartography, volume 7, pages186+, USA, 1967. ICA, Rand McNally & Co.[66] J. A Hartigan. Statistical theory in clustering. Journal of Classification, 2(1):63?76, 1985.[67] R. N. Dav? and R. Krishnapuram. Robust clustering methods: a unifiedview. Fuzzy Systems, IEEE Transactions on, 5(2):270?293, 1997.[68] P. Indyk and R. Motwani. Approximate nearest neighbors: towardsremoving the curse of dimensionality. In Proceedings of the thirtiethannual ACM symposium on Theory of computing, STOC ?98, pages604?613, New York, NY, USA, 1998. ACM.[69] H. D. Simon. Partitioning of unstructured problems for parallel processing.Computing Systems in Engineering, 2(2?3):135 ? 148, 1991.[70] B. Hendrickson and K. Devine. Dynamic load balancing in computationalmechanics. Computer Methods in Applied Mechanics and Engineering,184(2):485?500, 2000.[71] B. Smith, P. Bjorstad, and W. Gropp. Domain decomposition. CambridgeUniversity Press, 2004.[72] M. Wittmann, T. Zeiser, G. Hager, and G. Wellein. Domain decompositionand locality optimization for large-scale lattice boltzmann simulations.Computers and Fluids, 80(0):283 ? 289, 2013.87References[73] M. J. Berger and S. H. Bokhari. A partitioning strategy for nonuniformproblems on multiprocessors. Computers, IEEE Transactions on, C?36(5):570?580, 1987.[74] M. T. Jones and P. E. Plassmann. Computational results for parallelunstructured mesh computations. Computing Systems in Engineering, 5(4?6):297 ? 309, 1994.[75] J. R Gilbert, G. L Miller, and S. Teng. Geometric mesh partitioning:Implementation and experiments. SIAM Journal on Scientific Computing,19(6):2091?2110, 1998.[76] T. Asano, D. Ranjan, T. Roos, E. Welzl, and P. Widmayer. Space fillingcurves and their use in the design of geometric data structures. InLATIN?95: Theoretical Informatics, pages 36?48. Springer, 1995.[77] M. J. Aftosmis and M. J. Berger. Applications of space-filling curves toCartesian methods for CFD, 2004.[78] K. Schloegel, G. Karypis, and V. Kumar. Sourcebook of parallelcomputing. chapter Graph partitioning for high-performance scientificsimulations, pages 491?541. Morgan Kaufmann Publishers Inc., SanFrancisco, CA, USA, 2003.[79] T. N. Bui and C. Jones. Finding good approximate vertex and edgepartitions is np-hard. Information Processing Letters, 42(3):153 ? 159,1992.[80] C. Farhat and M. Lesoinne. Automatic partitioning of unstructured meshesfor the parallel solution of problems in computational mechanics.International Journal for Numerical Methods in Engineering, 36(5):745?764, 1993.[81] J. D. Teresco, K. D. Devine, and J. E. Flaherty. Partitioning and dynamicload balancing for the numerical solution of partial differential equations.In A. Bruaset and A. Tveito, editors, Numerical Solution of PartialDifferential Equations on Parallel Computers, volume 51 of Lecture Notesin Computational Science and Engineering, pages 55?88. Springer BerlinHeidelberg, 2006.[82] G. Zumbusch. On the quality of space-filling curve induced partitions. Z.Angew. Math. Mech, 81:25?28, 2000.88References[83] M. Bader. Space-Filling Curves: An Introduction with Applications inScientific Computing, volume 9. Springer, 2012.[84] B. Moon, H. V. Jagadish, C. Faloutsos, and J. H. Saltz. Analysis of theclustering properties of the hilbert space-filling curve. Knowledge andData Engineering, IEEE Transactions on, 13(1):124?141, 2001.[85] M. Parashar and J. C. Browne. On partitioning dynamic adaptive gridhierarchies. In System Sciences, 1996., Proceedings of the Twenty-NinthHawaii International Conference on,, volume 1, pages 604?613, 1996.[86] J. R. Pilkington and S. B. Baden. Dynamic partitioning of non-uniformstructured workloads with spacefilling curves. Parallel and DistributedSystems, IEEE Transactions on, 7(3):288?300, 1996.[87] C. Gotsman and M. Lindenbaum. On the metric properties of discretespace-filling curves. Image Processing, IEEE Transactions on, 5(5):794?797, 1996.[88] X. Liu and G. F. Schrack. An algorithm for encoding and decoding the 3-DHilbert order. IEEE transactions on image processing: a publication of theIEEE Signal Processing Society, 6(9):1333?7, 1997.[89] X. Liu. Four alternative patterns of the hilbert curve. Applied Mathematicsand Computation, 147(3):741?752, 2004.[90] J. K. Lawder and P. J. H. King. Querying multi-dimensional data indexedusing the hilbert space-filling curve. ACM Sigmod Record, 30(1):19?24,2001.[91] H. Sagan. Space-Filling Curves. Springer-Verlag, New York, 1994.[92] P. M. Campbell, K. D. Devine, J. E. Flaherty, L. G. Gervasio, and J. D.Teresco. Dynamic octree load balancing using space-filling curves.Williams College Department of Computer Science, Tech. Rep. CS-03-01,2003.[93] M. Connor and P. Kumar. Fast construction of k-nearest neighbor graphsfor point clouds. Visualization and Computer Graphics, IEEE Transactionson, 16(4):599?608, 2010.[94] A. Perez, S. Kamata, and E. Kawaguchi. Peano scanning of arbitrary sizeimages. In Pattern Recognition, 1992. Vol.III. Conference C: Image,Speech and Signal Analysis, Proceedings., 11th IAPR InternationalConference on, pages 565?568, 1992.89References[95] S. Tirthapura, S. Seal, and S. Aluru. A formal analysis of space fillingcurves for parallel domain decomposition. International Conference onParallel Processing, pages 505?512, 2006.[96] F. T. C. Yuen and ?. L. G?lder. Dynamics of Lean-Premixed TurbulentCombustion at High Turbulence Intensities. Combustion Science andTechnology, 182(4?6):544?558, 2010.[97] J. Kariuki, J. R. Dawson, and E. Mastorakos. Measurements in turbulentpremixed bluff body flames close to blow-off. Combustion and Flame, 159(8):2589 ? 2607, 2012.[98] C. Reynolds. Performance of a partially stratified-charge natural gasengine. PhD thesis, University of British Columbia, 2001.[99] OpenCFD. Openfoam software, 2013. URL http://www.openfoam.org/.[100] N. Chakraborty and R. S. Cant. Effects of Lewis number on scalar transportin turbulent premixed flames. Physics of Fluids, 21(3):035110, 2009.[101] D. Veynante, J. M. Duclos, and J. Piana. Experimental analysis of flameletmodels for premixed turbulent combustion. In Symposium (International)on Combustion, volume 25, pages 1249?1256. Elsevier, 1994.[102] O. Colin and A. Benkenida. The 3-zones extended coherent flame model(ECFM3Z) for computing premixed/diffusion combustion. Oil & gasscience and technology, 59(6):593?609, 2004.[103] V. Mittal and H. Pitsch. A flamelet model for premixed combustion undervariable pressure conditions. Proceedings of the Combustion Institute, 34(2):2995?3003, 2013.[104] Benjamin T. Zoller, M. L. Hack, and P. Jenny. A PDF combustion modelfor turbulent premixed flames. Proceedings of the Combustion Institute, 34(1):1421?1428, 2013.[105] H. Barths, C. Hasse, G. Bikas, and N. Peters. Simulation of combustion indirect injection diesel engines using a eulerian particle flamelet model.Proceedings of the Combustion Institute, 28(1):1161 ? 1168, 2000.[106] A. Odedra and W. Malalasekera. Eulerian particle flamelet modeling of abluff-body CH4/H2 flame. Combustion and Flame, 151(3):512?531, 2007.90References[107] N. Swaminathan and R. W. Bilger. Analyses of conditional momentclosure for turbulent premixed flames. Combustion Theory and Modelling,5(2):241?260, 2001.[108] J. Huang, P.G. Hill, W.K. Bushe, and S.R. Munshi. Shock-tube study ofmethane ignition under engine-relevant conditions: experiments andmodeling. Combustion and Flame, 136(1?2):25 ? 42, 2004.[109] I. Yimer, I. Campbell, and L.-Y. Jiang. Estimation of the turbulent schmidtnumber from experimental profiles of axial velocity and concentration forhigh-reynolds-number jet flows. Canadian Aeronautics and Space Journal,48(3):195?200, 2002.[110] L. Vervisch, R. Hauguel, P. Domingo, and M. Rullaud. Three facets ofturbulent combustion modelling: DNS of premixed v-flame, LES of liftednon-premixed flame and RANS of jet-flame. Journal of Turbulence, 5:1?36, 2004.[111] M. Thiele, S. Selle, U. Riedel, J. Warnatz, and U. Maas. Numericalsimulation of spark ignition including ionization. Proceedings of thecombustion institute, 28(1):1177?1185, 2000.[112] T. Kravchik and E. Sher. Numerical modeling of spark ignition and flameinitiation in a quiescent methane-air mixture. Combustion and flame, 99(3):635?643, 1994.[113] S. Pischinger and J. B. Heywood. A model for flame kernel development ina spark-ignition engine. In Symposium (International) on Combustion,volume 23, pages 1033?1040. Elsevier, 1991.[114] Z. Tan and R. D. Reitz. An ignition and combustion model based on thelevel-set method for spark ignition engine multidimensional modeling.Combustion and flame, 145(1):1?15, 2006.[115] B. Enaux, V. Granet, O. Vermorel, C. Lacour, C. Pera, C. Angelberger, andT. Poinsot. LES study of cycle-to-cycle variations in a spark ignitionengine. Proceedings of the Combustion Institute, 33(2):3115?3122, 2011.[116] O. Yas?ar. A new ignition model for spark-ignited engine simulations.Parallel computing, 27(1):179?200, 2001.[117] E. Mastorakos. Ignition of turbulent non-premixed flames. Progress inEnergy and Combustion Science, 35(1):57 ? 97, 2009.91References[118] C. K. Law. Heat and mass transfer in combustion: fundamental conceptsand analytical techniques. Progress in energy and combustion science, 10(3):295?318, 1984.[119] D. A. Eichenberger and W. L. Roberts. Effect of unsteady stretch onspark-ignited flame kernel survival. Combustion and flame, 118(3):469?478, 1999.[120] G. Rassweiler and L. Withrow. Motion pictures of engine flames correlatedwith pressure cards. Warrendale, PA : Society of Automotive Engineers,1980.[121] R. K. Hanson and S. Salimian. Survey of rate constants in the N/H/Osystem. In Combustion Chemistry, pages 361?421. Springer US, 1984.92


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items