Turbulent Premixed Combustion Simulation withConditional Source-term Estimation and Linear-EddyModel Formulated PDF and SDR ModelsbyHong P. TsuiB.Sc. (Honours), The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2016c© Hong P. Tsui, 2016AbstractComputational fluid dynamics (CFD) is indispensable in the development of com-plex engines due to its low cost and time requirement compared to experiments.Nevertheless, because of the strong coupling between turbulence and chemistry inpremixed flames, the prediction of chemical reaction rates continues to be a mod-elling challenge.This work focuses on the improvement of turbulent premixed combustion sim-ulation strategies requiring the use of presumed probability density function (PDF)models. The study begins with the development of a new PDF model that includesthe effect of turbulence, achieved by the implementation of the Linear-Eddy Model(LEM). Comparison with experimental burners reveals that the LEM PDF can cap-ture the general PDF shapes for methane-air combustion under atmospheric con-ditions with greater accuracy than other presumed PDF models. The LEM is ad-ditionally used to formulate a new, pseudo-turbulent scalar dissipation rate (SDR)model.Conditional Source-term Estimation (CSE) is implemented in the Large EddySimulation (LES) of the Gu¨lder burner as the closure model for the chemistry-turbulence interactions. To accommodate the increasingly parallel computationalenvironments in clusters, the CSE combustion module has been parallelised andoptimised. The CSE ensembles can now dynamically adapt to the changing flamedistributions by shifting their spatial boundaries and are no longer confined topre-allocated regions in the simulation domain. Further, the inversion calculationis now computed in parallel using a modified version of an established iterativesolver, the Least-Square QR-factorisation (LSQR). The revised version of CSEdemonstrates a significant reduction in computational requirement — a reductioniiAbstractof approximately 50% — while producing similar solutions as previous implemen-tations.The LEM formulated PDF and SDR models are subsequently implemented inconjunction with the optimised version of CSE for the LES of a premixed methane-air flame operating in the thin reaction zone. Comparison with experimental mea-surements of temperature reveals that the LES results are very comparable in termsof the flame height and distribution. This outcome is encouraging as it appearsthat this work represents a significant step towards the correct direction in devel-oping a complete combustion simulation strategy that can accurately predict flamecharacteristics in the absence of ad hoc parameters.iiiPrefaceThe content of Chapter 3 is published in:• H. P. Tsui and W. K. Bushe. Linear-eddy model formulated probability den-sity function and scalar dissipation rate models for premixed combustion.Flow, Turbulence and Combustion, 93:487–503, 2014.I was responsible for conducting all parts of the research and preparing the finalmanuscript in this publication. Dr. W. Kendal Bushe supervised the research, re-viewed the manuscript and provided feedback as appropriate.The content of Chapter 4 is published in:• H. P. Tsui, M. M. Kamal, S. Hochgreb and W. K. Bushe. Direct comparisonof PDF and scalar dissipation rates between LEM simulations and experi-ments for turbulent, premixed methane air flames. Combustion and Flame,165:208–222, 2016.I was responsible for conducting all parts of the numerical research and preparingthe final manuscript in this publication. Dr. Simone Hochgreb and Mr. M. M.Kamal of University of Cambridge provided the relevant experimental data andexpertise for the analysis of the measured results. Dr. W. Kendal Bushe supervisedthe research and provided feedback as appropriate. Dr. Hochgreb and Dr. Bushereviewed the manuscript and provided feedback as appropriate.The content of Chapter 5 is published in:ivPreface• H. P. Tsui and W. K. Bushe. Conditional source-term estimation using dy-namic ensemble selection and parallel iterative solution. Combustion, The-ory and Modelling, 2016. DOI: 10.1080/13647830.2016.1178811I was responsible for conducting all parts of the research and preparing the finalmanuscript in this publication. Dr. W. Kendal Bushe supervised the research, re-viewed the manuscript and provided feedback as appropriate.The content of Chapter 6 has been accepted for presentation in the CombustionSymposium (2016):• H. P. Tsui, M. M. Salehi and W. K. Bushe. LES of a turbulent premixed flamewith conditional source-term estimation and linear-eddy model presumedPDF. 36th Combustion Symposium, Seoul, Korea, 2016.I was responsible for conducting all parts of the research and preparing the finalmanuscript in this publication. Dr. M. M. Salehi provided feedback for the simu-lation methodologies as appropriate and reviewed the manuscript. Dr. W. KendalBushe supervised the research, reviewed the manuscript and provided feedback asappropriate.vTable of ContentsTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Literature Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Conservation Equations . . . . . . . . . . . . . . . . . . . . . . . 52.3 Transport Properties . . . . . . . . . . . . . . . . . . . . . . . . . 82.4 Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 102.5 Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.5.1 Chemical Kinetics . . . . . . . . . . . . . . . . . . . . . 122.5.2 Reduced Chemistry . . . . . . . . . . . . . . . . . . . . . 13viTable of Contents2.6 Laminar Premixed Flame . . . . . . . . . . . . . . . . . . . . . . 142.7 Scales of Turbulence . . . . . . . . . . . . . . . . . . . . . . . . 182.8 Regimes of Premixed Turbulent Combustion . . . . . . . . . . . . 222.9 Computational Fluid Dynamics (CFD) . . . . . . . . . . . . . . . 272.10 Turbulent Premixed Combustion Modelling . . . . . . . . . . . . 343 Pseudo-Turbulent Probability Density Function and Scalar Dissipa-tion Rate Models for Premixed Combustion . . . . . . . . . . . . . . 443.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.2 Probability Density Function Models . . . . . . . . . . . . . . . . 463.3 Linear-Eddy Model . . . . . . . . . . . . . . . . . . . . . . . . . 493.3.1 Evolutionary Equations . . . . . . . . . . . . . . . . . . . 493.3.2 Stochastic Advection . . . . . . . . . . . . . . . . . . . . 503.4 LEM Simulation Parameters . . . . . . . . . . . . . . . . . . . . 513.4.1 LEM Methodology for Experimental Study . . . . . . . . 513.4.2 LEM Methodology for DNS Study . . . . . . . . . . . . 533.5 Construction of PDF Models . . . . . . . . . . . . . . . . . . . . 533.5.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . . 563.6 Scalar Dissipation Rate . . . . . . . . . . . . . . . . . . . . . . . 573.7 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 593.7.1 Probability Density Functions . . . . . . . . . . . . . . . 593.7.2 Scalar Dissipation Rate . . . . . . . . . . . . . . . . . . . 633.8 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . 644 Direct Comparison of PDF and SDR between LEM Simulations andExperiments for Turbulent, Premixed Methane-Air Flames . . . . . 654.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.2 Numerical Conditions: Premixed Combustion . . . . . . . . . . . 694.2.1 Linear-Eddy Model . . . . . . . . . . . . . . . . . . . . . 694.2.2 LEM Simulation Methods . . . . . . . . . . . . . . . . . 714.3 Experimental Conditions: Stratified Swirl Burner . . . . . . . . . 754.4 Construction of PDF Models . . . . . . . . . . . . . . . . . . . . 774.4.1 LEM Flame Profiles . . . . . . . . . . . . . . . . . . . . 77viiTable of Contents4.4.2 LEM PDF Construction . . . . . . . . . . . . . . . . . . 794.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.5.1 Probability Density Function . . . . . . . . . . . . . . . . 824.5.2 Scalar Dissipation Rate . . . . . . . . . . . . . . . . . . . 864.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.6.1 Probability Density Function . . . . . . . . . . . . . . . . 914.6.2 Scalar Dissipation Rate . . . . . . . . . . . . . . . . . . . 934.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955 Conditional Source-term Estimation: Parallel Iterative Solution withDynamic Ensembles . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.2 Conditional Source-term Estimation . . . . . . . . . . . . . . . . 995.3 Optimisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.3.1 Dynamic Ensemble Selection . . . . . . . . . . . . . . . 1025.3.2 Matrix Inversion - Iterative Solver . . . . . . . . . . . . . 1105.4 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . 1135.4.1 Standalone Inversion Tests . . . . . . . . . . . . . . . . . 1135.4.2 LES of Axisymmetric Burner . . . . . . . . . . . . . . . 1145.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 1165.5.1 Standalone Inversion Tests . . . . . . . . . . . . . . . . . 1165.5.2 LES of Axisymmetric Burner . . . . . . . . . . . . . . . 1215.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1256 LES of Turbulent Premixed Flames with Optimised CSE and LEMSubmodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1266.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1266.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1286.2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . 1286.2.2 Chemistry Reduction . . . . . . . . . . . . . . . . . . . . 1286.2.3 Presumed PDF and SDR models . . . . . . . . . . . . . . 1296.3 Numerical Methods / Simulation Setup . . . . . . . . . . . . . . . 1326.3.1 Pre-Process . . . . . . . . . . . . . . . . . . . . . . . . . 132viiiTable of Contents6.3.2 Non-Reactive Simulation . . . . . . . . . . . . . . . . . . 1326.3.3 Reactive Simulation . . . . . . . . . . . . . . . . . . . . 1346.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 1356.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . 1387 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . 1397.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1397.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147ixList of TablesList of TablesTable 4.1 Relevant LEM simulation parameters: l0 and ηK are the integraland Kolmogorov scales. Other constant parameters are invari-ant between the cases, including δ f , sL, Cλ and Nη , which arerespectively, 588 µm, 0.214 m/s, 15.0 and 1.0. A minimum of1,000 computational cells are used for each simulation. . . . . 74Table 4.2 Operating Conditions: (1) SFR = ratio of split flow to swirlers tototal flow, (2) SN = measured swirl number, ratio of tangentialto axial momentum, (3) Maximum total u′/sL at the midpointof the flame brush at z = 30 mm [71, 194]. . . . . . . . . . . . 75Table 5.1 Relevant LES simulation parameters. The values of a and bare computed according to the approach of Billson [15]: a =exp(−∆t/τt) and b =√1−a2, where ∆t and τt are respectivelythe timestep size and the integral time scale. . . . . . . . . . . 116xList of TablesTable 5.2 Comparison of inversion times between LU-decomposition andLSQR at a number of different settings. The |e|, ~0 and x0 re-spectively denote the error tolerances in the solver, a zero ini-tial solution and a non-zero initial solution. The associated L2-norms are printed below the convergence times. Typical ma-trix sizes for one-condition CSE are represented by the first twocolumns while typical matrix sizes for two-condition CSE arerepresented by the last three columns. Failed cases are denotedby ‘×’, where lower error tolerances are required to providestable solutions. From the smallest to the largest matrix, 1,000,500, 100, 20 and 2 samples are used to obtain the average inver-sion times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118Table 5.3 Comparison of inversion times for single ensembles of varyingsizes with the MPI-LSQR implementation using a number ofdifferent processors. The numbers x × y in the first columnrepresent the number of nodes (x) and the number of processorsper node (y) used in the computation. The error tolerances arekept constant at 10−10 between all of the cases for consistency. 120Table 5.4 Summary of simulation and execution times. The executiontimes for LSQR with Semi-Dynamic Ensemble Selection aretypically 50% of the LU-decomposition method for larger sys-tems. The executions times for MPI-LSQR with Dynamic En-semble Selection are further decreased by 2–22% in comparisonwith the LSQR procedure. The ensemble sizes were decreasedfor the coarse grid (52,500) case in order to maintain localityfor cell information, leading to more comparable execution times.125Table 6.1 Summary of experimental conditions. Symbols: UB = bulk flowvelocity in the central jet; u′ = total turbulence intensity; ΛL= longitudinal integral length scale; Ret = turbulent Reynoldsnumber; Ka = turbulent Karlovitz number; Da = turbulent Damko¨hlernumber. The values are retrieved from [173]. . . . . . . . . . . 133xiList of TablesTable 6.2 Relevant synthetic turbulence parameters. Symbols: hg = gridlocation; K = turbulent kinetic energy; Λ = integral length; aand b = temporal correlation parameters; ∆t = timestep size. . . 134Table 6.3 Summary of simulation times. . . . . . . . . . . . . . . . . . . 138xiiList of FiguresList of FiguresFigure 2.1 One-dimensional laminar flame structure for premixed, stoi-chiometric, methane-air combustion at standard conditions. . . 15Figure 2.2 Laminar flame speeds and equilibrium temperatures of pre-mixed methane-air flames at different equivalence ratios in stan-dard ambient temperature (298.15 K) and pressure (100 kPa)conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 17Figure 2.3 A schematic diagram of the turbulent kinetic energy spectrumat high Reynolds numbers. . . . . . . . . . . . . . . . . . . . 21Figure 2.4 Effect of large-scale turbulent eddies on a premixed flame front. 23Figure 2.5 Modified Borghi diagram illustrating different turbulent pre-mixed combustion regimes. . . . . . . . . . . . . . . . . . . . 24Figure 2.6 A schematic representation of the resolved and modelled en-ergy spectrum in RANS. . . . . . . . . . . . . . . . . . . . . 29Figure 2.7 Schematic view of the resolved and subgrid-scale structures ina grid with spacing of ∆: Resolved structures (light gray) andsubgrid-scale structures (dark gray). . . . . . . . . . . . . . . 31Figure 2.8 A symbolic representation of the resolved and modelled turbu-lent kinetic energy spectrum in LES. . . . . . . . . . . . . . . 32Figure 2.9 CSE operational flowchart. . . . . . . . . . . . . . . . . . . . 43Figure 3.1 PDFs are exemplified at two axial locations of the DNS sim-ulation domain (Ld = 6 mm). Solid: DNS; dash: modifiedlaminar flamelet; dash dot: laminar flamelet; dot: β -PDFs [70]. 49xiiiList of FiguresFigure 3.2 Triplet map applied to a line. Cell values are transported withinthe eddy interval; interpolation is not required at any point [77].The before (left) and after (right) domain values are illustrated. 51Figure 3.3 Premixed laminar flame solution at an equivalence ratio of 0.73for the six-step reduced mechanism (line) and Cantera [58](symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Figure 3.4 (a) (x1,c1) and (x2,c2) mark the truncation limits where onlythe cells within the interval are retained; (b) Modified laminarflamelet PDF with c¯= 0.5 and c′2n = 0.33 is constructed usingthe truncated flame profile. . . . . . . . . . . . . . . . . . . . 54Figure 3.5 Individual flame realisations from LEM calculations. . . . . . 55Figure 3.6 (a) The portions to be retained are within the intervals markedby ’o’. The truncation positions are selected such that the re-sultant PDF would have c¯ = 0.5 and c′2n = 0.33. The trunca-tion boundaries are different for each realisation because theprofiles change with time. (b) The modified laminar flamelet(dash-circle) and LEM (solid) PDFs of similar mean and vari-ance are shown. . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 3.7 LEM generated PDFs at Ret = 50. Each line represents a dif-ferent number of LEM profiles used for the construction of thePDF. Solid: 5000 profiles; dash: 2000 profiles; dash-dot: 1000profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.8 Experimental PDFs conditioned at three different means of thereaction progress variable. Solid: experiment [167]; dot: mod-ified laminar flamelet; dash: LEM. The vertical axis representsthe probability at state c∗ conditioned by c¯, P(c)|c¯. . . . . . . 61Figure 3.9 PDF at x/Ld = 0.391 of the DNS simulation domain (Ld =6 mm). Solid: DNS; dash: modified laminar flamelet; dashdot: LEM. The DNS and modified laminar flamelet PDFs areobtained from [70]. . . . . . . . . . . . . . . . . . . . . . . . 62Figure 3.10 Scalar dissipation conditioned in c-space. Solid: experiment [167];dot: laminar flame; dash: LEM at Ret = 38; dash-dot: LEM atRet = 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . 63xivList of FiguresFigure 4.1 Premixed laminar flame solution at an equivalence ratio of0.73 for the six-step global mechanism (line) and Cantera [58](symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 4.2 Borghi diagram showing locations of the ten prototype LEMflames (1a to 4b). The LEM test cases are represented by ’◦’.The experimental flames are represented by triangles: SwB1(’C’), SwB2 (’5’), and SwB3 (’B’). It has been found thatincreasing the integral length beyond an order of magnitudeabove the laminar flame thickness while holding the turbulentfluctuations constant does not significantly alter the PDF pro-files for the LEM simulations. . . . . . . . . . . . . . . . . . 73Figure 4.3 Characteristic LEM temperature profiles of two prototype flamesat different Ret . The individual profiles on each graph are sep-arated by at least one large eddy turnover time. . . . . . . . . 78Figure 4.4 (a) (x1,c1) and (x2,c2) mark the truncation limits; only the cellswithin the interval are retained. (b) The modified laminar flameletPDF with c¯ = 0.50 and c′2n = 0.39 constructed using the trun-cated flame profile is shown. . . . . . . . . . . . . . . . . . . 80Figure 4.5 (a) The portions to be retained are within the intervals delim-ited by the circles. The four truncation positions are selectedsuch that the resultant PDF would have c¯ = 0.5 and c′2n = 0.39.The truncation boundaries are different for each temperatureprofile because of the transient effects. (b) The modified lami-nar flamelet (dash) and LEM (solid) PDFs of similar mean andvariance are shown. . . . . . . . . . . . . . . . . . . . . . . . 81Figure 4.6 PDF models at various c¯ and c′2n are shown; each row repre-sents one value of c¯ and three values of c′2n (ranging from leftto right: 0.25, 0.49 and 0.75). Vertical and horizontal axis oneach graph represent the probability and the progress variable,respectively. Solid: case 4b, dash: case 3b, dash dot: case 2b,dot: case 1b (notation is in accordance with Figure 4.2). . . . 83xvList of FiguresFigure 4.7 PDFs measured from the SwB1 flame at various axial locationsconditioned by different distribution means, c¯ (solid). The cor-responding LEM PDFs of similar c¯ and c′2n are superimposed(dash). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 4.8 PDFs measured from the SwB2 flame at various axial locationsconditioned by different distribution means, c¯ (solid). The cor-responding LEM PDFs of similar c¯ and c′2n are superimposed(dash). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85Figure 4.9 PDFs measured from the SwB3 flame at various axial loca-tions, conditioned by different distribution means, c¯ (solid).The corresponding LEM PDFs of similar c¯ and c′2n are super-imposed (dash). . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 4.10 Non-dimensionalized conditional average of the SDR modelsfrom the ten prototype flames. The left and right columns illus-trate changes in the turbulent fluctuations and integral lengthscales, respectively. The turbulent Reynolds number of eachcase is recorded in the parentheses. The laminar case is dis-played in light gray as a reference. . . . . . . . . . . . . . . . 88Figure 4.11 Non-dimensionalized conditional average of the SDR from thethree experimental flames. Each column represents a uniqueaxial position from the flame stabilisation point while each rowrepresents one swirling condition. The error bars indicate +/-one standard deviation from the mean. The laminar case isdisplayed in light gray as a reference. . . . . . . . . . . . . . 89Figure 4.12 Functional dependence of the scalar dissipation ( f (c∗)) at threeaxial positions downstream from the flame stabilisation point(10, 30 and 60 mm) with the LEM results superimposed. Theexperimental flames are represented by symbols: SwB1 (’+’),SwB2 (’’), and SwB3 (’◦’). The LEM result is representedby the solid line. The laminar case is displayed in light gray asa reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . 90xviList of FiguresFigure 4.13 Values of χ0 in relation to the distance downstream of the flamestabilisation point. The experimental flames are representedby symbols: SwB1 (’+’), SwB2 (’’), and SwB3 (’◦’). TheLEM (black line) and laminar flame (gray line) results are notdependent on the axial position. . . . . . . . . . . . . . . . . 90Figure 4.14 Non-dimensionalized unconditional mean scalar dissipation rates(χ¯c×δ f /sL) as predicted by (a) LEM and (b) experiments withLEM solutions superimposed at various combinations of c¯ andc′2n. Solid: LEM; dash: SwB1; dash-dot: SwB2; dot: SwB3. . 92Figure 5.1 CSE operational flowchart. . . . . . . . . . . . . . . . . . . . 101Figure 5.2 Computational domain incorporating different number of staticensembles. The volume of each ensemble decreases with in-creasing number of processors. . . . . . . . . . . . . . . . . . 104Figure 5.3 The three possible scenarios that one would encounter duringthe construction of dynamic ensembles are illustrated: N <M,N =M and N >M, where N and M are respectively the numberof processors available and the number of ensembles required. 105Figure 5.4 Semi-Dynamic Ensemble Selection processor map. . . . . . . 107Figure 5.5 Dynamic Ensemble Selection schematic. . . . . . . . . . . . 109Figure 5.6 MPI-LSQR iteration. Each n participating processor performs1/n of the total matrix-vector multiplications. . . . . . . . . . 113Figure 5.7 Comparison of inversion times between LU-decomposition andLSQR at a number of different error tolerance settings (|e|).This is a graphical representation of the data in Table 5.2. Eachof the five plots represents a different matrix size. For eachplot, the bars represent inversion times of (from left to right):LU, LSQR (~0) (|e| = 10−10), LSQR (x0) (|e| = 10−10), LSQR(x0) (|e| = 10−8), LSQR (x0) (|e| = 10−6), and LSQR (x0) (|e|= 10−4). The ‘×’ represents the failed cases. . . . . . . . . . 117xviiList of FiguresFigure 5.8 Comparison of inversion times with the MPI-LSQR imple-mentation using a number of different processors counts. Thisis a graphical representation of the data in Table 5.3. Each ofthe five plots represents a different matrix size. For each plot,the bars represent inversion times of (from left to right): LU,LSQR (~0), LSQR (x0), (1 × 2), (2 × 1), (1 × 3), (3 × 1), (1× 4), (4 × 1), (1 × 8), (8 × 1), (1 × 12), (12 × 1), (1 × 24),and (1 × 48). The numbers (x × y) represent the number ofnodes (x) and the number of processors per node (y) used in theMPI-LSQR computation. The inversion for the matrix of size100,000 × 2,500 is not performed by the LU-decompositiondue to an exceedingly exhaustive computational requirement;it is marked by ‘×’. . . . . . . . . . . . . . . . . . . . . . . . 119Figure 5.9 The radial temperature and species mass fraction profiles at sixaxial locations downstream of the inlet for the three grids areshown. The profiles at z/D = 1 are shown in black, then eachsuccessive jet diameter downstream is shown in a lighter shadeof gray, up to z/D = 6. The jet diameter, D, is 1.12·10−2 m.The symbols represent: LU–manual: ‘—’; LSQR–semi-dynamic:‘4’; PLSQR–dynamic: ‘•’. . . . . . . . . . . . . . . . . . . 123Figure 6.1 Burner cross-section. . . . . . . . . . . . . . . . . . . . . . . 133Figure 6.2 Average binarised temperature-based progress variable profileand characteristic flame height from the LES result (left) andexperimental measurements (right), where D = 11.2 mm. Thesolid line represents the half-burning surface, at which the av-erage progress variable takes on a value of 0.5. . . . . . . . . 136xviiiList of SymbolsList of SymbolsRoman SymbolsA PDF matrixA Pre-exponential factor in Arrhenius equationAc Cross-sectional area of streamtubeAL Laminar flame surface areaAT Turbulent flame surface areaa Coefficient for Billson’s temporal correlation (1)ai,α Viscosity fitting coefficientsb Mass fraction arrayb Coefficient for Billson’s temporal correlation (2)bi,α Thermal conductivity fitting coefficientsCBML Bray Moss Libby parameterCEBU Eddy Break-Up constantCEDM Eddy Dissipation Model constantCκ Universal constant for the energy spectrum in the inertial subrangeCλ LEM empirical parameter (1)c Progress variablec∗ Discretised progress variablec¯ Mean of the progress variablexixList of Symbolsc˜ Favre-averaged mean of the progress variablec′2 Variance of the progress variablec˜′2 Favre-averaged variance of the progress variablec′2n Normalised variance of the progress variableci NASA polynomial coefficientscp Mixture-averaged heat capacitycpα Heat capacity of species αD Molecular diffusivityDa Damko¨hler numberDα Molecular diffusivity of species αDαβ Binary diffusion coefficientdi,αβ Binary diffusion fitting coefficientsE Total EnergyEa Activation energyeˆi Unit vector in the ith direction~g Gravity or body force in generalh Mixture-averaged enthalpyhα Enthalpy of species α∆h0fα Enthalpy of formation of species αhg Grid locationhsα Sensible enthalpy of species αJα Mass diffusive flux of species αK Maximum number of processors in ensembleKa First Karlovitz numberKaδ Second Karlovitz numberk Turbulent kinetic energyxxList of SymbolskArr Arrhenious reaction rate constantL Characteristic flow length scalel0 Integral length scalelL Longitudinal integral length scalelδ Laminar flame thickness of inner layerM Number of ensemblesM Number of realisationsM˙ Mass flow rateM Mixture-averaged molar massMα Molar mass of species αmα Mass of species α in mixtureN Number of processorsN Number of speciesNη LEM empirical parameter (2)nα Mole number of species α in mixturep PressureQ˙ Heat source term~q Total energy fluxR Mixture-averaged gas constantRi j Two-point correlation, autocovarianceRu Universal gas constantRα Gas constant of species αRe Reynolds numberRet Turbulent Reynolds numberSct Turbulent Schmidt numbersd Displacement speed of iso-surface of the progress variablexxiList of SymbolssL Laminar flame speedsT Turbulent flame speedsα Entropy of species αT TemperatureTb Burnt gas temperatureTeq Temperature at equilibriumTre f Reference temperature for enthalpy of formationTu Unburnt gas temperaturet Time coordinate∆t Timestep sizeU Characteristic flow velocityUB Bulk flow velocity~u Velocityu Mean velocityu′ Turbulent intensityu′′ Velocity fluctuation around the Favre averageuK Kolmogorov velocityVα Diffusion velocity of species αXα Mole fraction of species αx Conditional mass fraction array~x Spatial coordinateYα Mass fraction of species αY ′2P Variance of the product mass fractionz Axial coordinateGreek SymbolsαBML Bray Moss Libby PDF parameter (1)xxiiList of Symbolsαc Thermal diffusivityβBML Bray Moss Libby PDF parameter (2)∆ Grid spacingδ f Laminar flame thicknessδi j Kronecker deltaε Dissipation of turbulent kinetic energyηK Kolmogorov lengthκ Wave numberΛ Thermal conductivityΛα Thermal conductivity of species αλ Oxidizer-to-fuel equivalence ratioλd Damping coefficient in Tikhonov regularisationλLEM LEM eddy event frequency per unit length of the domainµ Dynamic viscosityµα Dynamic viscosity of species αν Kinematic viscosityνT Coefficient of turbulent kinematic viscosityΞ Flame wrinkling factorΞ∆ SGS flame wrinkling factorρ Densityρu Unburnt mixture densityΣ SGS flame surface densityσi j Stress tensorσ sgsi j Subgrid-scale turbulent tress tensorτBML BML turbulent mixing time scaleτc Chemical time scalexxiiiList of Symbolsτi j Viscous stress tensorτK Kolmogorov time scaleτt Turbulent integral time scaleφ Fuel-to-oxidizer equivalence ratioχc Scalar dissipation rateω˙α Chemical source term of species αSuperscriptsB Temperature exponent in Arrhenius equationT Transpose of matrixSubscriptsα , β Species indexxxivList of AbbreviationsList of AbbreviationsBML Bray Moss LibbyCFD Computational Fluid DynamicsCMC Conditional Moment ClosureCSC Coherent Structure CaptureCSE Conditional Source-term EstimationDNS Direct Numerical SimulationEBU Eddy Break-UpEDM Eddy Dissipation ModelFDF Filtered Density FunctionFGM Flamelet-Generated ManifoldFPI Flame Prolongation of ILDMFSD Flame Surface DensityLEM Linear-Eddy ModelLES Large Eddy SimulationPCM Presumed Conditional MomentxxvList of AbbreviationsPDF Probability Density FunctionRANS Reynolds-Averaged Navier-StokesSDR Scalar Dissipation RateSDS Semi-Deterministic SimulationSFS Subfilter-ScaleSGS Subgrid-ScaleUHC Unburnt HydrocarbonsURANS Unsteady Reynolds-Averaged Navier-StokesVLES Very Large Eddy SimulationxxviAcknowledgmentsI would like to take this opportunity to express my appreciation to those who havehelped and supported my journey in the pursuit of knowledge in the field of turbu-lent combustion over the past five years. First and foremost, I would like extend mydeepest gratitude to my supervisor, Professor W. Kendal Bushe, for his patience,trust and guidance over the past few years. Without his knowledge and advice onthe myriad of problems, this research would not have been possible.I would also like to thank Professors Carl Ollivier-Gooch, Dana Grecov andJames Feng for dedicating their time to attend my proposal defense and to providevaluable feedback throughout different stages of this work. I am grateful to myfriends and colleagues in the combustion research group for their positive companyand valuable discussions; in particular, Dr. M. Mahdi Salehi invariably manages tospare time and share his insights and passion in the field of combustion research. Ialso thank Graham Hendra for the engaging discussions and thorough feedbacks,Girish Nivarti and Stefanie Asenbauer for enriching my graduate experience overthe years.I am sincerely thankful to my parents for the relentless encouragement andsupport, not only in the past five years of this work, but from the beginning. Iwould like to take a moment to extend my gratitude to all of my friends who havemade this journey even more enjoyable and memorable, especially Roland Hui,and of course, my juu, Yvonne Law.Last but not least, I would also like to acknowledge the financial support ofthe Natural Sciences and Engineering Research Council of Canada (NSERC) andRolls-Royce Canada. The majority of the computations in this work were per-formed on the Mammouth Paralle`le II supercomputer of Universite´ of Sherbrooke.xxviiChapter 1Introduction1.1 IntroductionAn abundance of energy is the prerequisite for the development of a flourishingcivilisation. Clean, renewable and sustainable energies, such as solar, wind, hydroand geothermal sources currently constitute approximately 4% of the total globalsupply [2]. Despite the benefits of the clean energies, the upfront cost of the relatedinfrastructures, relatively low power output and limited capacity for long-term stor-age hinder their general scalability. In consequence, fossil fuels will continue tobe the dominant source of energy for all modes of transportation, industrial powergeneration and domestic heating in the foreseeable future. As the global populationand standard of life continue to grow, it is evident that the accompanying demandfor energy will increase as well.The most convenient and effective approach to extract energy from fossil fuelsis through combustion, a process in which the chemical energy stored within themolecular bonds of the fuel is released primarily as heat. The problem, however, isthat combustion is not a perfect process. Pollutants such as carbon monoxide (CO),oxides of nitrogen (NOx), oxides of sulfur (SOx) and soot are inevitably formed.These substances are not only harmful to humans, but can also be detrimental to theenvironment [109]. Even in the most ideal scenario, the release of carbon dioxide(CO2) to the atmosphere is unavoidable as this greenhouse gas is the product of theoxidation of hydrocarbons with oxygen in the air [162].11.1 IntroductionWith these thoughts in mind, the world is confronted with two challenges re-lated to the usage of fossil fuels in the immediate future: (i) the conservation of fos-sil fuel supplies, which can be extended by developing more efficient combustiondevices; and (ii) the preservation of the atmosphere, which can be accomplishedby the reduction of pollutants generated from anthropogenic combustion processes.Both of these challenges spark an increasing interest in the research of sustainablecombustion technologies.Turbulent, premixed combustion is an important part of current heat engine de-signs; examples of these devices include the common automobile engine and sta-tionary gas turbines. The complex interactions between the turbulent flow, chem-ical kinetics and heat transfer at engine relevant pressures and temperatures es-tablish a highly non-linear problem spanning orders of magnitude in both lengthand time scales. The difficulty here is to introduce improvements to these highlycalibrated heat engines under such dynamic conditions.Computational fluid dynamics (CFD) is now indispensable in the developmentof complex engines due to its low cost and time requirement compared to exper-iments. While direct numerical simulation (DNS) is the most accurate approachto understanding the interactions between turbulence and flame, performing DNSfor engine relevant problems will remain impractical for the foreseeable future.With increasing computational capabilities, numerical simulation of these devicesvia the large eddy simulation (LES) paradigm is becoming more accessible to re-searchers. Nevertheless, due to the strong coupling between turbulence and chem-istry in premixed flames, the prediction of chemical reaction source terms continueto be a modelling challenge.The work in this thesis focuses on advancing two topics related to the simu-lation of turbulent, premixed combustion. First, new pseudo-turbulent probabilitydensity function (PDF) models of the reaction progress variable and scalar dis-sipation rate (SDR) models have been formulated. These PDF and SDR modelsare relevant to several approaches for premixed combustion; their applications arefurther explained in Chapters 3 and 4. Second, a modelling approach for the clo-sure of chemical source terms in turbulent premixed combustion simulations, theConditional Source-term Estimation (CSE), has been parallelised and improved forgreater efficiency and autonomy in complex geometries.21.2 Outline1.2 OutlineThis thesis is arranged such that the content of each chapter is self-contained; assuch, the relevant background information and the key acronyms are repeated inthe introduction of each chapter. In Chapter 2, a comprehensive introduction tothe mathematical tools and theoretical models necessary for the discussion of thenumerical simulation of turbulent, premixed combustion will be provided. Partic-ularly, this chapter focuses on the relevant aspects of chemistry, turbulence and theinteraction between these two highly coupled processes in premixed combustion.Furthermore, a number of modelling approaches for combustion simulation will bereviewed.In Chapter 3, a new method of constructing probability density function modelsfor premixed combustion is introduced. This method requires the use of a pseudo-turbulent, one-dimensional stochastic mixing model, referred to as the Linear-EddyModel (LEM). The methodologies are fully documented in this chapter and thenumerical results are compared to experimental data from a laboratory slot burner.In Chapter 4, the LEM has been implemented to not only formulate PDF mod-els for premixed combustion, but also to formulate the scalar dissipation rate mod-els; the SDR is an important quantity to both premixed and non-premixed com-bustion simulation strategies. The LEM-formulated PDF and SDR models arecompared to a number of turbulent, premixed, methane-air experimental flames,including swirling conditions. It is found that the PDF does not vary significantlywith changing turbulent intensities or swirl numbers. However, the SDR tendsto decrease in magnitude with increasing swirl and turbulent intensity. The LEMformulated models are well-matched with the experimental results in general.In Chapter 5, the Conditional Source-term Estimation combustion model, a clo-sure for chemical source terms in turbulent combustion simulations, has been fullyparallelised for optimised computations on parallel clusters. Large eddy simulationof a turbulent, premixed methane-air burner has been performed with updated CSEmodules. Comparison with the conventional CSE approach indicate that comput-ing times can be significantly reduced, in most cases, by approximately 50%, whileproducing nearly identical results.In Chapter 6, the LEM formulated PDF and SDR models are implemented31.2 Outlinein conjunction with the parallelised CSE module in a large eddy simulation of aturbulent, premixed methane-air, laboratory-scale burner. The LES result appearsto be well-matched with the experimental data, measured via Mie scattering, interms of the flame height and overall temperature field distribution. This modellingapproach is one step closer to being truly predictive as all of the tunable parametersare nearly eliminated.In Chapter 7, a summary of conclusions from each of the preceding chapterswill be highlighted, including brief discussions of the key accomplishments. Inaddition, future research topics related to the current work will be considered here.4Chapter 2Literature Survey2.1 BackgroundTurbulent combustion is a multi-scale phenomenon that requires a broad under-standing of fluid mechanics, chemical kinetics, heat transfer, molecular diffusionand other physical principles [181]. To perform numerical simulations of this com-plex problem, further knowledge regarding combustion models and ComputationalFluid Dynamics (CFD) is necessary. In this chapter, the fundamental mathematicaldescription of turbulence and flame will first be provided, including the governingequations for Newtonian fluid motion, thermodynamic properties, chemical kinet-ics and a discussion of the scales of turbulence. This will be followed by a review ofthe different paradigms of fluid simulation and an overview of combustion modelsfor turbulent, premixed combustion.2.2 Conservation EquationsThe derivation of the conservation equations for mass, momentum, energy andspecies mass fractions are well-established and can be found in many fluid mechan-ics textbooks [122, 132, 186]. These differential equations assume the continuumapproximation, whereby fluid particles can be viewed as continuously distributedwithin the corresponding volume without spatial voids [136]. For methane-aircombustion under atmospheric conditions, it is possible to further simplify the52.2 Conservation Equationsequations and assume negligible Soret1, Dufour2 and radiation effects [192]. Intensor notation, the time-dependent conservation equations for a reactive multi-component gaseous mixture involvingN species are commonly written as [132]:• Conservation of mass (commonly referred to as the continuity equation):∂ρ∂ t+∂ (ρui)∂xi= 0, (2.1)where ρ , t, ~u and ~x are respectively the mixture density, time coordinate,velocity and spatial coordinate. The Einstein notation is applied to the sub-scripts, i, j and k, in all of the conservation equations described in this sec-tion3.• Conservation of species:∂ (ρYα)∂ t+∂ (ρuiYα)∂xi+∂Ji,α∂xi= ω˙α , α = 1 . . .N , (2.2)where Yα and ω˙α are the mass fraction and chemical source term of speciesα . The mass diffusive flux of species α is,Ji,α =−ρDα ∂Yα∂xi , (2.3)where Dα is the molecular diffusivity of species α . Observing the conser-vation of mass and species equations, there are a total of N + 1 equationsfor the N species, resulting in an over-determined system. This is usuallyresolved by removing one species equation, typically an inert species witha large mass fraction in the mixture. For methane-air oxidation, a fittingcandidate is N2.1 Soret effect refers to the differential diffusion of particles with different molecular masses in thepresence of a temperature gradient within the system.2 Dufour effect refers to an enthalpy flux induced by a concentration gradient.3k represents a directional component of the Cartesian coordinate system when it is applied to asubscript. In other sections of this thesis, k represents the turbulent kinetic energy.62.2 Conservation Equations• Conservation of momentum:∂ (ρu j)∂ t+∂ (ρuiu j)∂xi+∂σi j∂xi= ρgi, (2.4)where σi j and ~g represent the stress tensor and the force of gravity, or bodyforce in general. The stress tensor is defined as a combination of the pressureand viscious stress tensors, σi j = pδi j− τi j. The viscous stress tensor reads,τi j = µ(∂ui∂x j+∂u j∂xi)− 23µδi j∂ukxk, (2.5)where µ , p and δi j are the dynamic viscosity, pressure and Kronecker delta,respectively.• Conservation of energy:∂ (ρE)∂ t+∂ (ρuiE)∂xi+∂qi∂xi+∂ (σi jui)∂x j= ρgiui+ Q˙, (2.6)where E, ~q and Q˙ are respectively the total energy, the total energy flux andthe heat source term. The heat source term is an external source, independentof the chemical energies. The energy flux is further defined as,qi =−Λ∂T∂xi −ρN∑α=1hαDα∂Yα∂x j, (2.7)where T , Λ and hα are respectively the temperature, thermal conductivity ofthe mixture and enthalpy of species α .To close the system of equations and obtain the pressure, the usual ideal gaslaw is required,p = ρRT = ρRuT/M, (2.8)where R, Ru and M are the mixture-averaged gas constant, the universal gas con-stant and the mixture-averaged molar mass.For completeness, the definition of the mass fraction of species α , Yα , and the72.3 Transport Propertiesmole fraction for species α , Xα , in a mixture ofN species are respectively,Yα =mα∑Nβ=1 mβ, (2.9)Xα =nα∑Nβ=1 nβ, (2.10)where mα is the total mass content of species α and nα is the total number of molesof species α within the mixture. A logical extension from the definitions of Yα andXα is that ∑Nα=1Yα = 1 and and ∑Nα=1 Xα = 1. In general, the mass fraction willbe written in the conservation equations while the mole fraction will be used in thenumerical computations by chemical software packages, for example, CHEMKIN-II [74]. Through basic manipulations, the mass and mole fractions can be convertedfrom one to the other from the following equations,Yα = XαMαM= XαMα∑Nβ=1 XβMβ, (2.11)Xα = MYαMα=Yα/Mα∑Nβ=1Yβ/Mβ, (2.12)where Mα is the molar mass of species α .2.3 Transport PropertiesTo correctly account for the transport of momentum, species and energy in a gaseousmixture, one must properly tabulate and evaluate the viscosities, thermal conduc-tivities, diffusion coefficients and thermal diffusion coefficients; these four quan-tities constitute the properties of transport in a mixture. In general, there are twocommon practical approaches to calculating these transport characteristics [73],the mixture-averaged and the multi-component methods. For many applications, itis possible to use the basic mixture-averaged properties while maintaining a highdegree of accuracy. In other cases, multi-component properties may be required.While some of these properties can be computed by first principles, such as thestandard kinetic theory expressions, vibrational modes, Lennard-Jones collisiondiameters and so forth [73], the quantities are typically tabulated as functions of82.3 Transport Propertiesthe temperature and pressure using high-order polynomials. The viscosity (µα )and thermal conductivity (Λα ) of each species are pressure-independent and aregenerally fitted with third order polynomials as functions of the temperature,ln µα =4∑i=1ai,α(ln T )i−1, (2.13)ln Λα =4∑i=1bi,α(ln T )i−1, (2.14)where ai,α and bi,α are the corresponding fitting coefficients. The binary diffusioncoefficients (Dαβ ) are similarly fitted,ln Dαβ =4∑i=1di,αβ (ln T )i−1, (2.15)where di,αβ are the corresponding fitting coefficients. However, this property ispressure-dependent. As such, these coefficients are computed at unit pressure andthe correct values are determined by dividing the evaluated coefficient from the fitby the actual pressure in the control volume.To evalutate the mixture-averaged viscosity, one can invoke the semi-empiricalformula of Wilke [188] and later modified by Bird et al. [16],µ =N∑α=1Xαµα∑Nβ=1 XβΨαβ, (2.16)where,Ψαβ =1√8(1+MαMβ)− 12 (1+(µαµβ) 12(MβMα) 14)2. (2.17)The mixture-averaged thermal conductivity can be computed with the combinationaveraging formula [108],Λ=12(N∑α=1XαΛα +1∑Nα=1 Xα/Λα). (2.18)92.4 ThermodynamicsThe mixture diffusion coefficient for species α , Dα , is,Dα =1−Yα∑Nβ 6=α Xβ/Dβα. (2.19)Equations 2.16, 2.18 and 2.19 provide the three necessary quantities for the mixture-averaged formulation. The thermal diffusion coefficients are only included for thespecies with low molecular mass, H to He.The detailed mathematical derivation of the multi-component formulation isbeyond the scope of this section, but its defining characteristics will be high-lighted. The multi-component formulation offers two advantages over the mixture-averaged transport properties. First, it does not have inherent limitations to thecomposition of the mixture qualities, such as the Fickian assumption, and thereforewill be more accurate. Second, multi-component transport ensures the conserva-tion of mass without the necessity of introducing correction factors. In exchangefor the improved accuracy, this model requires a substantial increase in computa-tional cost. An N ×N matrix must be inverted to obtain the multi-componentdiffusion coefficients. Further, to compute the thermal conductivity and thermaldiffusion coefficients, a 3N ×3N system of equations must be solved [73]. Un-less otherwise stated, mixture-averaged properties are used in the computation ofthe transport mechanics in this research.2.4 ThermodynamicsThe thermodynamic properties of chemical species are another important com-ponent to combustion modelling. In particular, these include the enthalpy (hα ),constant pressure heat capacity (cpα ) and entropy (sα ) of species α . Typicallythe thermodynamic properties are fitted through high order polynomials as func-tions of temperature over two separate ranges for improved accuracy. These coef-ficients are then compiled into lookup tables; examples of these tables include theChemkin Thermodynamic Data Base [75] and the database compiled by Gordonand McBride [59, 110]. For each species, there are seven tabulated coefficients at102.4 Thermodynamicseach of the two temperature ranges:cpαRα= c1+ c2T + c3T 2+ c4T 3+ c5T 4, (2.20)hαRαT= c1+c22T +c33T 2+c44T 3+c55T 4+c6T, (2.21)sαRα= c1ln(T )+ c2T +c32T 2+c43T 3+c54T 4+ c7. (2.22)This particular format is often referred to as the NASA polynomials [25].The mixture heat capacity and enthalpy can subsequently be calculated throughlinear combinations of the species specific properties, either by mass or mole frac-tions,cp =N∑α=1cpαYα , (2.23)h =N∑α=1hαYα . (2.24)The Chemkin Thermodynamic Data Base [75] has been used for all of the thermo-dynamic computations performed by this work unless specified otherwise.Aside from tabulating the thermodynamic properties for the purposes of numer-ical simulation, it is often useful to analytically relate the properties, particularly,hα =∫ T ∗Tre fcpα (T )dT +∆h0fα = hsα +∆h0fα , (2.25)where ∆h0fα is the species standard enthalpy of formation taken at a reference tem-perature, Tre f , typically 298 K. The integration of the species constant pressureheat capacity from Tre f to the current temperature, T ∗, yields the sensible enthalpyof the species, hsα . Futher, the total energy of the mixture (Equation 2.6) can becomputed from the mixture enthalpy as follows,E = h− pρ+12uiui. (2.26)112.5 Chemistry2.5 ChemistryOne can obtain valuable information related to the combustion process from ther-modynamics alone; however, this limits the calculations to global values such asenthalpy changes and adiabatic flame temperatures. For modern engine design,where the prediction of thermal efficiency and pollutant formation are equally crit-ical, thermodynamics alone is not sufficient to characterise the system. This iswhere a detailed understanding of the underlying chemical phenomenon allowsone to compute a crucial aspect of combustion — the time scales or rates at whichthe reactions occur. In the most accurate representation, the chemical reactions arefully coupled to the fluid flow, leading to the finite rate chemistry representation.2.5.1 Chemical KineticsFor methane-air combustion, finite rate chemistry was first introduced as a one-step, global reaction,CH4+2(O2 + 3.76N2)→ CO2+2H2O+7.52N2. (2.27)This representation, however, is not complete as a large number of minor speciesand radicals are formed and consumed during the combustion process via elemen-tary reactions. The GRI-Mech 3.0 [156] is a well-known example of a detailedfinite rate mechanism that represents methane-air oxidation. For this basic hydro-carbon, the GRI-Mech 3.0 chemistry mechanism prescribes 325 reactions and 53species to model the underlying phenomenon.The rate at which each of these elementary reactions occurs is dependent on theactivation energy required to change the reactants to products, the local concentra-tion of the reactants and the local temperature. In general, it has been establishedthat the Arrhenius law can adequately capture the reaction rate constant, kArr, formost elementary chemical reactions,kArr = A ·T Be−Ea/(RuT ), (2.28)where A, B and Ea are respectively the pre-exponential factor, temperature expo-nent and activation energy; these three parameters are determined experimentally.122.5 ChemistryTo obtain the overall reaction rate, one can multiply the reaction rate constant, kArr,with the concentrations of the reactants raised to the power of their stoichiometriccoefficient. Using a simple two-reactant reaction as an example,Reactant(A)+Reactant(B)→ Product(A)+Product(B), (2.29)the overall consumption rate of Reactant(A) would be,d[Reactant(A)]dt=−kArr · [Reactant(A)]1[Reactant(B)]1, (2.30)where [] represents the concentration of the species. For more complex mech-anisms, the overall rate of consumption or production of a species would be alinear combination of all elementary reactions in which the species participates;this follows the basic principles of the law of mass action. Indeed, there are anumber of specialised reactions that deviate from the standard Arrhenius law forcertain chemical mechanisms, including three-body reactions, fall-off reactions,chemically-activated reactions and pressure-dependent reactions [74].2.5.2 Reduced ChemistryWhile a detailed mechanism can accurately compute the reaction rates and de-scribe the underlying physics of the combustion process, it is much too expensivefor numerical simulations using the current generation of computational hardwarefor industrially relevant problems. As a result, a number of different modellingapproaches have been developed to overcome this obstacle.Detailed mechanisms can be reduced by the use of quasi steady-state andpartial-equilibrium assumptions [32, 125]. The resultant skeletal mechanisms wouldinclude fewer species and elementary reactions compared to the full chemistry,which would then decrease the demand on the computational hardware. In ex-change, however, the skeletal mechanisms are typically accurate only for a re-stricted range of equivalence ratios and flame conditions. Moreover, it is stillpossible for the time scales to differ greatly between reactions, leading to a stiffsystem [102]. Nevertheless, a large number of reduced mechanisms are avail-able for premixed, methane-air combustion due to their low computational require-132.6 Laminar Premixed Flamement [10, 11, 29, 114, 118, 146].A more sophisticated class of techniques for chemistry reduction are low-dimensional manifold approaches [104]. In these methods, the underlying conceptis to eliminate the slow time scales such that the system can be generalised bya lower-dimensional manifold in the composition space. Low-dimensional man-ifolds are well-suited for simulation purposes as the chemistry is fully tabulatedprior to the flow computation. This tabulation only needs to be done once foreach combination of thermodynamic and chemical composition. Both the chemi-cal source terms and species mass fractions are stored as a function of a few con-trolling variables for which transport equations are usually solved. Examples ofthese variables for methane-air combustion include H2O and CO2 mass fractions.The chemical properties are then retrieved from the table in real time, directingcomputing resources to the flow calculation instead of the reaction kinetics.The Flamelet-Generated Manifold (FGM) [121] is one such approach in whicha prototype flame is used to build a lower-dimensional manifold. In premixedcombustion, a one-dimensional laminar flame can be used to generate these tra-jectories. A combination of skeletal mechanism and FGM approaches has beenused to complete the work in this thesis; the skeletal mechanisms will be primar-ily used to construct the Probability Density Function (PDF) models for premixedcombustion in a pre-processing manner, whereas the FGM will be coupled directlyto the three-dimensional reactive flow calculations. The reduction techniques willbe discussed in greater detail in Chapters 4 and 5, respectively.2.6 Laminar Premixed FlameBefore the discussion of turbulent premixed flames, it is important to briefly reviewthe general structure of a laminar premixed flame. A premixed flame can be viewedas a moving front in which the exothermic reaction is sustained by the continuousconsumption of the nearby reactants. In the ideal planar laminar premixed flame,the flame front propagates steadily from the burnt mixtures towards the unburntreactants; the speed at which this front propagates divides the premixed flame intotwo categories, the subsonic deflagration and the supersonic detonation waves [91,96]. In this thesis, the discussion will be limited to the subsonic deflagration waves.142.6 Laminar Premixed FlameFrom a classical standpoint, as first proposed by Mallard and Le Chatelier andlater developed by Peters [128], the structure of a laminar premixed flame containsthree distinct layers: the preheat layer, the inner layer and the oxidation layer. Inthe preheat layer, the reaction rates are relatively low and the temperature increaseis caused primarily by conduction from the inner layer. This layer tends to bethicker than the other two due to the lack of chemical reactions in the governingmechanism. Next, the inner layer is where the fuel is consumed, radicals are de-pleted and the bulk of the heat of combustion is released. This layer tends to be thethinnest as the chemical time scales are relatively short. Last, the final products,accompanied by minor heat release, are formed through the slow reactions in theoxidation layer. Figure 2.1 illustrates the structure of a one-dimensional laminarpremixed flame. The three layers combine to form the flame front, which separatesthe unburnt fresh gas mixture from the burnt product mixture. The temperature ofthe mixtures increases across the flame front while the mixture density decreasessuch that the pressure changes are minimal across the flame.0 0.004 0.008 0.012 0.016 0.0200.10.20.30.40.50.60.70.80.91x (m)NormalisedTemperature|MassFraction PreheatlayerOxidationlayerT/TeqCH4O2CO2COInner layerFigure 2.1: One-dimensional laminar flame structure for premixed, stoichio-metric, methane-air combustion at standard conditions.152.6 Laminar Premixed FlameTwo important parameters are typically used to characterise a laminar premixedflame: the flame front propagation speed, commonly referred to as the laminarflame speed or the laminar burning velocity, sL, and the laminar flame thickness,δ f . The laminar flame speed describes the flame front propagation normal to itself.This property is dependent primarily on the equivalence ratio and temperature ofthe premixed reactants and secondarily on the overall pressure of the system. Thelaminar flame thickness is usually defined by either the diffusive or thermal thick-nesses. In this thesis, the thermal thickness is used exclusively unless otherwisespecified. The thermal thickness takes on the following form,δ f =Tb−TudTdx |max, (2.31)where Tb and Tu are respectively the burnt and unburnt gas temperatures; dTdx |maxrefers to the maximum gradient of the temperature across the flame front. Anotherparameter that characterises a premixed flame, which is commonly found in theliterature, is the thickness of the inner layer, lδ . This parameter is typically oneorder of magnitude smaller than the thermal thickness [128].For a laminar premixed flame at standard ambient temperature (298.15 K) andpressure (100 kPa), the ratio of fuel-to-oxidizer notably alters the laminar flamespeed, the overall flame thickness, and the equilibrium flame temperature. Thisratio is typically expressed as a non-dimensional quantity, referred to as the equiv-alence ratio (φ ), which compares the current fuel-to-oxidizer ratio to the stoichio-metric fuel-to-oxidizer ratio; stoichiometric combustion is achieved if the fuel andoxidizer are both perfectly consumed and converted to stable products in the reac-tion. Mathematically, the equivalence ratio is defined as,φ =mfuel/moxidizer(mfuel/moxidizer)stoich=nfuel/noxidizer(nfuel/noxidizer)stoich, (2.32)where mx and nx respectively represent the mass and number of moles of the quan-tity designated by x. Three possible scenarios can occur: (i) fuel-rich combustion,φ > 1; (ii) stoichiometric combustion, φ = 1; and (iii) fuel-lean combustion, φ < 1.Another commonly found dimensionless quantity is the oxidizer-to-fuel equiva-lence ratio, λ . This is simply the inverse of φ , where λ = φ−1. Figure 2.2 illustrates162.6 Laminar Premixed Flamethe laminar flame speeds and equilibrium temperatures of premixed methane-airflames at different equivalence ratios in standard ambient temperature and pressureconditions. The values are calculated by Cantera [58] using the detailed GRI 3.0mechanism [156] and mixture-averaged transport.0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.500.040.080.120.160.20.240.280.320.360.4φLaminar Flame Speed (m/s)13001400150016001700180019002000210022002300Equilibrium Flame Temperature (K)Figure 2.2: Laminar flame speeds and equilibrium temperatures of premixedmethane-air flames at different equivalence ratios in standard ambient temper-ature (298.15 K) and pressure (100 kPa) conditions.To further expand on the equivalence ratios, lean premixed combustion is gen-erally the most efficient as the excess oxidizer serves to exhaust the fuel supplywhile driving the reaction to completion. If the process is controlled appropriately,it can also reduce pollutant formation, such as NOx, due to the reduced combustiontemperatures. However, the shortcoming is a reduced flame speed, which can leadto a lower power output for the same engine dimensions and intake conditions.Stoichiometric combustion is typically found in spark-ignited automobile engines,as it yields a suitable balance between power, efficiency and pollutant formation.Rich combustion is not frequently desired as it results in Unburnt Hydrocarbons(UHC) due to the excess fuel. This leads to a reduction in thermal efficiency andthe UHC would become an additional source of pollutants. It is important to notethat the flow conditions would be turbulent for engine-relevant combustions ratherthan laminar, though the rich, stoichiometric and lean combustion generalisationsstill stand.172.7 Scales of Turbulence2.7 Scales of TurbulenceIn practical combustion applications, the fluid motion predominately exists in aturbulent rather than a laminar regime. A flow is considered turbulent if it can becharacterised as unsteady, three-dimensional, chaotic, vortical and dissipative [136,186]. This highly non-linear, multi-scale phenomena is, in principle, fully de-scribed by the equations shown in Section 2.2. The conservation of momentumequations are referred to as the Navier-Stokes equations, which were first derivedin the early 19th century by G. Stokes and M. Navier, independently. They ele-gantly describe the coupling between the density, velocity, temperature and pres-sure within the fluid. Coincidentally, due to the strong non-linearity of this systemof equations, analytical solutions currently do not exist; such analytical solutionscan only be derived for laminar conditions or highly simplified systems. Indeed,simply providing the mathematical proof for the existance and uniqueness of a so-lution to the Navier-Stokes equations is one of the most challenging theoreticalproblems of the 21st century [48].In the absence of a robust analytical framework, turbulence research has beendependent on experiments and numerical simulations. Through observations inthe late 19th century by Reynolds [138], a fundamental non-dimensional parameterthat characterises the intensity of turbulent flows emerged, the Reynolds number,Re. It is the ratio of inertial forces to viscous forces and is defined as,Re≡ U ·Lν, (2.33)where U , L and ν are respectively the characteristic velocity, the characteristiclength and kinematic viscosity of the fluid. The bulk flow velocity would typi-cally be used as U , while a suitable length scale determined by the geometry ofthe flow field would be used as L. In the internal pipe flow experiments performedby Reynolds [138], the transition between laminar and turbulent states occured at aReynolds number between 2,300 to 4,000. For external boundary layer flows overflat plates, however, the transition between laminar and turbulent regimes typicallyoccur at a Reynolds number of 300,000 to 500,000. Nevertheless, despite havingrather extensive analytical understanding of laminar flows and empirical relationsfor turbulent flows, the transition between the two regimes is still not well charac-182.7 Scales of Turbulenceterised [5].Two parameters that are typically used in the literature to specify the charac-teristics of a turbulent flow field are the turbulent kinetic energy, k, and the integrallength scale, l0. By the Reynolds decomposition, which will be discussed in greaterdetail in Section 2.9, the instantaneous velocity, ui, can be separated into a meanvelocity, ui, and a fluctuating component, u′i. The turbulent kinetic energy can thenbe defined as,k =12((u′x)2+(u′y)2+(u′z)2), (2.34)where u′x, u′y, u′z are respectively the fluctuating components of the velocity in eachdirection of the Cartesian coordinate system. For isotropic4, homogeneous5 turbu-lent fields, it is possible to further simplify the turbulent kinetic energy to,k =32u′2, (2.35)where u′ is the root-mean-square of the velocity fluctuations, also commonly knownas the turbulent intensity. The integral length scale indicates the size of the turbu-lent motions, or turbulent eddies, in which the majority of the turbulent kineticenergy manifests [136]. It is usually considered the large-scale motion of the tur-bulent field and is typically on the same order of magnitude as the characteristiclength, L. To formally characterise the spatial structure of a random field, one canintroduce a two-point correlation, or autocovariance, defined as [136],Ri j(~r,~x, t)≡ 〈ui(~x, t)u j(~x+~r, t)〉, (2.36)where ~r is some distance away from ~x. The integral length scale related to thex1-coordinate direction can then be written as,L11(~x, t)≡ 1R11(0,~x, t)∫ ∞0R11(eˆ1r,~x, t)dr, (2.37)where eˆ1 is the unit vector in the x1-coordinate direction. Indeed, one can definean integral length scale for each direction of the coordinate system. For the special4Statistically invariant under rotations and reflections of the coordinate system.5Statistically invariant under translations of the coordinate system.192.7 Scales of Turbulencecase of isotropic homogeneous turbulence, one can associate the integral lengthscale, l0, with any one of the Lii. From these parameters, one is able to defineanother non-dimensional number, the turbulent Reynolds number, Ret ,Ret =u′l0ν. (2.38)As a first approximation, the turbulent Reynolds number is generally an order ofmagnitude lower than the Reynolds number for the same flow field.Aside from the non-dimensional numbers, it has also been observed that theproduction mechanism driving the turbulent motions, or turbulent eddies, trans-fers its kinetic energy from the largest characteristic scales to the smallest scalesuntil it is finally dissipated by viscous action; this concept of energy transfer is re-ferred to as the energy cascade [139]. Kolmogorov later formulated two similarityhypotheses pertaining to the different ranges of the energy cascade. The first ofthese postulated that the statistics of the small-scale motions can be ubiquitouslydetermined by the kinematic viscosity of the fluid, ν , and the rate of dissipationof the turbulent kinetic energy, ε , provided the turbulent flow is at a sufficientlyhigh Reynolds number [89]. Through dimensional analysis, three parameters canbe derived from the postulation describing the smallest scales of turbulence; theseare respectively named the Kolmogorov length (ηK), time (τK) and velocity (uK),ηK ≡ (ν3/ε)1/4, (2.39)τK ≡ (ν/ε)1/2, (2.40)uK ≡ (εν)1/4. (2.41)The second of Kolmogorov’s hypotheses postulated that the statistics of the mo-tions of scale l in the range ηK l l0 can be uniquely determined by the rate ofdissipation of the turbulent kinetic energy, ε , and independent of the kinematic vis-cosity of the fluid, ν , provided the turbulent flow is at a sufficiently high Reynoldsnumber [89]. This range of length scales is referred to as the inertial subrange. Theinertial subrange exhibits a unique quality where the energy carried by a turbulenteddy of size l decreases with increasing wavenumber, κ (κ = 2pi/l), with a known202.7 Scales of Turbulencefunctional form. Particularly, the energy spectrum in the inertial subrange is,E(κ) =Cκε2/3κ−5/3, (2.42)where Cκ is a universal constant [136]. Figure 2.3 illustrates a typical energy spec-trum observed in turbulent motion at high Reynolds numbers.Figure 2.3: A schematic diagram of the turbulent kinetic energy spectrum athigh Reynolds numbers.A very informative relationship can be established between the integral scalesand the Kolmogorov scales through an estimation of the dissipation rate, ε , interms of the large-scale features of the flow. The ε is approximately equal to thekinetic energy production rate because of the transfer of energy through the energycascade [136]. The kinetic energy is proportional to U2 while the time scale of theenergy transfer is estimated to be on the order of the large-scale motion, l0/U . Thedissipation rate can thus be scaled as,ε ≈ U2l0/U≈ U3l0. (2.43)212.8 Regimes of Premixed Turbulent CombustionSubstituting this relation into Equation 2.39, one obtains,ηK =(ν3ε)1/4≈(ν3l0U3)1/4. (2.44)Finally, the ratio between the integral length scale and Kolmogorov length scalecan be approximated,l0ηK≈ l0U3/4ν3/4l1/40≈(Ul0ν)3/4≈ Re3/4. (2.45)Following a similar derivation, one can additionally obtain the relationship betweenthe integral time scale, τt , and the Kolmogorov time scale, τK ,τtτK=l0/U(ν/ε)1/2≈ l0/U(ν l0/U3)1/2≈(Ul0ν)1/2≈ Re1/2. (2.46)Equations 2.45 and 2.46 suggest that the disparity between the largest and small-est turbulent scales of the flow field increases rapidly as the Reynolds number in-creases.2.8 Regimes of Premixed Turbulent CombustionThe overall behaviour of the premixed flame would differ under various turbu-lent conditions as a result of the changing turbulence-chemistry interactions. Ithas been recognised that these interactions can be categorised into large and smallscales [36]. Figure 2.4 provides a conceptual overview to the behaviour of a pre-mixed flame front exposed to large-scale turbulence. At very low turbulent inten-sities, the eddies typically do not alter the overall shape and thickness of the pre-mixed flame. The flame front retains its planar geometry and propagates towardsthe reactants at approximately the laminar flame speed, sL. As the turbulent inten-sity increases, the premixed flame would become corrugated due to the interactionwith the vortical motion of the eddies. The total surface area of the flame frontwould increase due to the corrugation, resulting in a higher rate of consumption ofthe reactants. This increased effective flame speed is referred to as the turbulent222.8 Regimes of Premixed Turbulent Combustionflame speed, sT . As proposed by Damko¨hler [36],sTAT=sLAL, (2.47)where AT and AL are the turbulent flame surface area and the laminar flame surfacearea. The ratio of AT /AL is referred to as the flame wrinkling factor, Ξ. As the tur-bulent Reynolds number increases, the kinetic energy from the large-scale motionswould be transported to smaller and smaller eddies; such eddies are then able to in-teract with the flame front more effectively. Past some threshold, it is even possiblefor the smallest eddies to not just wrinkle the flame front, but the Kolmogorov-scale turbulent motion could also affect the transport mechanisms within the flame— this is the regime of small-scale turbulent interactions.Increasing turbulent intensityFigure 2.4: Effect of large-scale turbulent eddies on a premixed flame front.This complex coupling between the chemistry and turbulence leads to a sep-aration of combustion characteristics dependent on the physical parameters thatare representative of the system. The length and time scales of the turbulent flowfield and the chemical reactions then provide a basis for the identification of dif-ferent regimes for turbulent premixed combustion; the analysis of which is bestsummarised by regime diagrams, as originally proposed by Borghi [19]. A numberof such diagrams have been proposed to characterise the turbulence-chemistry in-teractions for turbulent premixed combustion, including the work of Abdel-Gayedand Bradley [3], Borghi [19], Bray [23], Du¨sing et al. [46], Peters [126, 127],Pitsch and Duchamp de Lageneste [131] and Poinsot et al. [133]. Figure 2.5 illus-trates a commonly referenced regime diagram in the literature, the modified Borghidiagram, as proposed by Peters [127]. The ratio of the turbulence intensity to the232.8 Regimes of Premixed Turbulent Combustionlaminar flame speed is plotted against the ratio of the integral length scale over thelaminar flame thickness on logarithmic scales.10−1 100 101 102 10310−1100101102103u′/sLl0/δfBroken reaction zonesThin reaction zoneCorrugated flameletsWrinkled flameletsLaminarflamesKaδ = 1Ka= 1Da= 1Ret = 1Figure 2.5: Modified Borghi diagram illustrating different turbulent premixedcombustion regimes.Prior to the detailed discussion of the combustion regimes illustrated by themodified Borghi diagram, it is important to acknowledge a few additional dimen-sionless parameters aside from the turbulent Reynolds number. Scaling argumentsbetween the ratios of turbulence and chemistry lead to the definition of the follow-ing entities:• The Damko¨hler number, Da, which is the ratio between the turbulent integraltime scale to the chemical time scale,Da≡ τtτc=l0/u′δ f /sL. (2.48)• The First Karlovitz number, Ka, which is the ratio between the chemical242.8 Regimes of Premixed Turbulent Combustiontime scale to the Kolmogorov time scale,Ka≡ τcτK=δ f /sLηK/uK=δ 2fη2K, (2.49)where by dimensional analysis,ηK =(ν3ε)1/4, sL ≈ νδ f , uK = (εν)1/4. (2.50)• The Second Karlovitz number, Kaδ , which is the ratio between the chemicaltime scale, based on the inner layer, to the Kolmogorov time scale,Kaδ ≡l2δη2K. (2.51)Assuming the thickness of the inner layer is approximately 10% of the lam-inar flame thickness [128], then Kaδ ≈ 0.01 ·Ka.• The turbulent Reynolds number, Ret , which has been previously defined inEquation 2.38, can be rewritten in terms of the Damko¨hler and Karlovitznumbers,Ret = Da2Ka2. (2.52)Armed with the preceding dimensionless parameters, one can observe five dis-tinct regimes of turbulent premixed combustion within the modified Borghi dia-gram:• The laminar flame regime: Ret < 1. The turbulent intensity is at a negligiblelevel. The flame behaviour can be fully encapsulated by the description ofthe laminar premixed flame in Section 2.6.• The wrinkled flamelet regime: Ret > 1, Ka< 1 and u′ < sL. Since the chem-ical time scale is shorter than the integral turbulent time scale, the internalflame structure remains unaffected by the turbulent motion. More precisely,the thickness of flame front is smaller than the smallest eddies, such that252.8 Regimes of Premixed Turbulent Combustionthe laminar structure remains unchanged. The turbulent motion would onlyserve to transport the flame on a large scale.• The corrugated flamelet regime: Ret > 1, Ka< 1 and u′ > sL. This is similarto the wrinkled flamelet regime, except the turbulent intensity is at a higherlevel. The increased turbulent motion may lead to the formation of isolatedpockets of fresh and burnt gases, hence the name. The smallest eddies arestill larger than the flame front thickness, allowing the internal flame struc-ture to remain mostly laminar.• The thin reaction zone: Ret > 1 and 1 < Ka < 100. The Kolmogorov lengthscale is shorter than the overall flame thickness such that the smallest ed-dies are able to penetrate into the preheat layer, though not into the innerlayer [132]. Part of the internal flame structure can be modified by the ed-dies and therefore should no longer be modelled with locally laminar as-sumptions.• The broken reaction zone: Ka > 100. The Kolmogorov length scale isshorter than the thickness of the inner layer such that the smallest eddiesare able to penetrate into the inner layer and modify the flame structure. Theturbulent mixing is extremely fast and the reaction rate is limited by chem-istry [132].While the turbulence-chemistry interaction is often described from the per-spective of eddies penetrating into the various layers of the flame, it is importantto realise that the flame correspondingly modifies the turbulent field. For example,the temperature across the premixed flame front can increase by a factor of eightor so (Figure 2.2), whereas the kinematic viscosity of air varies with temperatureapproximately as ν ∼ T 1.7. This results in a significant decrease in the Reynoldsnumber in the product mixture, which may lead to relaminarisation of the locallyturbulent field [132]. Moreover, a rapid increase in gas temperatures in the thinfront leads to a rapid increase in gas velocities in that region. This accelerationin turn modifies the turbulent field, which results in an effect referred to as theflame-generated turbulence [132]. The intricate feedback between turbulence and262.9 Computational Fluid Dynamics (CFD)chemistry is one of the fundamental reasons why turbulent combustion modellingwill continue to be a challenge in the foreseeable future.Practical combustion devices, such as gas turbines and automotive internalcombustion engines, often operate within the corrugated flamelet regime and thinreaction zone, resulting in an increased interest for models that are representative ofthese regimes. A brief review to the turbulent combustion models available todaywill be discussed in Section 2.10.2.9 Computational Fluid Dynamics (CFD)The complexity of the governing conservation equations suggests that the analyt-ical solutions to any practical fluid dynamics problems involving turbulence maynot be derived for some time [48]. This limits researchers to two possible ap-proaches in the quest of understanding the nature of turbulence: experimental ob-servations and computation of the governing equations using numerical solvers,which is the study of computational fluid dynamics. Though the results from awell-designed experiment are highly accurate and reliable, these are often costlyand the desired metrics may not be readily measured in combustion devices. More-over, simply introducing a measuring tool to the system may alter the character-istics of the device. Computational fluid dynamics, thus, emerges as an impor-tant tool for understanding the effects of turbulence and flame from an theoret-ical perspective and for predicting the power and pollutant outputs in practicalcombustion devices. Three general paradigms of CFD currently exist: Direct Nu-merical Simulation (DNS), Large Eddy Simulation (LES) and Reynolds-AveragedNavier-Stokes (RANS). Each of these strategies has certain benefits over one an-other [34, 141] and their characteristics will be briefly summarised in this section.DNS, as the name implies, directly resolves all of the spatial and temporalscales of the turbulent field without subjecting the simulation to any modellingapproaches. This necessarily requires the computational grid to resolve the threedimensions of space and one dimension of time at the Kolmogorov scales. Moreprecisely, the grid spacing in each direction of the simulation domain, ∆6, mustbe on the order of ηK , while the size of the timesteps must be on the order of τK .6 Assuming the grid spacing is uniform in x, y and z directions.272.9 Computational Fluid Dynamics (CFD)Referring to Equations 2.45 and 2.46, one can roughly estimate the computationalrequirement of DNS to scale approximately as Re11/4 [136]. The tremendous com-putational effort renders DNS unfeasible for practical applications as the turbulentReynolds number can exceed well beyond 105 [112]. Another difficulty is the im-plementation of the correct initial and boundary conditions; if these are inaccurate,then even a perfect grid would be unable to provide an accurate solution. DNSis presently a useful research tool for studying the nature of turbulence, for bothreacting and non-reacting flows, in relatively low Reynolds number environments.Moreover, the detailed data obtained from the DNS can also be used to validateother simulation models [99, 112].The restricted applicability of DNS motivates researchers to introduce meth-ods that can decrease the computational cost of fluid simulations. This can oftenbe achieved by reducing the degrees of freedom of the problem by solving foraveraged quantities. The most intuitive approach is to compute the time-averagedsolution of the Navier-Stokes — this method is the well-known Reynolds-AveragedNavier-Stokes simulation paradigm. By the Reynolds decomposition, the instan-taneous velocity, ui, can be separated into a mean velocity, ui, and a fluctuatingcomponent, u′i,ui = ui+u′i. (2.53)The mean velocity is defined as either,ui =1t∗∫ t∗0uidt, (2.54)or,ui =∑Mj=1 ui, jM. (2.55)In the former equation, t∗ is the period over which the average should be performed,where t∗ must be larger than the time scale of the large-scale motion, l0/U7. In thelatter equation, ui, j is the velocity in the jth realisation in an ensemble ofM realisa-tions, whereM is taken to be sufficiently large. This suggests that while the meanfield can be obtained from RANS, all of the turbulent features are unresolved and7This time scale is also commonly referred to as the large-eddy turnover time.282.9 Computational Fluid Dynamics (CFD)require modelling. Figure 2.6 illustrates a symbolic representation of the resolvedand modelled energy spectrum of the solution in RANS. The statistical character-istics of RANS prevent the simulation from isolating transient turbulent structuresand therefore it is not a useful tool for studying the fundamental behaviour of tur-bulence. However, the general thermodynamic performance of combustion devicescan be obtained if the turbulent fluctuation is properly modelled. As the computa-tional cost of RANS is significantly lower than that of DNS and LES, it is still apreferred method for industrial applications.E(κ)κResolved+E(κ)κModelled=E(κ)κTotalFigure 2.6: A schematic representation of the resolved and modelled energyspectrum in RANS.Prior to deriving the RANS governing equations, it is useful to introduce adensity-weighted averaging method referred to as Favre averaging, which servesthe purpose of simplifying the governing equations by eliminating the fluctuatingterms related the density in compressible flows. It is defined as,u˜i ≡ ρuiρ . (2.56)The Reynolds decomposition can then be restated in terms of the Favre notation,ui = u˜i+u′′i , (2.57)with u′′i being the fluctuations around the Favre average. Substituting into Equa-292.9 Computational Fluid Dynamics (CFD)tion 2.4, one can obtain the Favre-averaged conservation of momentum equation8,∂ρ u˜ j∂ t+∂∂xi[ρ u˜iu˜ j + ρu′′i u′′j︸ ︷︷ ︸unclosed+pδi j− (τ˜i j + τ ′′i j︸︷︷︸unclosed)]= ρgi. (2.58)In general, τ ′′i j can be safely neglected for most flows as |τ˜i j| |τ ′′i j|. The ρu′′i u′′jrepresents the Reynolds stress tensor, which remains unclosed. Indeed, providingsuitable models for this quantity is the fundamental study of RANS. Many closuremodels for the Reynolds stress tensor have been proposed in the past; each hasexcelling features in certain applications in comparison with the others. Thesemodels include: (i) the one-equation formulation of Spalart-Allmaras [159]; (ii)the two-equation formulations, such as the k-ε model of Launder [94] and the k-ωmodel of Wilcox [187]; (iii) the stress-transport models of Launder [95]; and (iv)the algebraic models of Baldwin-Lomax [6].As computational hardware became more powerful in the past decades, a mod-ification to the traditional RANS approach, referred to as the Unsteady Reynolds-Averaged Navier-Stokes (URANS), had been proposed. A change in formulation inthis method allows for some of the low-frequency modes, usually in the range of afew hundred hertz, to be resolved [141]. It has a number of different names, includ-ing, Semi-Deterministic Simulation (SDS), Very Large Eddy Simulation (VLES)and Coherent Structure Capture (CSC) [9, 179]. Though this unsteady formula-tion contains more information than the usual RANS, the majority of the wavemodes are still modelled; as a result, this approach shares many of the underlyingdeficiencies as RANS.It appears, then, a shift in the fundamental averaging process is necessary tocapture the turbulent motion in fluid flows. A solution is to resolve the large-scale structures and to model the small-scale structures in the flow field by theapplication of a spatial filter — this is the underlying principle of large eddy sim-ulation [141]. In turbulent flows, features such as mixing, instabilities and certaintransient effects that are governed by the large-scale structures can be computedexactly by LES. Furthermore, modelling the fine-scale structures is justified bythe first hypothesis of Kolmogorov, where the small turbulent motions tend to have8 This process can be repeated similarly for the continuity and energy equations.302.9 Computational Fluid Dynamics (CFD)universal characteristics dependent only on ν and ε of the fluid, especially at higherReynolds numbers [89]. The unresolved flow features that are too small for the gridare commonly described as Subgrid-Scale (SGS) or Subfilter-Scale (SFS) struc-tures. A two-dimensional schematic view of the resolved and subgrid-scale flowstructures in a grid with spacing of ∆ is shown in Figure 2.7.Figure 2.7: Schematic view of the resolved and subgrid-scale structures ina grid with spacing of ∆: Resolved structures (light gray) and subgrid-scalestructures (dark gray).The grid spacing, ∆, acts as an implicit filter that eliminates all of the high-frequency wave modes from direct computation. For anisotropic grids, the filterwidth is calculated as [37],∆= (∆1∆2∆3)1/3. (2.59)In practice, the filter width is typically assumed to be proportional to the local cellsize and does not need to be exactly equal. Figure 2.8 symbolically demonstratesthe separation between the resolved and modelled scales in the energy spectrumfor LES. The boundary is solely established by the filter width; as one decreases∆, one can resolve more high-frequency modes at the expense of increasing thecomputational cost. Further discussion on the topic of implicit filtering can befound in [61, 141, 150, 151].312.9 Computational Fluid Dynamics (CFD)E(κ)κResolved+E(κ)κModelled=E(κ)κTotalFigure 2.8: A symbolic representation of the resolved and modelled turbulentkinetic energy spectrum in LES.While the implicit filter is commonly found in LES implementations becauseof its convenience and its ability to take complete advantage of the grid resolution,it includes the unavoidable drawbacks of truncation errors and a filter functionthat cannot be exactly determined. This may be of issue as wave modes are notidentically filtered in different regions of the grid. For explicit LES, the spatialfiltering is defined as [97],〈ui(~x, t)〉=∫VG(~x−~x′)ui(~x′, t)d~x′, (2.60)where G(~x) and 〈ui〉 are respectively the low-pass filter and the filtered velocitycomponent. The three classical filters commonly found in the literature are: (i) thespectral filter, where all of the wave numbers κ > pi/∆ are eliminated; (ii) the top-hat filter, where a volume average over a sphere of radius ∆/2 is performed [136];and (iii) the Gaussian filter, where it represents the joint normal distribution with amean of zero and a covariance of ∆2/12. Each of these filter functions satisfies thefollowing normalisation condition,∫G(~x)d~x = 1. (2.61)An additional complexity arises in the spatial filtering introduced by LES thatdoes not appear in RANS. The differentiation and filtering operations do not nec-essarily commute if a non-uniform filter is applied to the governing equations. Ithas been found that the commutation errors are of second order, O(∆2) [54]. Assuch, it may be unnecessary to employ discretisation schemes higher than order322.9 Computational Fluid Dynamics (CFD)two because the commutation errors would dominate.One can also derive the momentum equation in the context of spatial filteringand Favre averaging for LES by first redefining u˜i, which now represents the Favre-filtered velocity component,u˜i ≡ 〈ρui〉〈ρ〉 . (2.62)Substituting into Equation 2.4, one can obtain the Favre-filtered conservation ofmomentum equation9,〈ρ〉∂ u˜ j∂ t+∂∂xi〈ρ〉u˜iu˜ j + 〈p〉δi j− τ˜i j + σ sgsi j︸︷︷︸unclosed= 〈ρ〉gi, (2.63)where σ sgsi j ≡ 〈ρ〉(u˜iu j− u˜iu˜ j) is the subgrid-scale turbulent stress tensor and τ˜i j isthe viscous stress tensor defined as,τ˜i j = 〈µ〉(∂ u˜i∂x j+∂ u˜ j∂xi)− 23〈µ〉δi j ∂ u˜kxk . (2.64)To provide closure for the SGS turbulent stress tensor, two main classes ofmodelling approaches are available, the functional modelling and the structuralmodelling. In the functional approach, the unclosed, subgrid scales are modelledaccording to the resolved scales. Some examples include the one-equation eddy-viscosity model [149]10, the Smagorinsky model [155], the dynamic Smagorinskymodel [53] and a number of variants of the preceding models [54, 100, 111]. Thesemodels are based on solving equations for a turbulent viscosity, νT , which can thenbe related to σ sgsi j . In the structural approach, the subgrid scales are computedwithout directly relating the resolved fields to the SGS quantities. Examples in-clude the scale-similarity model [7], the differential stress equation model [38], thenon-linear anisotropic model [103] and so forth [141].9 This process can be repeated similarly for the continuity and energy equations. The closure ofthe additional fluctuating terms can be found in [141].10 The one-equation eddy-viscosity model is used for all of the LES calculations in this thesis.332.10 Turbulent Premixed Combustion Modelling2.10 Turbulent Premixed Combustion ModellingFirst moment closure, whereby the filtered or averaged reaction rates are directlycomputed by the filtered or averaged values of the scalar fields via the Arrheniusrate equations, can lead to considerable discrepencies compared to the fully re-solved values, due to the highly non-linear interactions between chemical sourceterms and the fluctuating scalar fields. The chemically reactive regions are con-fined to the inner layer, which exists at a smaller length scale than the average LESgrid resolution. In consequence, combustion models are required for LES, which,by their nature, cannot fully resolve the turbulence-chemistry interactions. A num-ber of different combustion models have been proposed for LES; many of theseare derived from concepts developed for RANS approaches. While the remainderof this section will present brief descriptions of some of the most commonly citedcombustion models, it is important to mention that the modelling of the filteredreaction rates remains a topic of significant research in the combustion community.Eddy Break-Up ModelOne of the simplest closures for chemical source terms is an algebraic approachreferred to as the Eddy Break-Up (EBU) model, first introduced by Spalding [160].This model assumes fast chemistry and thus the rate determining process for thechemical source terms would be completely governed by the rate of mixing. Thechemical source term under such turbulent conditions for the products (ω˙P) areestimated to be proportional to the turbulent time scale and the variance of massfraction of the product (Y ′2P ),ω˙P =CEBUρεkY ′2P , (2.65)where CEBU is the Eddy Break-Up constant.Although the EBU model can be successfully applied to some RANS simu-lations, more accurate closures are desired for LES [181]. The problem here isthat the EBU model does not take into account any aspect of chemistry; it can nei-ther distinguish between different equivalence ratios nor identify different types offuels. Indeed, the reaction rate is linked directly to the tunable Eddy Break-Up con-342.10 Turbulent Premixed Combustion Modellingstant, which must be adjusted on a case-by-case basis. As a result, the chemistryeffects are completely ad hoc and non-physical results may ensue [181].An improvement to the EBU model was made subsequently by Magnussen andHjertager [105], termed the Eddy Dissipation Model (EDM). The modification re-places the variance of mass fraction of the product by the mean mass fraction ofthe deficient species, for example, fuel in lean mixtures or oxygen in rich mixtures.This allows the model to distinguish between rich and lean mixtures; however, itshares the same drawbacks as the EBU in that the chemistry effects are completelyad hoc and the constant, CEDM for this model, would require adjustment for eachnew application.Bray Moss Libby (BML) ModelThe Bray Moss Libby (BML) model, similar to the EBU and EDM models, isan algebraic approach first introduced by Bray, Moss and Libby [20, 98]. Whileit has been improved a number of times since its introduction [21, 24], the BMLmodel fundamentally assumes the flame to be infinitely thin, where the progressvariable would intermittently take on a value of 0 or 1. This model generally usesa temperature-based progress variable,c =T −TuTb−Tu , (2.66)where T , Tb and Tu are respectively the local temperature, the burnt gas temperatureand the unburnt gas temperature; this formulation is also referred to as the reducedtemperature. Under the thin-flame assumption, the PDF of the progress variablecan be reduced to the following form with two Dirac delta functions,P(c) = αBMLδ (c)+βBMLδ (1− c), (2.67)where αBML +βBML = 1. Combining this PDF and the scalar dissipation rate for-mulation derived by Bray and Moss [20], it is possible to obtain the average reac-352.10 Turbulent Premixed Combustion Modellingtion rate of the progress variable,ω˙c =ρχ˜c2CBML−1 , (2.68)where χ˜c is the scalar dissipation rate of the progress variable. Furthermore, theparameter, CBML, is defined as [20],CBML =cω˙cω˙c. (2.69)One method to obtain an approximation to the scalar dissipation rate would be asimple balance equation [106],ρχ˜c = ρχc ≈ ρ c˜′2τBML, (2.70)where c˜′2 and τBML are the variance of the progress variable and the turbulentmixing time scale in the BML model, respectively. Finally, if one assumes theτBML to be related to the turbulent kinetic energy, k, and the dissipation rate, ε , asτBML = k/ε , one can recover the following equation for the production rate of theprogress variable,ω˙c =12CBML−1ρεkc˜′2. (2.71)Indeed, a brief comparison of the preceding equation with Equation 2.65 wouldreveal that the BML model is closely related to the EBU model; some view theBML model as a theoretical derivation of the EBU model, except it is superior inthat the CBML parameter partially captures the effects of chemistry. However, muchlike its counterpart, the BML model is not able to account for the effects of chem-istry accurately.Flame Surface Density (FSD) ModelThe Flame Surface Density (FSD) model is a flamelet approach that assumesthe flame front is a surface that can be convected, diffused and strained by theturbulent interactions [64]. Assuming c is the usual progress variable, where it can362.10 Turbulent Premixed Combustion Modellingbe based on the temperature or a mass fraction:c =T −TuTb−Tu , (2.72)c =YCO2−Y uCO2Y bCO2−Y uCO2, (2.73)where YCO2 , YuCO2 and YbCO2 are respectively the local, unburnt and burnt CO2 massfractions11. By definition, c= 0 in reactants and c= 1 in products, and the transportequation for c can be established as,∂ (ρc)∂ t+∂ (ρuic)∂xi=∂∂xi(ρDc∂c∂xi)+ ω˙c = ρsd |∇c|. (2.74)Here, sd is the local displacement speed of the iso-surface of the progress variable.One can filter the preceeding equation and obtain a term related to the displacementof the iso-surface, ρsd |∇c|, which requires closure. The concept behind the FSDmodel is to relate this displacement term to the flame surface density,ρsd |∇c| ≈ ρusLΣ= ρusLΞ∆|∇c|, (2.75)where ρu, Σ and Ξ∆ are respectively the unburnt mixture density, the SGS flamesurface density and the SGS wrinkling factor. The closure of Σ has been proposedin the form of transport equations [17, 64], algebraic formulations [17] and a scale-similarity model [87]. Furthermore, the closure of Ξ∆ can be found in the form oftransport equations [50, 184] and algebraic expressions [30, 31, 51]. Understand-ably, the accurate closure of the SGS terms and the modelling of the effects ofstrain on the flame surface constitute some of the challenges in applying the FSDmodel in LES calculations [47, 130, 132].11 The mass fraction CO2 can be used as the progress variable for lean to stoichiometric methane-air combustion under atmospheric conditions.372.10 Turbulent Premixed Combustion ModellingG-Equation ModelThe G-Equation model was introduced by Markstein [107] and is based on alevel set approach in the calculation of a G-field. G typically takes on a value ofless than 0 in the unburnt mixture, 0 at the flame front and greater than 0 in theburnt mixture. The transport equation for G is,∂G∂ t+ u˜i∂G∂xi= sT |∇G|, (2.76)where sT is the turbulent flame speed. While there exists several empirical corre-lations for sT [44] in the RANS context, this parameter is associated with someuncertainties in experimental measurements [132]. A detailed discussion of theG-equation model in the RANS context can be found in [128].In the LES context, a number of G-equation approaches have been proposed [69,83, 131]. Interestingly, the G-field, as originally proposed, exhibits a special sym-metry which makes it incompatible with LES [119]. Only recently has Pitsch [129]achieved consistency in the LES filtering of the G-field by the addition of a termwhich includes the effects of curvature on the flame propagation. Despite the sim-ple formulation, the G-equation approach has been used in a number of simula-tions, including LES of a premixed swirl-stabilised combustor [68], low-swirl pre-mixed flame [117] and other flames [69, 83, 88]. As with the problem found in theRANS formulation, closure for the turbulent flame speed is required; there doesnot appear to be a universal concensus on a suitable model or correlation for thisparameter [14].Artificially Thickened Flame ModelThe artificially thickened flame model, proposed by Butler and O’Rourke [27],was introduced to overcome the problem of unresolved flame fronts on coarse com-putational grids. The premise of the model, as the name implies, is to artifically in-crease the thickness of the flame while preserving the flame speed. This is achievedby manipulating the basic scaling laws from the theory of laminar premixed flames.In general, the laminar flame thickness, δL, can be related to the molecular diffu-382.10 Turbulent Premixed Combustion Modellingsivity, D, and the laminar flame speed, sL,δL ∝DsL. (2.77)Furthermore, the laminar flame speed is proportional to the molecular diffusivityand mean reaction rate, ω˙ , as follows,sL ∝√Dω˙. (2.78)By increasing the flame thickness by a factor of Q, from δL to Q ·δL, one must alsoincrease the diffusivity from D to Q ·D and decrease the mean reaction rate fromω˙ to ω˙/Q. An appropriate value of Q would be in the range of 5–30 to properlyresolve the artifically thickened flame front on typical LES grids. This model wasfirst applied to LES by Veynante and Poinsot [1] and has been developed furtherthereafter [35]. Since its introduction, the thickened flame model has been used tosimulate combustors with complex geometry [152], to study flashback and blow-off effects in a swirl burner [158], and to simulate the ignition sequence of a gasturbine engine [18].While the artificially thickened flame model demonstrates merit in accountingfor ignition and near-wall interaction, the factor, Q, is not complete and the physicsof the turbulence-chemistry interaction is not represented precisely. Indeed, an ad-ditional adjustable parameter, termed the efficiency factor, is required in the speciestransport equations to negate the changes introduced to the Damko¨hler number byQ.Presumed Conditional Moment (PCM) ModelThe Presumed Conditional Moment (PCM) Model was originally proposed byVervisch et al. [180] for non-premixed combustion in the RANS paradigm. Thismodel was later extended to LES as a SGS model by Domingo et al. for the sim-ulation of a premixed flame [39]. The model assumes a locally laminar structurefor a progress variable in the turbulent flames; the mean reaction rates are simply392.10 Turbulent Premixed Combustion Modellingcalculated by a convolution with the most appropriate PDF model,˜˙ωc = ∫ 10ω˙c|c∗P˜(c∗;~x, t)dc∗, (2.79)where ω˙c represents the source term for the selected progress variable. The pre-sumed PDF model can be the usual β -PDF for non-premixed combustion or themodified laminar flamelet PDF for premixed combustion [70, 153]. The ω˙c|c∗ istypically obtained from canonical, one-dimensional laminar flame calculations andtabulated prior to the LES or RANS simulations.One of the primary advantages of PCM is that it can be conveniently combinedwith tabulated chemistry reduction techniques, such as the Flame Prolongation ofILDM (FPI) model [55] described in Section 2.5.2. In this way, the PCM-FPI sim-ulations can emulate the effects of detailed chemistry at a fraction of the compu-tational cost of strategies requiring reaction mechanisms. The PCM-FPI approachhas demonstrated good agreement with experimental results in the simulations ofa lifted jet flame [40] and a lean premixed turbulent swirl flame [52]. Furthermore,this strategy has also been successfully applied to premixed Bunsen flames for alaboratory-scale burner [143, 145, 153]. The PCM-FPI appears to be a promisingtechnique for the simulations of premixed and non-premixed combustion, allowingfor detailed chemistry at a reasonable computational cost. Its main disadvantage isthat it is only formally valid in the flamelet regime, which could restrict its valuefor predictions of flames in industrial applications.Transported PDF ModelThe transported PDF model, proposed by Pope [134], solves a transport equa-tion for the joint PDF of the velocity and composition. Unlike the previous mod-els, this one does not assume a locally laminar structure for the turbulent flamesand is applicable to combustion problems outside of the flamelet regime. Thechemical source terms are perfectly closed at the expense of unclosed diffusionterms in the transport equation for the joint PDF. This tradeoff is of particularconcern for premixed combustion where the turbulent mixing and chemical kinet-ics are intricately linked. Nevertheless, the transported PDF model has been em-402.10 Turbulent Premixed Combustion Modellingployed in a number of RANS simulations for turbulent premixed laboratory-scaleflames [28, 101, 115, 163, 182].The transported PDF model has also been extended to LES by the introductionof a Filtered Density Function (FDF) [56, 135] and successfully implemented forthe LES of turbulent non-premixed [154] and premixed flames [189]. However,this model is not ideally suited for industrial problems as the Lagrangian MonteCarlo methods are required for the solution of the transport equation of the jointPDF. The additional computational complexity and cost hinder the transported PDFmodel from practical applications.Conditional Moment Closure (CMC)Conditional Moment Closure (CMC) was first proposed by Klimenko [85] andBilger [13] independently in the context of non-premixed combustion for RANS.CMC has a theoretical advantage over flamelet models in that it does not assumethe turbulent flame to be spatially thin. The RANS-CMC approach has demon-strated reliability in the chemical closure for non-premixed combustion [47, 132]and for diesel spray combustion in engines [124]. The successes led to the exten-sion of CMC to premixed combustion [84, 166]. More recently, CMC has beenapplied in the simulation of stoichiometric pilot-stabilised Bunsen flames [4]. Fur-thermore, the LES formulations of CMC have also been developed for both non-premixed [82, 116] and premixed [174] combustion, though the LES-CMC imple-mentation for industrial applications appears to be less frequent.The underlying hypothesis of CMC assumes that the conditionally-averagedchemical source terms (ω˙α |c∗) are closed by evaluating the reaction rate expres-sions with the conditionally-averaged scalars,ω˙α |c∗ ≈ ω˙α(T |c∗,Yα |c∗,ρ|c∗), (2.80)where c∗ is the discretised conditioning variable. T |c∗, Yα |c∗ and ρ|c∗ are the con-ditional temperature, mass fraction of species α and density, respectively. Theunconditionally-filtered chemical source term is subsequently calculated by in-tegrating the conditional source term with the PDF of the conditioning variable,412.10 Turbulent Premixed Combustion ModellingP˜(c∗;~x, t), ˜˙ωα = ∫ 10ω˙α |c∗P˜(c∗;~x, t)dc∗. (2.81)There are three main challenges in applying CMC to practical problems forturbulent premixed combustion. First, the conditionally-averaged scalars are cal-culated from transport equations that include the conditioning variable as an addi-tional independent variable; this added dimension translates to significantly highercomputational cost. Second, these transport equations contain a number of un-closed terms that are not well understood and require modelling [12]. As a result,the closure for these terms are typically not accurately represented. Finally, forpremixed combustion, the solution is extremely sensitive to the model for the con-ditional scalar dissipation. Models for this term, which permit CMC to be usedas a chemical closure strategy for premixed combustion, have only recently beendeveloped.Conditional Source-term Estimation (CSE)Conditional Source-term Estimation, first proposed by Bushe and Steiner [26],is a combustion model that provides closure for the mean chemical source termsfor turbulent flames. This model has been shown to produce promising results fordiffusion flames [67, 93, 161, 183]. Further, Salehi et al. has demonstrated thatCSE can also be applied to premixed flames and produce meaningful results [42,143–145, 153]. Recent effort has extended CSE to lifted turbulent flames [43]and flames operating in the Moderate and Intense Low Oxygen Dilution (MILD)combustion mode [92].CSE invokes the first-order Conditional Moment Closure hypothesis such thatthe conditionally-averaged chemical source terms are determined by evaluating thereaction rate expressions with the conditionally-averaged scalars (Equation 2.80).The unconditional mean chemical source term is calculated in the same manner, byintegrating the conditional source term with the PDF of the conditioning variable(Equation 2.81). The Favre-averaged presumed PDF models adopt the functionalform:P˜(c∗;~x, t)≈ P˜(c∗; c˜, c˜′2), (2.82)422.10 Turbulent Premixed Combustion Modellingwhere c˜ and c˜′2 are the Favre-averaged mean and variance of the progress variable.The usual spatial and temporal coordinates are represented by ~x and t. SeveralPDF models are tabulated in this format; a PDF model commonly implemented inconjunction with CSE is the modified laminar flamelet PDF [70, 143–145, 153].The central concept behind CSE is that there is more information available ifthe filtered field is considered in its entirety rather than as unrelated data points.Each filtered volume, governed by the same underlying physical processes, thencontributes additional information about the localised reactive flow. Specifically,the relevant conditional scalar averages are obtained by the inversion of the follow-ing integral equation [26],Y˜α(~x, t) =∫ 10P˜(c∗;~x, t)Yα |c∗(~x, t)dc∗, (2.83)where Y˜α(~x, t) is the Favre-averaged mass fraction. Figure 2.9 illustrates the inter-action between the CFD module, the CSE module, the chemistry module (FGM inthis example) and the PDF. The mass fractions, means and variances of the progressvariable are calculated by the CFD module via transport equations. These meansand variances of the progress variable are then used to obtain the PDFs, whichin conjunction with the mass fractions, serve as the input to the CSE module viaEquation 2.83. The computed conditional mass fractions are subsequently used todetermine the mean reaction rates through the FGM and the PDF model.OpenFoam CFDỸ CO2Ỹ H2Oω˙k |c*Y H2O |c*Y CO2 |c*CSẼ˙ωkc̃ c̃ '2PDFFGMFigure 2.9: CSE operational flowchart.43Chapter 3Pseudo-Turbulent ProbabilityDensity Function and ScalarDissipation Rate Models forPremixed Combustion3.1 IntroductionThere are several different models for turbulent premixed combustion that rely onhaving an accurate probability density function (PDF) of a reaction progress vari-able [39, 55, 70, 121]; these PDFs are conventionally expressed as a function ofthe mean and variance of that progress variable. Such presumed PDF approachesare used to account for finite rate chemistry, with tabulated data used to include de-tailed chemistry into simulations of turbulent combustion. Presumed PDF modelsare being developed for both the Reynolds-Averaged Navier-Stokes (RANS) andLES paradigms. However, the accuracy of these methods inevitably depends tosome considerable degree on the accuracy of the function presumed for the PDF ofthe progress variable.The oft-used β -PDF is able to recover the extreme properties expected of thetrue PDF, such as having δ s appear at 0 and 1 for maximal variance or a single δ443.1 Introductionappear at the mean for a zero variance; however, it fails to reproduce the shape ofthe true PDF [70]. Here, the issue is the shape of the PDF is strongly dependent onhow the chemical reaction rates vary as a function of the progress variable. Hence,different chemical kinetics will lead to a different shape of the PDF of progressvariable and the distribution of the β -PDF will, for most cases, lead to significantinaccuracies. Most critically, the inaccuracies can lead to a bias, as one typicallyfinds that the error is consistently always at the same values of progress variable inany particular flame.In an attempt to account for the effects of chemistry on the shape of the PDF,Bray et al. [22] proposed using a premixed laminar flame to model the functionaldependence of the PDF on progress variable. Their original formulation only pro-vided coverage for flames with very high variances. Jin et al. [70] proposed amodification of the Bray PDF that extends the original formulation to cover allpossible mean and variance combinations. This is accomplished by truncating thefunction as needed to match the required mean and variance. It was found to sig-nificantly improve the fit of the PDF to that extracted from DNS results. Howeverthere remains the issue that, at the point of truncation, the model PDF has a sharpdrop to zero; the true PDF tends to be more rounded.In this chapter, a one-dimensional pseudo-turbulent method is proposed to re-place the laminar flame that is conventionally used to generate PDF models for thereaction progress variable for turbulent reacting fields; the hope is that this strat-egy will smooth out these unphysical peaks in the distributions. Furthermore, forflames that are not in the flamelet regime, the pseudo-turbulent calculation mightcapture whatever effects small-scale eddies penetrating the reaction zone mighthave on the shape of the distribution.Hence, the objective of this study is to verify the ability of the Linear-EddyModel (LEM) to capture the overall PDF distributions for two test scenarios inthe premixed combustion mode. Particularly, PDFs generated by LEM simula-tions are compared to both experimental and DNS results. For the experimentalcomparison, LEM is used to simulate a premixed methane-air v-flame exposed tolow intensity turbulence. Temperature-based PDF models for the reaction progressvariable are subsequently constructed. In addition, the scalar dissipation rates arealso investigated for this v-flame. For the DNS comparison, LEM is used to sim-453.2 Probability Density Function Modelsulate a premixed methane-air flame in an isotropic, homogeneous turbulent field;the product-based PDF models are then constructed.3.2 Probability Density Function ModelsProduct-based and temperature-based progress variables can be formulated respec-tively as follows,c(x, t) =Yp(x, t)Ypeq, (3.1)c(x, t) =T (x, t)−TuTb−Tu , (3.2)where Ypeq , Tb and Tu are respectively the mass fraction of the products at completecombustion, burnt gas mixture temperature and unburnt gas temperature. In eithercase, one can formulate a model PDF with the following form,P(c∗;x, t)≈ P(c∗; c¯,c′2), (3.3)where c¯ and c′2 are the mean and variance of the progress variable. Such a form iscompatible with several current models for the turbulence-chemistry interactionsin turbulent premixed combustion.Four different PDF models will be considered here, three of which have beendiscussed previously: the β -PDF, the Bray PDF and a modified laminar flameletPDF. The β -PDF is typically a bimodal distribution widely used in non-premixedcombustion.The Bray, or laminar flamelet PDF, was originally formulated to model flamesin the thin flamelet regime. The Bray PDF is defined as [22]:P(c∗; c¯,c′2) = Aδ (c∗)+B f (c∗)+Cδ (1− c∗), (3.4)where f (c∗) is calculated using the solution of an unstrained laminar flame:f (c∗) =1|∇c∗| . (3.5)The coefficients A, B and C are obtained based on the first three moment equations463.2 Probability Density Function Modelsof c: ∫ 10P(c∗; c¯,c′2)dc∗ = 1 (3.6)∫ 10c∗P(c∗; c¯,c′2)dc∗ = c¯ (3.7)∫ 10c∗2P(c∗; c¯,c′2)dc∗ = c′2. (3.8)The Bray PDF is problematic in that it cannot produce PDFs for all values of meanand variance [143].Jin et al. [70] proposed a modification to this PDF whereby the function f (c∗)is truncated to zero over a range of c∗ and one or both of the δ s is eliminatedfrom the PDF in order to expand the applicable range of the PDF to all possiblevalues of mean and variance of c. In doing this, it was first recognized that thePDF will ultimately be stored in a look-up table for discrete intervals in c∗. Thatdiscretisation is used to provide bounds for using the laminar flame solution toassemble the PDF. A lower bound at cmin is defined which would be the maximumc∗ in the first discrete interval (that is to say that this discrete interval in c∗ would bebounded by 0 and cmin) and an upper bound at cmax is defined to be the minimumc∗ in the last discrete interval (this discrete interval in c∗ would be bounded by cmaxand 1).A region bounded by [x1,x2] — with corresponding values of the progress vari-able [c1,c2] — is located in the laminar flame such that the mean and variance of csampled within [x1,x2] match a given combination of c¯ and c′2; x1 and x2 then arecompared to xmin and xmax, which are the points in the laminar flame solution thatcorrespond to cmin and cmax. The functional form of the modified PDF changesdepending on the locations of x1 and x2 relative to xmin and xmax. Four cases arethen possible:• Case 1: xmin < x1 < x2 < xmaxP(c∗;~x, t) ={0 if c∗ < c1 or c∗ > c2B f (c∗) if c1 ≤ c∗ ≤ c2473.2 Probability Density Function Models• Case 2: x1 < xmin < x2 < xmaxP(c∗;~x, t) ={Aδ (c∗)+B f (c∗) if c∗ ≤ c20 if c∗ > c2• Case 3: xmin < x1 < xmax < x2P(c∗;~x, t) ={0 if c∗ < c1B f (c∗)+Cδ (1− c∗) if c∗ ≥ c1• Case 4: x1 < xmin < xmax < x2P(c∗;~x, t) = Aδ (c∗)+B f (c∗)+Cδ (1− c∗). (3.9)As in the Bray PDF, each of these cases has three unknown coefficients (among A,B, C, x1 and x2). The first three moment equations of c (Eqs. 3.6 through 3.8) areused to evaluate these. In practical implementation [144], the PDF is calculated in apre-processing operation where the laminar premixed flame is first computed, thenit is used with the above rules to fill in the complete range of mean and variance ofc and store the PDF in a table for later retrieval.Figure 3.1 shows PDFs collected from DNS data compared to the three differ-ent model PDFs described above. The β -PDF overestimates the probability over awide range of progress variable for c < 0.6, then utterly fails to capture the peakaround c = 0.85 and instead rises to a δ at c = 1. The original laminar flame PDFof Bray underestimates the PDF virtually over the entire range of progress vari-able and also includes a δ at c = 1 that is not exhibited by the DNS. The modifiedlaminar flamelet PDF clearly does the best job of representing the DNS PDF butit has too high a peak towards the product side of the distribution and exhibits asharp cut-off not evident in the DNS. To try to smooth off the sharp cut-off, a one-dimensional pseudo-turbulent flame calculation obtained using the Linear-EddyModel of Kerstein [77] is proposed to replace the usual laminar flame calculationfor the function f in the formulation above. The four PDF cases described for themodified laminar flamelet PDF will be the same, including the need to identifylocations in space, x1 and x2, at which to truncate the solution so as to fit the first483.3 Linear-Eddy Model(a) x/Ld = 0.234 (b) x/Ld = 0.391Figure 3.1: PDFs are exemplified at two axial locations of the DNS simula-tion domain (Ld = 6 mm). Solid: DNS; dash: modified laminar flamelet; dashdot: laminar flamelet; dot: β -PDFs [70].three moments of c.3.3 Linear-Eddy ModelThe Linear-Eddy Model is selected due to its minimal empiricism and its abil-ity to replicate turbulent flow conditions inexpensively [77–81]. It can be fun-damentally divided into two modules. The deterministic component consists ofone-dimensional gas dynamics evolutionary equations. The stochastic componentconsists of random eddy events. The two processes are implemented simultane-ously during the simulation to achieve a pseudo-turbulent effect. The couplingbetween the stochastic advective process and deterministic evolutionary equationson a one-dimensional computational domain permits the LEM to simulate turbulentflows with higher Reynolds numbers than multi-dimensional models. The currentLEM variant is adapted from the original formulation to accommodate turbulentpremixed reacting flows.3.3.1 Evolutionary EquationsThe following set of evolutionary equations is derived for an unsteady, isobaric,quasi-one-dimensional flame [76]. Small changes in pressure across the premixed493.3 Linear-Eddy Modelflame reaction zone permit the assumption of an isobaric system. The continuity,energy, species transport equations and equation of state are respectively,M˙ = ρuxAc, (3.10)ρAc∂T∂ t+ M˙∂T∂x− 1cp∂∂xΛAc∂T∂x+AccpN∑α=1ρcpαYαVα∂T∂x+AccpN∑α=1ω˙αhαMα = 0,(3.11)ρAc∂Yα∂ t+ M˙∂Yα∂x+∂∂x(ρAcYαVα)−Acω˙αMα = 0, (3.12)ρ =pMRuT. (3.13)The M˙, ux, Ac and Vα respectively represent the mass flow rate, streamwise velocity(assumed to be in the x direction), the cross-sectional area of the streamtube1 andthe diffusion velocity of species α . The other symbols are previously defined inSection 2.2.3.3.2 Stochastic AdvectionLEM postulates a random process that rearranges fluid elements along a line in or-der simulate the chaotic vortices that appear in turbulent flows. The inversion gen-erates discontinuous fluid motions, which lead to random walks of fluid elements.Each inversion requires three random variables: epoch, position and length. Thesevariables are determined from the turbulent Reynolds number (Ret) governing theflow via Kolmogorov scaling [77]. For the purposes of the LEM simulations to beperformed here, the turbulent Reynolds number should correspond to the turbulentReynolds number of the flame to be modeled. The eddy inversions, illustrated byFigure 3.2, are now referred to as triplet maps.1 Ac is assumed to be unity for all of the LEM simulations in this thesis.503.4 LEM Simulation ParametersFigure 3.2: Triplet map applied to a line. Cell values are transported withinthe eddy interval; interpolation is not required at any point [77]. The before(left) and after (right) domain values are illustrated.This mapping strategy is able to conserve the scalar fields because the contentof each control volume within the selected interval is only rearranged to a neigh-bouring position. Interpolation is not required at any point. In this way, inadvertentchanges in the solution vectors due to continuous implementation of the stochasticadvection will not occur.3.4 LEM Simulation Parameters3.4.1 LEM Methodology for Experimental StudyA reduced, six-step mechanism designed by Chen et al. [29] to simulate premixedmethane-air combustion is implemented to replicate the experimental flame. Thepropagation speed and overall temperature profile for the unstrained laminar flamehas been calibrated to the reference solution from Cantera [58] by modifying thereaction rates of the mechanism. This mechanism is selected because it is able toreproduce the temperature profile sufficiently well compared to the full chemistry;other reduced mechanisms with fewer reaction steps tend to be less satisfactoryin this regard. For the purposes of this chapter, the temperature gradient must befairly accurate as the functional form of the PDF model depends on the normalizedinverse of this parameter.All of the species properties are calculated through CHEMKIN-II [74], in-cluding specific heats, diffusion coefficients, thermal conductivities and enthalpies.513.4 LEM Simulation ParametersMixture-averaged transport is adopted to reduce computational time. Figure 3.3 il-lustrates the flame solutions from Cantera and the six-step mechanism. The overallflame thickness, propagation speed, equilibrium temperature and mass fractions ofmajor species are well matched.Figure 3.3: Premixed laminar flame solution at an equivalence ratio of 0.73for the six-step reduced mechanism (line) and Cantera [58] (symbols).The inflow mixture is set to atmospheric pressure and 294 K at an equivalenceratio of 0.73. The integral scale and turbulent Reynolds number are respectively1.9 mm and 38. The LEM simulations are performed with a mixed first orderupwind and second order centered spatial scheme. Explicit timesteps are taken inorder to cope with the discontinuities introduced by the stochastic component ofthe model. Each realisation is a snapshot of the instantaneous temperature profilerecorded at 1/100th large eddy turnover time intervals. 10,000 realisations are usedto generate the PDFs and scalar dissipation curves at each of the specified turbulentReynolds numbers.523.5 Construction of PDF Models3.4.2 LEM Methodology for DNS StudyA simplified, two-step, five-species mechanism is used to simulate a methane-airpremixed flame in a one-dimensional, homogeneous turbulent field. The two reac-tions are chosen to be as economical as possible while retaining crucial characteris-tics of premixed methane-air flames. Further information regarding the mechanismand inlet conditions can be found in [62]. The integral scale and turbulent Reynoldsnumber are respectively 4.5 mm and 51. All species have equal and constant heatcapacities of 1250 J/kg·K. The binary diffusion coefficients and thermal conduc-tivities are equal for all species and vary only with temperature. The mechanismand thermochemical parameters are selected to emulate the DNS flame as closelyas possible.Similar numerical methods prescribed in the previous subsection are applied;however, each realisation is a snapshot of the instantaneous product mass fractionrecorded at 1/10th large eddy turnover time intervals. 10,000 realisations are aver-aged to generate the final PDF models.3.5 Construction of PDF ModelsThe method used to construct the LEM PDF models will first be discussed in thecontext of the modified laminar flamelet PDF then extended to the LEM formula-tion. As per Section 3.2, the relevant parameters governing the distributions are themean (c¯) and variance (c′2) of the PDFs. Truncations (in space) applied to the one-dimensional flame profile would lead to changes in both c¯ and c′2. The mean of thedistribution will increase for every point removed from the unburnt mixture bound-ary while the contrary is true for points removed from the burnt mixture boundary.The variance can be adjusted in a similar manner. Figure 3.4(a) shows an exampleof such truncations. The symbols correspond to the truncation limits, (x1,c1) and(x2,c2), where only the cells within the interval are retained. The arbitrary valuesof mean and normalised variance for illustrative purposes are set to c¯ = 0.5 andc′2n = 0.33; the normalised variance is determined by the limiting variance of a533.5 Construction of PDF Modelsperfectly segregated mixture of hot and cold gases:c′2n =c′2c¯(1− c¯) . (3.14)The modified laminar flamelet PDFs with desired means and variances can thenbe constructed with the truncated flame profiles according to methods outlined inSection 3.2. An example is shown in Figure 3.4(b). Crucially, it can be deducedthat there exists one unique modified laminar flamelet PDF for every mean andvariance combination.(a) Truncated laminar flame profile (b) Modified laminar flamelet PDFFigure 3.4: (a) (x1,c1) and (x2,c2) mark the truncation limits where only thecells within the interval are retained; (b) Modified laminar flamelet PDF withc¯ = 0.5 and c′2n = 0.33 is constructed using the truncated flame profile.The LEM PDFs are constructed analogously; however, the LEM flame profilesare not steady with time. Figure 3.5 illustrates five sample realisations obtainedfrom LEM calculations where each profile is slightly different but equally valid.Thus, it is necessary to average the contributions from each flame realisation in or-der to arrive at the final model. The truncation strategy prescribed for the modifiedlaminar flamelet profiles to obtain the desired mean and variance is applied to eachLEM result.Figure 3.6(a) illustrates two truncated LEM profiles with c¯ = 0.5 and c′2n =0.33, similar to the laminar flame shown in Figure 3.4(a). The truncation positionsare shifted for individual realisations because the distributions are inherently dif-543.5 Construction of PDF ModelsFigure 3.5: Individual flame realisations from LEM calculations.ferent. The final PDF model is acquired by averaging numerous PDFs built fromsuch truncated LEM flame profiles. Figure 3.6(b) illustrates an example of theconverged solution constructed by averaging 10,000 LEM PDFs at a relatively lowturbulent Reynolds number. It is apparent the LEM PDF displays peaks of lowermagnitude at each boundary, thereby redistributing the probability of finding c.This smoothing effect is typical of PDFs constructed with LEM.553.5 Construction of PDF Models(a) Truncated LEM flame profiles (b) Modified laminar flamelet and LEMPDFsFigure 3.6: (a) The portions to be retained are within the intervals marked by’o’. The truncation positions are selected such that the resultant PDF wouldhave c¯= 0.5 and c′2n = 0.33. The truncation boundaries are different for eachrealisation because the profiles change with time. (b) The modified laminarflamelet (dash-circle) and LEM (solid) PDFs of similar mean and variance areshown.3.5.1 ConvergenceThe number of LEM profiles required to construct a statistically converged PDF so-lution appears to depend on two criteria; the first of these is the turbulent Reynoldsnumber governing the stochastic advection. A flow of higher turbulent intensityrequires additional LEM profiles for a smooth probability distribution. The secondcriterion is the desired resolution in the progress variable (c) space. The numberof LEM profiles needed to construct a converged PDF will increase with the reso-lution; typically 5000 realisations is sufficient for Ret < 200 combined with 0.02resolution in c-space. This behaviour is illustrated by Figure 3.7. It can be seen thatthe PDFs constructed from 2000 and 5000 profiles are nearly identical at Ret = 50.563.6 Scalar Dissipation Rate(a) c-space resolution = 0.02 (b) c-space resolution = 0.02 (magnified toboxed region)Figure 3.7: LEM generated PDFs at Ret = 50. Each line represents a differentnumber of LEM profiles used for the construction of the PDF. Solid: 5000profiles; dash: 2000 profiles; dash-dot: 1000 profiles.The computational time required to fully construct a LEM PDF model for pre-mixed combustion is directly proportional to the number of LEM profiles neededfor a converged solution. A typical quad-core computer can produce on the orderof thousands of flame profiles per 24-hour period, depending on the complexityof the chemistry and the large eddy turnover time of the system. However, thecomputational time can be reduced linearly with the number of processors as theLEM simulation on each processor runs independently. This allows for quick andstraightforward implementation on large computing clusters. If sufficient proces-sors are available, one can effectively amass the flame profiles on the order ofhours. After collecting the required number of realisations, the PDF model can bebuilt in a matter of minutes.3.6 Scalar Dissipation RateFor premixed combustion, the Scalar Dissipation Rate (SDR), χc, of a temperature-based reaction progress variable is defined as,χc(T,φ) = αc(T,φ)|∇c||∇c|, (3.15)573.6 Scalar Dissipation Ratewhere αc(T,φ) and ∇c are the thermal diffusivity at the local temperature andequivalence ratio and the gradient of the progress variable, respectively. The scalardissipation rate is an important quantity to both non-premixed and premixed com-bustion modeling [45, 181]. It is an unclosed term that appears in the transportequation for the variance of progress variable; χc directly measures the decayingspeed of fluctuations through turbulent micromixing [181]. Since the burn rate ofmany combustion processes depends on the contact area and local gradient betweenreactants, it is inevitable that the mean burning rate of the flame either explicitlyand implicitly depends on the scalar dissipation in most combustion models. Forexample, Conditional Moment Closure (CMC) uses scalar dissipation rate condi-tioned by the progress variable to calculate the micromixing term [181]. Coinci-dentally, modeling χc|c∗ emerges as one of the main difficulties in applying CMCto turbulent premixed combustion problems [4].A relationship can be established between the conditional scalar dissipationrate (χc|c∗) and the PDF of the reaction progress variable (P(c∗)) [181],χc =∫ 10(χc|c∗) ·P(c∗)dc∗. (3.16)Equation 3.16 is particularly interesting to the analysis because LEM realisationscould be used to construct both χc|c∗ and P(c∗). The methods applied to obtainthe PDF model has been described in the previous section; the conditional scalardissipation rates can be retrieved similarly. First, the thermal diffusivity at the lo-cal temperature and equivalence ratio, αc(T,φ) is generated from laminar flamesimulations by Cantera [58] and stored prior to calculations. Second, the progressvariable gradients (∇c) at the corresponding progress variables (c) are obtainedfrom LEM flame realisations. This is then combined with the stored values of αcaccording to Equation 3.15 at each discrete interval of c to form the conditionalscalar dissipation rate. Third, the temperature gradients are arranged into 50 binsof size 0.02, varying between 0 to 1 in c-space. The values within each bin are av-eraged and additional LEM realisations are introduced until statistical convergenceis achieved. It is worth noting that, for a given chemistry and set of LEM parame-ters, there would be one pseudo-invariant conditional average of scalar dissipationrate predicted. Practically, it appears that the addition of an entry for the uncon-583.7 Results and Discussionditional mean of scalar dissipation in the PDF look-up table that would containthat pseudo-invariant conditional scalar dissipation rate already combined with thecompanion PDF using Equation 3.16 would be appropriate; this would be done forevery combination of mean and variance of progress variable stored in the tablefor easy look-up when needing the closure for this term in the progress variablevariance transport equation.3.7 Results and Discussion3.7.1 Probability Density FunctionsLEM – Experiment ComparisonThe premixed methane-air experiment was performed by a multi-slot burner devel-oped at the University of Cambridge. The equivalence ratio (φ ), inlet temperatureand pressure were respectively 0.73, 294 K and one atmosphere. The flame struc-ture, temperature and relevant species concentrations were measured via a com-bination of line imaging of spontaneous Raman scattering, Rayleigh scattering,two-photon laser-induced fluorescence (LIF) of CO and crossed planar imaging ofOH-LIF. Further information on the Cambridge slot burner and multiscale laserdiagnostics used for the experiment could be found in [8].PDF models are built with LEM using the temperature-based reaction progressvariable defined by Equation 3.2 at three preconditioned means. Experimental,modified laminar flamelet and LEM PDFs are shown in Figure 3.8. For a precon-ditioned mean of c¯ = 0.05, it is evident from Figure 3.8(a) that the LEM formu-lated PDF captures both the position and magnitude of the experimental result ina superior manner than the modified laminar flamelet PDF towards the left bound-ary. The modified laminar flamelet PDF well under-predicted the distribution for0.02 < c∗ < 0.2 and over-predicted the result close to c∗ = 0. However, it is ableto capture the slight increase in probability near c∗ = 0.9 whereas LEM cannot.The small compromise in the LEM PDF within this vicinity is compensated by thesignificant improvement for small c values.The more curious observation is shown by Figure 3.8(b), where both the mod-593.7 Results and Discussionified laminar and LEM PDFs differ from the experimental result. This could becaused by at least two different possible sources of error. First, the underlying re-action properties from the reduced mechanism may be incorrect; it may not be ableto fully capture the reaction rates, thus causing the temperature gradients whichgovern the PDF to be inaccurate. Second, the temperature gradients for the experi-ment are measured 25-30 mm downstream of the slot burner [168]. This could po-tentially skew the distribution such that low c values are inadequately representedcomparing to high c values. Typically one would expect to find slight to moder-ate peaks towards the boundaries of the PDF (for c¯ ≈ 0.5) because reactants andproducts tend to dominate the total mass fraction (excluding inert species) at anygiven time. This means there should be a higher probability to find the reactioneither near the beginning or the finish. The modified laminar flamelet and LEMPDFs both predict this result, opposite of the experimental data. One could specu-late measurements at an axial location closer to the slot exit may reveal a differentPDF even at the same preconditioned mean. Nevertheless, it can be seen that theLEM PDF achieves better overall resemblance to the experimental PDF than theprevious formulation.603.7 Results and Discussion(a) c¯ = 0.05 (b) c¯ = 0.53(c) c¯ = 0.95Figure 3.8: Experimental PDFs conditioned at three different means of thereaction progress variable. Solid: experiment [167]; dot: modified laminarflamelet; dash: LEM. The vertical axis represents the probability at state c∗conditioned by c¯, P(c)|c¯.For a preconditioned mean of c¯ = 0.95, it can be seen from Figure 3.8(c) thatLEM closely captures the experimental distribution in peak magnitude and overallshape. The modified laminar flamelet PDF tends to mildly over-predict the prob-abilities for 0.7 < c∗ < 0.9 and is not able to acquire either the peak location orthe magnitude for c∗ > 0.9. Considering the three preconditioned means, it can beconcluded that LEM PDFs tend to be more comparable to the experimental resultsthan the modified laminar flamelet PDFs.613.7 Results and DiscussionLEM – DNS ComparisonThe product-based PDF model constructed with LEM is superimposed with theDNS and modified laminar flamelet PDFs in Figure 3.9. It can be seen that theLEM formulated PDF is able to capture the gradual decay of probabilities near0.85< c∗ < 0.95 exhibited by the DNS result; the peak values near complete com-bustion are more closely approximated than the modified laminar flamelet PDF.LEM is able to achieve this smoothing effect because of the stochastic nature ofthe model. Each realisation is slightly different; when averaged, they lead to amore realistic distribution of flame speeds and temperature gradients. Midrangevalues of c∗ reveal excellent agreement for all three PDF models. However, it isevident that LEM under-predicts the probability for 0.02< c∗ < 0.2. This could bean indication of over-mixing in the preheat layer, leading to a higher rate of productformation. Such a process inadvertently increases the gradient of the progress vari-able, thus lowering the probability of finding it at those states. The colder regionof the preheat layer is particularly susceptible to this problem because the flamegenerally does not have enough time to recover before successive triple maps. Ingeneral, it can be seen that LEM can acquire the positions, distributions and peakvalues of the DNS PDF in a superior manner than the modified laminar flameletmethod.Figure 3.9: PDF at x/Ld = 0.391 of the DNS simulation domain (Ld = 6mm). Solid: DNS; dash: modified laminar flamelet; dash dot: LEM. TheDNS and modified laminar flamelet PDFs are obtained from [70].623.7 Results and Discussion3.7.2 Scalar Dissipation RateEach scalar dissipation result is constructed using 10,000 LEM realisations. Fig-ure 3.10 illustrates the scalar dissipation rate conditioned in c-space from experi-mental data [167] and LEM at two turbulent Reynolds numbers.Figure 3.10: Scalar dissipation conditioned in c-space. Solid: experi-ment [167]; dot: laminar flame; dash: LEM at Ret = 38; dash-dot: LEMat Ret = 200.It is apparent that the model is able to reflect the decrease in peak amplitudesof the distributions relative to the laminar flamelet result. This amplitude decreaseis likely attributed to thickening of the flame, also seen in previous DNS studies[147]. Another scalar dissipation rate generated from a higher turbulent Reynoldsnumber is also investigated; the tendency for the peak temperature gradients todecrease with Ret coincides with current understanding of flame broadening asturbulent intensity increases. This investigation also reveals that LEM may not besufficiently mixing the flow field as Ret = 200 still does not appear to lower thetemperature gradients to the level of the experimental result.Moreover, it may indicate the reduced mechanism is reacting too rapidly, caus-ing an increased temperature gradient. Though the peak amplitudes are generallyhigher than the experimental values for 0.5 < c∗ < 0.8, it serves to validate thegeneral physics of turbulent flames is appropriately represented by the stochasticadvection introduced by LEM. The model has demonstrated the ability to acquirethe qualitative behaviour of the temperature gradients and scalar dissipation rates633.8 Concluding Remarkssatisfactorily even with a six-step, reduced mechanism.3.8 Concluding RemarksOne objective of this chapter has been to examine an alternate PDF model for pre-mixed combustion which uses the LEM paradigm to mimic the effect of turbulentmixing on the PDF. LEM has the advantages of low computational cost and mini-mal empiricism. PDFs generated by LEM have demonstrated the ability to quali-tatively match DNS and experimental results with greater success than previouslysuggested PDF models. LEM has also been used to predict the decrease in con-ditional scalar dissipation rate caused by flame thickening due to turbulent effectsin premixed combustion, albeit over-predicting the experimental results in peakmagnitudes. The stochastic nature of LEM enables turbulent reactions in a simplegeometry to be captured inexpensively and qualitatively. It can be concluded thatLEM is suitable for data pre-processing of simulation strategies that require the useof PDFs and scalar dissipation rates.64Chapter 4Direct Comparison of PDF andSDR between LEM Simulationsand Experiments for Turbulent,Premixed Methane-Air FlamesThe work in this chapter is an extension of the LEM study in Chapter 3. Therange of turbulent conditions simulated by the LEM has been increased and thecomparison with experimental results is more complete. In addition, a number ofimprovements to the software has been introducted, including the parallelisationof the code that constructs the probability density function (PDF) and scalar dis-sipation rate (SDR) models. As this chapter is an extension of the previous, someintroductory material and acronyms are inevitably repeated with the thought ofimproving readability.4.1 IntroductionSeveral practical models for turbulent premixed combustion rely on an accuraterepresentation of the probability density function of a reaction progress variable,which is often parameterised by the mean and variance of that progress variable [39,55, 70, 121, 143, 145]. These presumed PDF approaches are often implemented in654.1 Introductionconjunction with tabulated chemical variables to achieve detailed chemistry calcu-lations in turbulent combustion simulations. Such models have been developed forboth the Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy Simulation(LES) paradigms. Previous work has shown that the accuracy of these methodsdepends to a considerable extent on the accuracy of the function presumed for thePDF of the progress variable [70].A number of different presumed PDF models have been previously investigatedfor premixed combustion. The often used β -PDF does recover the extreme prop-erties expected of the true PDF, such as δ functions at the zero and unity extremesof reaction progress for maximal variance, and single δ functions at the mean forzero variance. However, it fails to reproduce the shape of the true PDF in moregeneral cases [70]. The issue is related to the fact that the shape of the true PDFappears to be a function of how the chemical reaction rates vary as a function of theprogress variable; hence, different chemical kinetics lead to different shapes of thePDF of progress variable. The form of the β -PDF — which is of course entirelyindependent of the chemical kinetics — can lead to significant inaccuracies. Mostcritically, there can be a biasing error, as the discrepancies tend to occur at the samevalues of the progress variable for any particular flame.One of the primary concerns in the design of modern engines is the reduction ofharmful pollutants; specifically, the current generation of numerical models mustbe able to predict both the thermodynamic properties of the reacting mixture andthe formation of these minor species with sufficient accuracy to resolve the partsper million produced. A prominent example is the prediction of prompt flame NOxvia the Fenimore pathway [49]. The predicted NOx values from this mechanismare strongly affected by the choice of the PDF model, as this pathway is sensitiveto the predicted temperatures and flame profile in the reactive regions.In an attempt to account for the effects of chemistry on the shape of the PDF,Bray et al. [22] proposed using a premixed laminar flame to model the functionaldependence of the PDF on progress variable. The proposed probability dependenceof the flame existing at any given state is inversely proportional to the magnitude ofthe gradient of the temperature. Their original formulation only provided coveragefor flames with very high variance. Jin et al. [70] then proposed a modification tothe Bray PDF that extends the original formulation to cover all possible mean and664.1 Introductionvariance combinations. This is accomplished by truncating the PDF shape functionas needed to match the mean and variance parameters. It was found to significantlyimprove the fit of the PDF to that extracted from Direct Numerical Simulation(DNS) results. A shortcoming of this method is that, at the point of truncation, themodel PDF has a sharp drop to zero, whilst the true PDF tends to be more rounded.To address this issue, a one-dimensional turbulent method was proposed to takethe place of the typical laminar flame calculation to tabulate pseudo-turbulent PDFmodels for RANS and LES closures. The Linear-Eddy Model (LEM), an inex-pensive one-dimensional stochastic mixing model, has demonstrated the ability tocapture important effects from the interaction between chemistry and turbulenceon the PDF distributions sufficiently well [148, 157, 175]. This model providesus with a mechanism to investigate the general flame characteristics at very highturbulence intensities, much beyond the capability of current DNS strategies. Inturn, it permits us to analyze the behavior of the PDF constructed at these highlyturbulent states.While the LEM has been implemented to investigate the shape of the PDFdistributions [148, 157], such simulations were performed to analyze the PDF forspecific flames. The current study is primarily interested in the tabulation of aPDF lookup table useful for subsequent RANS and LES flame computations (apre-processing operation not unlike pre-calculating the β -function and storing thatin a lookup table). There was an initial suggestion that one ought to try to matchthe turbulence statistics in the LEM calculations to those that one expects to findin the later turbulent flame calculation. Indeed, a possibility is that one mightneed to add a dimension to the lookup table (something like the local turbulentReynolds number) to account for variations in local turbulence properties in theturbulent flame calculation and their effect on the shape of the PDF. A large partof the motivation for the LEM work presented here is to answer the followingquestions: Do the turbulence properties affect the shape of the PDF? If so, how?How important is it to match the LEM turbulence properties used in generating thePDF lookup table to those that will be found in the turbulent flow to be calculatedlater?A related question can be asked about the local gradient of progress variablein a premixed turbulent flame, which is closely related to the SDR. For premixed674.1 Introductioncombustion, the SDR, χc, of a temperature-based reaction progress variable is,χc(T,φ) = αc(T,φ)|∇c| · |∇c|, (4.1)where αc(T,φ) and ∇c are the thermal diffusivity at the local temperature andequivalence ratio and the gradient of the progress variable, respectively.The SDR is an important quantity to both non-premixed and premixed com-bustion modeling [45, 181]. It is an unclosed term that appears in the transportequation for the variance of progress variable, where χc directly measures the de-cay rate of fluctuations through turbulent micromixing [181]. Since the burn rateof many combustion processes depends on the contact area and local gradient be-tween the reactants, it is reasonable for most combustion models to assume thatthe mean burning rate of the flames either explicitly or implicitly depends on thescalar dissipation rate. For example, Conditional Moment Closure (CMC) uses thescalar dissipation rate conditioned by the progress variable to calculate micromix-ing [181]; not surprisingly, modeling the conditional scalar dissipation term χc|c∗conditioned on the local and global progress of reaction emerges as one of the maindifficulties in applying CMC to turbulent premixed flames [4].The same LEM method can be used to generate a model for the conditionalscalar dissipation, where one simply conditionally averages the scalar dissipation inthe LEM temperature profiles that are used to construct the PDF model. The uncon-ditional mean SDR can be obtained by taking the inner product of the conditionalSDR with the model PDF. This allows for the construction of pseudo-turbulentPDF and SDR models that are perfectly consistent with one another. While it ispossible to obtain PDF and SDR distributions for premixed combustion from DNS,the associated cost is generally prohibitive for flows with relevant turbulent condi-tions [70]. More importantly, studies typically tend to focus on the analysis of thePDFs and SDRs at specific points within the domain. This leads to the problemthat although DNS and experiments can provide valuable insight to the behaviorof the PDF and SDR [168, 170, 193], they cannot provide usable input models forsubsequent RANS and LES combustion calculations, which is the primary motiva-tion behind the current work. The incorporation of turbulence characteristics in thePDF distributions could provide a solution to the deficiencies seen in a number of684.2 Numerical Conditions: Premixed Combustioncurrent PDF models, such as the ad hoc β -pdf or the laminar approaches. Havingmore accurate PDF and SDR models would be beneficial to a number of RANSand LES strategies as closures for turbulent premixed combustion typically rely onsome variant of the PDF or SDR model for the reaction progress variable as theinput [45, 181].Recent detailed measurements of species and temperature have been made bythe Cambridge-Sandia swirl burner [169–171, 194]. This swirl burner was de-signed specifically to explore the influence of stratification on the flame. However,the very detailed nature of the scalar and velocity measurements have made the dataset attractive as a target for premixed flame model validation as well [72, 137]. Inparticular, the comprehensive database allows conditioning on a number of differ-ent variables, including equivalence ratio (for the stratified flames), temperature (orprogress of reaction) or any other suitable scalar.In this chapter, the experimental dataset is used to obtain detailed PDFs ofthe progress of reaction. The temperature is used to characterise the extent ofreaction, for direct comparison with the PDFs generated from LEM simulations forthree different swirl (and turbulence) levels. The measured and computed variancesare considered, as well as the detailed shape of the PDFs, in the comparison. Inaddition, the unconditional mean SDR values obtained from experiment and LEMsimulations are also directly evaluated. In the following sections, the numericalapproach will first be discussed, followed by a summary of the experiments andthe data treatment used.4.2 Numerical Conditions: Premixed Combustion4.2.1 Linear-Eddy ModelThe Linear-Eddy Model has been demonstrated to replicate the flow statistics forsimple turbulent conditions with acceptable accuracy [77–81]. Given the one-dimensional nature of the model, the computational costs remain relatively lowfor most practical cases. Here, the LEM is used in a pre-processing manner forthe tabulation of discrete PDF and SDR models, which can be implemented insubsequent RANS and LES applications.694.2 Numerical Conditions: Premixed CombustionThe LEM can be divided into two modules. The deterministic component con-sists of the usual one-dimensional gas dynamics evolution equations, whereas thestochastic component consists of random eddy events. The turbulence conceptof LEM postulates a random process that rearranges fluid elements along a linein order simulate the chaotic vortices that appear in turbulent fields. These one-dimensional vortices, known as triplet maps, generate discontinuous fluid motions,which lead to a random walk of fluid elements. Three quantities control the behav-ior of each eddy event: the frequency, the size and the position. Firstly, the eddyevent frequency per unit length of the domain, λLEM, is governed by [79],λLEM =545νRetCλ l30(l0/ηK)5/3−11− (ηK/l0)4/3, (4.2)where ν , Ret , l0 and ηK are the kinematic viscosity, turbulent Reynolds number,the integral length scale and the Kolmogorov scale, respectively. The empiricalparameter Cλ should typically be tuned to the flame in question for LEM stud-ies [148]; however, as the primary interest here is to investigate the changes in thePDF models associated with variations in the turbulent fluctuations and integrallength scales, Cλ is held constant for consistency between test cases. The value of15.0 for Cλ is adopted from [148]. Similarly, the parameter used to scale the Kol-mogorov length, Nη , which typically requires tuning to case specific conditions,is also held at a constant value of 1.0 for all of the LEM cases. This parameter isapplied to the inertial scaling law in the following form,ηK = Nη l0Re−3/4t . (4.3)Secondly, the size of the turbulent eddies are sampled from a PDF1,f (l) =53l−8/3η−5/3K − l−5/30, (4.4)that varies within the interval of ηK < l < l0, derived from inertial scaling [79].Thirdly, the position of the eddies are determined from a random uniform distribu-1 For clarity, it should be explicitly mentioned that this PDF of eddy sizes is different than thePDF of the reaction progress variable.704.2 Numerical Conditions: Premixed Combustiontion within the domain.The deterministic and stochastic modules are implemented simultaneously dur-ing the simulation to achieve a pseudo-turbulent effect; this coupling betweenthe stochastic advective process and deterministic evolution equations in a one-dimensional computational domain permits the LEM to simulate turbulent flowswith higher Reynolds numbers than multi-dimensional models. The current LEMvariant follows the formulation of [120] to accommodate premixed turbulent react-ing flows.To a large extent, the applicability of LEM depends on the homogeneity ofthe turbulence and the degree of which the problem can be characterised in onedimension within the simulation domain. This suggests that the method may failwhere there is significant coupling between turbulence non-homogeneity at thesmall scales and the flow field. For many problems in premixed combustion, thespatial scales where reactions take place tend to be of the order of the laminarflame thickness, and in many cases, experiments show that at such small scalesof turbulence, its characteristics have decayed to conditions that are reasonablyisotropic and uniform.4.2.2 LEM Simulation MethodsA global, six-step mechanism by Chang et al. [29] designed to simulate premixedmethane-air combustion is implemented to generate the flame profiles. The sixglobal reaction rates have been reduced by approximately 8% to achieve the cor-rect unstrained laminar flame propagation speed observed in the reference solutionfrom Cantera’s GRI MECH 3.0 [156] calculation. This extra calibration procedureensures increased accuracy for this value of equivalence ratio of 0.732.All of the species properties are calculated using CHEMKIN-II [74], includingspecific heats, diffusion coefficients, thermal conductivities and enthalpies. Thethermodynamic coefficients are based on the CHEMKIN Thermodynamic DataBase [75]. Figure 4.1 illustrates the flame solutions from Cantera and the Chang2Reduced mechanisms are typically designed to operate over a range of thermodynamic condi-tions. However, due to a more restricted degree of freedom from the reduced number of steps, thelevel of accuracy may not be identical over the range of operations. For this numerical study, theprimary interest is with premixed flames at φ = 0.73; thus, the reduced mechanism is optimised forthis particular equivalence ratio.714.2 Numerical Conditions: Premixed Combustion00.020.040.060.080.10.12MassFractions(Major)00.511.522.53MassFractions·10−3(Minor)0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1020406080100120140Velocity(cm/s)x (cm)0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 120040060080010001200140016001800Temperature(K)x (cm)00.150.30.450.60.750.91.051.21.35Density·10−3(g/cm3)CH4H2OCO2COH2HOHODensity TemperatureFigure 4.1: Premixed laminar flame solution at an equivalence ratio of 0.73for the six-step global mechanism (line) and Cantera [58] (symbols).mechanism. Mixture-averaged transport is adapted to reduce computational time.The inflow mixture was set to atmospheric pressure and 294 K at an equivalenceratio of 0.73. The calculated laminar flame speed, sL, under these conditions is0.214 m/s, and the calculated laminar flame thickness, δ f , is 588 µm. The overallflame thickness, propagation speed, equilibrium temperature and mass fractions ofmajor and minor species are sufficiently well-matched between the reduced Changand full GRI mechanisms.The LEM simulations are performed with a mixed first order upwind and sec-ond order centered spatial scheme. Explicit time steps are taken in order to copewith the stochastic nature of the model. The instantaneous temperature profiles arerecorded at regular time intervals on the order of 1% of a large eddy turnover timefor the freely propagating flames. Some 3,500 to 16,000 temperature profiles arestored for post-simulation construction of the LEM PDF models. The exact num-ber of profiles required to create a statistically converged PDF model is dependentprimarily on the turbulence intensity of the flame. Ten prototype turbulent pre-mixed flames with various integral length scales and turbulent Reynolds numbersare tested. The flow parameters are varied from case to case such that the prototypeflames adhere to predetermined locations on the Borghi diagram, as illustrated in724.2 Numerical Conditions: Premixed Combustion10−1 100 101 102 10310−1100101102103u′/sLl0/δflaminarflamesbrokenreactionzones4b3a 3b 3c2a 2b 2c1a 1cRe= 1lk = 0.01δflk = δf1bFigure 4.2: Borghi diagram showing locations of the ten prototype LEMflames (1a to 4b). The LEM test cases are represented by ’◦’. The exper-imental flames are represented by triangles: SwB1 (’C’), SwB2 (’5’), andSwB3 (’B’). It has been found that increasing the integral length beyond anorder of magnitude above the laminar flame thickness while holding the tur-bulent fluctuations constant does not significantly alter the PDF profiles forthe LEM simulations.Figure 4.2. The LEM parameters used to simulate the freely propagating flames atthe prescribed turbulent conditions are summarised in Table 4.1.Throughout the chapter, a temperature-based reaction progress variable is used:c =T −TuTeq(φ)−Tu , (4.5)where Teq(φ) is the equilibrium temperature at the equivalence ratio φ and Tu isthe unburnt gas temperature. From the simulations or experiments, it is possible to734.2 Numerical Conditions: Premixed CombustionCase Ret l0 (cm) ηK (cm) u′/sL Domain (cm) Min. cells/cm Cells used1a 5 0.02992 8.95·10−3 9.826 1.0 670 1,0001b 50 0.2992 1.59·10−2 9.826 1.0 337 1,0001c 500 2.992 2.83·10−2 9.826 4.0 318 1,2722a 50 0.02992 1.59·10−3 98.26 1.0 3,771 3,7722b 500 0.2992 2.83·10−3 98.26 1.0 2,121 2,1222c 5000 2.992 5.03·10−3 98.26 4.0 1,193 4,7723a 500 0.02992 2.83·10−4 982.6 1.0 21,204 21,2043b 5000 0.2992 5.03·10−4 982.6 1.0 11,924 11,9243c 50000 2.992 8.95·10−4 982.6 4.0 6,706 26,8244b 50000 0.2992 8.95·10−5 9826 1.0 67,053 67,054Table 4.1: Relevant LEM simulation parameters: l0 and ηK are the integraland Kolmogorov scales. Other constant parameters are invariant between thecases, including δ f , sL, Cλ and Nη , which are respectively, 588 µm, 0.214 m/s,15.0 and 1.0. A minimum of 1,000 computational cells are used for eachsimulation.construct a PDF with the following form,P(c∗;x, t)≈ P(c∗; c¯,c′2), (4.6)where c¯ and c′2 are the mean and variance of the progress variable, and c∗ is thediscretized variable representing the continuous c-space. In another words, thePDF at any space and time within the simulation domain can be approximated bythe mean and variance of the distribution. Such a formulation is compatible withseveral current models for the turbulence and chemistry interactions for turbulentpremixed combustion [39, 55, 70, 121]. A detailed description of the LEM PDFformulation can be found in [175], while a brief overview of the construction ofthe LEM PDF is provided in Section 4.4. Furthermore, a normalised variance isdefined, c′2n, determined by the limiting variance of a perfectly segregated mixtureof hot and cold gases as:c′2n =c′2c¯(1− c¯) . (4.7)This formulation is more suitable for the discrete tabulation of the PDF table, asthe boundary values conveniently vary between zero and unity.744.3 Experimental Conditions: Stratified Swirl Burner4.3 Experimental Conditions: Stratified Swirl BurnerThe flames and experimental techniques are summarised briefly in this section3as they have been described in previous papers [170–172]. A turbulent flame isstabilised on a bluff body, with reactants fed through two concentric streams atspecified inner and outer equivalence ratios (φ ). In the present chapter, only pre-mixed cases are considered. Swirl can be added to the flame by splitting the outerstream through a tangential inlet.The global equivalence ratio for the flame is nominally 0.75 (measured viaspecies is 0.73 +/- 0.018), and the range of local equivalence ratios spans 0.375-1.125. Scalar data obtained from Rayleigh/Raman/CO-LIF line measurements at103 µm resolution allows the behaviour of key combustion species — CH4, CO2,CO, H2, H2O and O2 — to be probed within the instantaneous flame front. Simul-taneous cross-planar OH-PLIF is used to determine the orientation of the instanta-neous flame normal in the scalar measurement window, allowing real gradients ofthe temperature and hence the progress variable to be obtained.The operating conditions considered in the present analysis are listed in Ta-ble 4.2. Measurements were taken at six heights along the centerline of the flame,by collecting data in 6 mm linear segments. The notation SwBNz is used to denotecases at different locations, where N is the case number and z is the axial coordi-nate. Three premixed flames under different swirl conditions are considered.Case SFR(1) SN(2) u′/s(3)LSwB1 0 0 10.2SwB2 0.25 0.26 12.2SwB3 0.33 0.45 14.3Table 4.2: Operating Conditions: (1) SFR = ratio of split flow to swirlersto total flow, (2) SN = measured swirl number, ratio of tangential to axialmomentum, (3) Maximum total u′/sL at the midpoint of the flame brush at z= 30 mm [71, 194].The bulk velocity in the outer annulus, Uo = 18.7 m/s, is set at more than3It is important to explicitly mention that the experimental expertise and analysis are provided byMr. M. Kamal and Professor S. Hochgreb of the University of Cambridge. They are the co-authorsof the work documented in this chapter.754.3 Experimental Conditions: Stratified Swirl Burnertwice the value of the velocity in the inner annulus, Ui = 8.3 m/s, so as to generatesubstantial levels of shear, and thus turbulence, between the two flows. Co-flow airwas supplied around the outer annulus with a bulk velocity Uco = 0.4 m/s to providewell-characterised boundary conditions. The Reynolds numbers derived from thebulk velocities at the exit geometry are Rei = 5,960 for the inner flow and Reo =11,500 for the outer flow.Multi-scalar laser diagnostics were applied at the Turbulent Combustion Lab-oratory in Sandia National Laboratories, and extensively described in [170–172].The diagnostics setup allows for the line measurement of temperature (Rayleighscattering) and major species (Raman scattering and CO-LIF) at 103 µm projectedpixel resolution with simultaneous cross planar OH-PLIF at 48 µm projected pixelresolution. Signal to noise ratios are of the order of 150 for temperatures, and about60 for equivalence ratio, with estimated accuracies of 2% and 5%, respectively. Asthe optical resolution of the Raman-Rayleigh-LIF measurements is smaller than thespatial sampling rate, the resolution of the temperature and major species measure-ments is limited by the sampling resolution (103 µm) and the laser beam diameter(0.22 mm, 1/e2). The optical resolution (1/e2) of the OH-PLIF measurement isbetween 98 and 144 µm and therefore the spatial resolution of the OH-PLIF mea-surements is limited by the optical resolution rather than the sampling resolution.The OH-PLIF cross planar technique allows the flame normal in 3D to be assessedrelatively to the measurement line, thus allowing the real gradients to be obtainedby projection. Radial profiles were obtained by moving the burner horizontallyin 4 mm steps, producing overlapping steps in the relative position of the 6 mmwide measurement window, with 300 laser shots taken at each step. Radial profileswere taken at axial increments of 10 mm above the burner exit to capture changesin flame structure with axial distance. Substantially larger datasets (5,000-30,000shots) were taken at the intersection of the mean flame brush and the mixing layerfor stratified flames only. Here, only the premixed cases confined to within 2.5%of the nominal value of φ = 0.75 are considered.764.4 Construction of PDF Models4.4 Construction of PDF Models4.4.1 LEM Flame ProfilesFigure 4.3 illustrates a few selected LEM temperature profiles for two prototypeflames of varying turbulence intensities. It can been seen that the influence of thetriplet maps becomes more pronounced and modify a larger portion of the flame asthe turbulent Reynolds number increases. Consequently, the turbulent temperatureprofiles become more distinguishable from the laminar counterpart, with largereffective flame thicknesses. These perturbations propagate through the flame andreduce in strength until they either blend into the flame profile or get transportedout of the domain.774.4 Construction of PDF Models0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1200400600800100012001400160018002000Temperature(K)x (cm)(a) Case 2b (Ret = 500).0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1200400600800100012001400160018002000Temperature(K)x (cm)(b) Case 4b (Ret = 50,000).Figure 4.3: Characteristic LEM temperature profiles of two prototype flamesat different Ret . The individual profiles on each graph are separated by at leastone large eddy turnover time.At first inspection, it appears that the average temperature gradients shouldincrease with the number of triplet maps implemented per unit time; however, thediffusion mechanism from the one-dimensional evolution equations quickly dimin-ishes the sharp discontinuities introduced in the slope of the temperature field. As784.4 Construction of PDF Modelsa consequence, the flame begins to broaden as it recovers from this perturbed state.The net effect is to decrease the overall conditional averages of the gradients whenconsidered across the entire width of the flame. The temperature profiles depictedin Figure 4.3 can then be transformed into the progress variable space via Equa-tion 4.5.4.4.2 LEM PDF ConstructionThe method used to construct the LEM PDF models is first discussed in the contextof the modified laminar flamelet PDF, which is then extended to the LEM formu-lation. The relevant parameters governing the behavior of the distributions are themean (c¯) and variance (c′2), as previously mentioned by Eq. 4.6. Truncations ofthe profiles in the spatial domain are applied to the one-dimensional flame profiles,leading to changes in both c¯ and c′2. The mean of the distribution increases forevery point removed from the unburnt mixture (c = 0) boundary, whist the oppositeis true for points removed from the burnt mixture (c = 1) boundary. Moreover, thevariance decreases for every point removed from either boundary.Figure 4.4(a) demonstrates an example of these truncations. The symbols cor-respond to the truncation limits, (x1,c1) and (x2,c2), where only the cells withinthe interval are retained. Arbitrary values of mean and variance are set to c¯ =0.50 and c′2n = 0.39 for illustrative purposes. The modified laminar flamelet PDFswith the prescribed means and variances can then be constructed with the truncatedflame profiles according to methods described by previous work [70]. An exam-ple of the modifield laminar flamelet PDF of c¯ = 0.50 and c′2n = 0.39 is shown inFigure 4.4(b). Crucially, it can be deduced that there exists one unique modifiedlaminar flamelet PDF for every mean and variance combination.794.4 Construction of PDF Models0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91cx (cm)(x1, c1)(x2, c2)(a) Truncated laminar flame profile0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.3P(c∗)c∗(b) Modified laminar flamelet PDFFigure 4.4: (a) (x1,c1) and (x2,c2) mark the truncation limits; only the cellswithin the interval are retained. (b) The modified laminar flamelet PDF with c¯= 0.50 and c′2n = 0.39 constructed using the truncated flame profile is shown.The LEM PDFs are constructed in a similar manner compared to the modifiedlaminar flamelet PDFs. Instead of one unique, steady, laminar temperature profile,the LEM generates transient temperature profiles; however, each LEM temperatureprofile is equally valid and must be equally weighted in an ensemble average toarrive at the final PDF model. Consequently, the truncation strategy prescribed forthe modified laminar flamelet profiles to obtain the desired mean and variance isapplied to each LEM temperature profile. Here, when a (x1, x2) pair is fixed in realspace, averaging together the many LEM realizations leads to a particular meanand variance of progress variable. In principle then, one could vary the positionsof x1 and x2 for each temperature profile until the desired mean and variance forthe PDF is achieved. In practice, all possible discrete (x1, x2) combinations areused to populate the PDF lookup table, which is stored as a function of the meanand variance. During subsequent RANS and LES calculations, this table is calledand interpolated as needed to obtain the PDF for any required mean and variancecombination.Figure 4.5(a) illustrates two truncated LEM profiles that correspond to PDFswith c¯ = 0.5 and c′2n = 0.33, similar to the laminar flame shown in Figure 4.4(a).The truncation positions are shifted for individual temperature profiles because the804.4 Construction of PDF Modelsflame changes with time. The final PDF model is acquired by averaging a num-ber of PDFs built from such truncated LEM flames. Figure 4.5(b) illustrates anexample of the converged solution constructed by averaging 10,000 LEM PDFs ata relatively low turbulent Reynolds number. It is apparent the LEM PDF displayspeaks of lower magnitude at each boundary than the modified laminar flamelet PDFmodel. This smoothing effect is a direct consequence of the application of tripletmaps to the temperature profiles. The mapping introduces variations in the temper-ature gradients at a given temperature (Figure 4.3(b)); as a result, the probabilityof the flame existing at any given state varies with different temperature profiles.When averaged, the PDF decreases to zero at each boundary more gradually [175].0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91cx (cm)(a) Truncated instantaneous LEM flame pro-files0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.3P(c∗)c∗(b) Modified laminar flamelet and LEMPDFsFigure 4.5: (a) The portions to be retained are within the intervals delimitedby the circles. The four truncation positions are selected such that the resul-tant PDF would have c¯ = 0.5 and c′2n = 0.39. The truncation boundaries aredifferent for each temperature profile because of the transient effects. (b) Themodified laminar flamelet (dash) and LEM (solid) PDFs of similar mean andvariance are shown.814.5 Results4.5 Results4.5.1 Probability Density FunctionThe LEM PDF models spanning a range of distribution means and variances areillustrated by Figure 4.6. This study reveals that the characteristics of the PDF mod-els for all ten prototype flames display remarkable similarities across the scope oftest scenarios despite a three-order change of magnitude in turbulent fluctuation in-tensity and length scale on the Borghi diagram. The PDFs tend to become slightlyless bivariate as the turbulence intensities increase. This effect is particularly no-ticeable for cases with means and normalised variances having values close to 0.5.The PDF distributions captured at various axial positions for the three exper-imental flames (SwB1-3) operating at different swirling conditions are illustratedby Figures 4.7-4.9. The LEM PDFs tabulated as a function of c∗ from one par-ticular prototype flame has been superimposed for direct comparison between themodel and the swirl burner data4. In general, a high degree of similarity can be ob-served between the experimental and LEM PDFs in overall features, magnitudesand peak positions. The 3 × 5 panels correspond to different distances from thestabilisation point of the flame and radial positions across the flame brush, whichare characterised by various combinations of c¯ and c′2n.Starting with SwB1 (Figure 4.7), the model PDFs are well-matched for allnormalised variance values below a distribution mean of approximately 0.7. Thereactant edge of the flame brush can be found at low c¯ (left columns); and thePDF peaks at c∗ = 0. This agrees with the calculated LEM PDFs, but the measuredPDFs are somewhat wider than the calculations suggest. At the intermediate valuesof c¯, two low peaks are found, one at zero and one around 0.8-0.9, which are well-captured by the LEM model, even at low variance levels. Closer to the productend of the flame brush (right columns), the PDF peak appears not at unity, but atslightly lower values, between 0.8 and 0.9, much like the experimental values. Thepeak at lower values than full reaction progress is attributed to the limited time forcompletion of reaction due to turbulent mixing of cold and hot gases. However,4LEM prototype flame 1b is selected for this demonstration. However, it is evident from Fig-ure 4.6 that the location on the Borghi diagram does not significantly influence the resulting LEMPDF for a given variance and mean combination.824.5 Resultsheat transfer to the base of the flame can also contribute to the lowering of theextent of reaction.00.10.20.3P(c∗)|c¯=0.11c¯′2n= 0.25 c¯′2n= 0.49 c¯′2n= 0.7500.10.20.3P(c∗)|c¯=0.3100.10.20.3P(c∗)|c¯=0.5100.10.20.3P(c∗)|c¯=0.710 0.2 0.4 0.6 0.8 000.10.20.3P(c∗)|c¯=0.910.2 0.4 0.6 0.8 0c∗0.2 0.4 0.6 0.8 1Figure 4.6: PDF models at various c¯ and c′2n are shown; each row representsone value of c¯ and three values of c′2n (ranging from left to right: 0.25, 0.49and 0.75). Vertical and horizontal axis on each graph represent the probabilityand the progress variable, respectively. Solid: case 4b, dash: case 3b, dashdot: case 2b, dot: case 1b (notation is in accordance with Figure 4.2).834.5 Results00.10.20.3P(c∗)c′2n = 0.618c ≈ 0.11c′2n = 0.664c ≈ 0.31c′2n = 0.652c ≈ 0.51c′2n = 0.609c ≈ 0.71z=60mmc′2n = 0.479c ≈ 0.9100.10.20.3P(c∗)c′2n = 0.560 c′2n = 0.586 c′2n = 0.562 c′2n = 0.509z=30mmc′2n = 0.2680 0.2 0.4 0.6 0.8 000.10.20.3P(c∗)c′2n = 0.2700.2 0.4 0.6 0.8 0c′2n = 0.3450.2 0.4 0.6 0.8 0c∗c′2n = 0.3460.2 0.4 0.6 0.8 0c′2n = 0.2870.2 0.4 0.6 0.8 1z=10mmc′2n = 0.024Figure 4.7: PDFs measured from the SwB1 flame at various axial locationsconditioned by different distribution means, c¯ (solid). The correspondingLEM PDFs of similar c¯ and c′2n are superimposed (dash).Proceeding to the first swirl flame (Figure 4.8), it can be seen that the LEM andexperimental PDFs share as many similarities as with the no swirl case. The modeland experiment agree fairly well for most values of variances at a distribution meanof less than approximately 0.7. The experimental PDFs tend to be slightly lessbivariate in comparison to the no swirl case, reflected by the decreasing distributionvariance. An interesting difference can be observed for the results having c¯ ≈0.7 and c′2n < 0.4. The probability of c∗ continues to be non-zero towards the 0boundary for the experimental PDFs whereas the LEM PDFs tend to go to zero.This could suggest that the swirling flow may prolong the time during which theflame spends in the preheat layer, effectively increasing the probability of findingthe flame at a low c∗ state and decreasing the gradient of the progress variablearound this region of the flame.844.5 Results00.10.20.3P(c∗)c′2n = 0.489c ≈ 0.11c′2n = 0.623c ≈ 0.31c′2n = 0.613c ≈ 0.51c′2n = 0.562c ≈ 0.71z=60mmc′2n = 0.449c ≈ 0.9100.10.20.3P(c∗)c′2n = 0.573 c′2n = 0.593 c′2n = 0.563 c′2n = 0.509z=30mmc′2n = 0.2730 0.2 0.4 0.6 0.8 000.10.20.3P(c∗)c′2n = 0.3900.2 0.4 0.6 0.8 0c′2n = 0.4350.2 0.4 0.6 0.8 0c∗c′2n = 0.4120.2 0.4 0.6 0.8 0c′2n = 0.3410.2 0.4 0.6 0.8 1z=10mmc′2n = 0.112Figure 4.8: PDFs measured from the SwB2 flame at various axial locationsconditioned by different distribution means, c¯ (solid). The correspondingLEM PDFs of similar c¯ and c′2n are superimposed (dash).For the SwB3 flame (Figure 4.9), the LEM and experimental PDFs are typicallywell-matched for distributions with 0.3 < c¯ < 0.7. For PDFs with low values of c¯(c¯ ≈ 0.11), it can be seen that the experimental PDF peaks have lower magnitudesand in some cases, are shifted towards a non-zero c∗ value (c∗ = 0.05). More so, thedistributions appear to be slightly wider than the LEM counterparts. This suggeststhat increasing the swirl number of the flame may have greater impact on the colderregions of the flame, where the conditional mean is low — such a behavior is notapparent in the SwB1 and SwB2 cases. Overall, the PDFs for all cases, with andwithout swirl, are still captured for the corresponding combination of means andvariances. However, the agreement is generally better for the intermediate valuesof c¯ than for the extremes.854.5 Results00.10.20.3P(c∗)c′2n = 0.394c ≈ 0.11c′2n = 0.500c ≈ 0.31c′2n = 0.578c ≈ 0.51c′2n = 0.581c ≈ 0.71z=60mmc′2n = 0.367c ≈ 0.9100.10.20.3P(c∗)c′2n = 0.496 c′2n = 0.567 c′2n = 0.585 c′2n = 0.556z=30mmc′2n = 0.3540 0.2 0.4 0.6 0.8 000.10.20.3P(c∗)c′2n = 0.4500.2 0.4 0.6 0.8 0c′2n = 0.4980.2 0.4 0.6 0.8 0c∗c′2n = 0.4890.2 0.4 0.6 0.8 0c′2n = 0.4430.2 0.4 0.6 0.8 1z=10mmc′2n = 0.317Figure 4.9: PDFs measured from the SwB3 flame at various axial locations,conditioned by different distribution means, c¯ (solid). The correspondingLEM PDFs of similar c¯ and c′2n are superimposed (dash).4.5.2 Scalar Dissipation RateA relationship can be established between the conditional scalar dissipation rate(χc|c∗) and the PDF of the reaction progress variable (P(c∗)) [4],χc =∫ 10(χc|c∗)P(c∗) dc∗, (4.8)where χc is the unconditional mean scalar dissipation rate, or simply the meanscalar dissipation rate. Equation 4.8 is particularly interesting to the analysis be-cause the LEM temperature profiles can be used to construct both χc|c∗ and P(c∗).A more detailed description of the methodology used in constructing the PDF andSDR models can be found in Chapter 3.There are two observable behaviors from the conditionally averaged SDR mod-els from the LEM. First, as the characteristic flame properties are shifted upwards864.5 Resultsin the Borghi diagram (increasing u′/sL), the peak magnitudes of the dissipationrates tend to decrease. Second, if the characteristic flame properties are shiftedrightward in the Borghi diagram (increasing l0/δ f ), the peak magnitudes of theSDRs tend to increase towards the limit of a laminar flame. It is difficult to deter-mine from the current simulation results whether the flame broadening (increasingu′/sL) or the integral scale increment (increasing l0/δ f ) effect is dominant withina specific combustion regime. It is clear, however, that the processes will competewith one another in general, leading to relatively unchanged conditionally averagedSDR distributions for the special case of isotropic, homogeneous flames.For a given chemical mechanism and set of LEM parameters, there should beone pseudo-invariant conditional scalar dissipation rate (cSDR) averaged from thetemperature profile datasets. This cSDR can then be multiplied with the compan-ion PDF through an inner product to arrive at the unconditional mean SDR viaEquation 4.8. To reduce simulation time during practical implementation for com-plex combustion models, the values of the unconditional mean SDR can be pre-tabulated at desired combinations of means and variances of the reaction progressvariable for efficient retrievals. The same LEM datasets from the PDF study canbe used to generate one conditionally averaged SDR model for each of the tenprototype flames. The results are presented in Figure 4.10.Figure 4.11 illustrates the pseudo-invariant conditional SDR distributions fromthe experiments. The conditional SDR typically decreases with increasing swirland appears to be most similar to the laminar limit at z = 30 mm. This axial locationroughly coincides with the intersection of the mixing and shear layers, promotingthe development of a highly turbulent and dissipative region, after which the SDRquickly diminish. The comparison with the LEM model is striking: in the model,there is little sensitivity to the local turbulence characteristics, whereas the experi-ments indicate that there are significant changes in the SDR values; the magnitudetends to be more sensitive to the effects of swirl. Nevertheless, the shape of thedistributions are well-captured, with the maximum in the vicinity of c∗ = 0.65 to0.75. The maximum generally moves rightward, towards c∗ = 0.75, as the SDRmagnitude decreases.874.5 ResultsFigure 4.10: Non-dimensionalized conditional average of the SDR modelsfrom the ten prototype flames. The left and right columns illustrate changesin the turbulent fluctuations and integral length scales, respectively. The tur-bulent Reynolds number of each case is recorded in the parentheses. Thelaminar case is displayed in light gray as a reference.As mentioned in the Introduction, several models for turbulent combustion re-late the local rate of micromixing to the conditional scalar dissipation. A commonapproach to modelling scalar dissipation [113, 164, 165] is to extract the functionaldependence on the conditioning variable and model the conditional scalar dissipa-tion as,χ|c∗= χ0× f (c∗), (4.9)where χ0 is then calculated by some other means and is not a function of the condi-tioning variable. As usual, the integration of f (c∗) over the progress variable spaceresults in the value of unity. Here, the same decomposition is applied to the resultsto see if the underlying functional dependence of scalar dissipation on the condi-tioning variable changes — that is, to test whether or not such a decompositionmight be sensible.884.5 ResultsFigure 4.11: Non-dimensionalized conditional average of the SDR from thethree experimental flames. Each column represents a unique axial positionfrom the flame stabilisation point while each row represents one swirling con-dition. The error bars indicate +/- one standard deviation from the mean. Thelaminar case is displayed in light gray as a reference.Figure 4.12 depicts the functional dependence of scalar dissipation ( f (c∗)) atthree axial positions from the flame stabilisation point. This figure confirms theobservation that the cSDRs tend to shift rightward with increasing swirl intensity,particularly at the beginning and end of the flame brush, corresponding to the ax-ial positions of 10 mm and 60 mm. The reason for such changes in the shape ofthe normalised cSDR is not immediately clear from the current results; rather thanspeculate on the possible causes, this question can be an opportunity for futurework. In general, the LEM results tend to remain close to that of the laminar flameas the model is unable to emulate the effects of swirl. Furthermore, Figure 4.13illustrates the values of the normalization constant, χ0, in relation to the axial po-sition for the experimental flames. It is interesting that the maximum values of χ0for all of the experimental flames occur at the intersection of the mixing and shearlayers near z = 30 mm.894.5 Results0 0.2 0.4 0.6 0.8 000.511.522.53c∗χ|c∗/χ010 mm0.2 0.4 0.6 0.8 0c∗30 mm0.2 0.4 0.6 0.8 1c∗60 mmFigure 4.12: Functional dependence of the scalar dissipation ( f (c∗)) at threeaxial positions downstream from the flame stabilisation point (10, 30 and60 mm) with the LEM results superimposed. The experimental flames arerepresented by symbols: SwB1 (’+’), SwB2 (’’), and SwB3 (’◦’). TheLEM result is represented by the solid line. The laminar case is displayed inlight gray as a reference.10 20 30 40 50 6000.20.40.60.81z (mm)χ0=∫1 0χc|c∗dc∗×δ f/SLFigure 4.13: Values of χ0 in relation to the distance downstream of theflame stabilisation point. The experimental flames are represented by sym-bols: SwB1 (’+’), SwB2 (’’), and SwB3 (’◦’). The LEM (black line) andlaminar flame (gray line) results are not dependent on the axial position.The inner product of the pseudo-invariant conditional SDR with the appropri-ate PDF distribution via Equation 4.8 provides the solution for the unconditionalmean SDR, χc, or simply the mean SDR. The modelling of this term remains asone of the final challenges in applying some combustion models to premixed flamecalculations [181]. Figure 4.14(a) illustrates the behavior of the mean SDR as a904.6 Discussionfunction of the mean and variance of the progress variable as calcuated by LEM.The mean SDR peaks around c¯ = 0.5 and c′2n ≈ 0 as the PDF and cSDR distri-butions have coinciding maxima for these values of c¯ and c′2n. The experimentalvalues can be seen in Figure 4.14(b), where the magnitudes of the mean SDR areevidently lower than the predicted values. The LEM predictions are best when theswirl intensity of the burner is at its lowest. This suggests that swirling flows in-duce a decrease in the magnitude of the SDR, both conditional and unconditional,which cannot be replicated by a one-dimensional turbulence model. Moreover, re-calling the PDF results from the previous section, one can conclude that swirl playsa more prominent role in the determination of the SDR than the PDF distributions.4.6 Discussion4.6.1 Probability Density FunctionThe LEM formulated PDF models display minimal variations in distribution shapeand magnitude with changing turbulent fluctuations and integral lengths. Typically,the PDF distributions become slightly less bivariate with increasing turbulent fluc-tuations. Towards the reactant boundary (c∗ = 0), the turbulent fluctuations appearto increase the diffusive mixing in the preheat layer, leading to slightly faster flamedevelopment in this region. In turn, this causes a decrease in the duration of theflame residing at low values of c∗. Towards the products boundary (c∗ = 1), themodel suggests that the eddy interactions impede the flame from fully oxidizingtowards the end of the reaction. This does not imply a reduction in the flame speed,but rather a barrier to reaching full chemical equilibrium. This effect can be ob-served in Figure 4.6, where the PDF models constructed from flames with higherturbulent fluctuations have the rightward boundaries shifted slightly towards theleft. The combination of these phenomena causes the LEM PDF to deviate fromthe bivariate distribution seen in laminar PDF models as turbulence intensity in-creases, in agreement with the experimental results. The model of course neglectsany potential heat transfer issues, which are present in the experiments, and mayinfluence the approach to equilibrium at the base of the flame.The swirl burner results (Figures 4.7-4.9) reveal that the shape of the PDFs914.6 Discussion(a) LEM00.250.50.751c¯′2 n=0.300.250.50.75c¯′2 n=0.4χ¯c00.250.50.75c¯′2 n=0.50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.250.50.75c¯c¯′2 n=0.6(b) Experiments and LEMFigure 4.14: Non-dimensionalized unconditional mean scalar dissipationrates (χ¯c× δ f /sL) as predicted by (a) LEM and (b) experiments with LEMsolutions superimposed at various combinations of c¯ and c′2n. Solid: LEM;dash: SwB1; dash-dot: SwB2; dot: SwB3.924.6 Discussionis dependent on the distance from the base of the flame and secondarily on theamount of swirl in the flame. The dependence on distance from the base arisesfrom the different turbulence characteristics across the flame brush as the flowevolves downstream. The variance increases with distance from the base, and withswirl, and this directly affects the corresponding PDFs. The experimental PDFsare perhaps surprisingly well-captured by the two input parameters of c¯ and c′2n.Nevertheless, once the variance and mean are determined, there are only so manysensible ways for the PDF to be arranged within the LEM context, so perhaps thisjustifies the good agreement.The data suggests that variations in the turbulent Reynolds number and integrallengths induce a rather modest effect on the shape of the PDF; it appears that it maynot be necessary to add an additional dimension related to the turbulence intensityto the lookup table. Indeed, while it would be ideal to construct a table using repre-sentative conditions, using a model generated from a flame with lower turbulenceintensities (which is the less computationally intensive option) may suffice. Thenumber of temperature profiles required to construct converged LEM PDF mod-els can be markedly reduced while the additional error contribution remains small.This is especially true for premixed flames with turbulent Reynolds numbers andintegral lengths within an order of magnitude of the targeted flame on the Borghidiagram. Further discretion may be required when applying different PDF modelsfor flames exposed to more intense swirling conditions. The experimental resultsreveal that the overall flame gradients and hence the flame surface density functionreduce in magnitude under swirl. However, such changes may not greatly alter theshape of the PDF of the progress variable because of the normalization procedurerequired in the construction of the LEM PDF models.4.6.2 Scalar Dissipation RateAs previously mentioned, there are two observable behaviors from the condition-ally averaged SDRs computed by the LEM. First, as the intensity of turbulence in-creases, the peak magnitudes of the mean dissipation rates tend to decrease. Suchan effect has been previously observed in DNS studies and is attributed to flamebroadening caused by the decrease in peak temperature gradients [147]. This ef-934.6 Discussionfect is emulated by the LEM via the interactions between the turbulent eddies andthe scalar fields. These small eddies penetrate into the reaction zone, increasingthe local gradients of the scalar fields. This is followed by the dissipation of thesharp discontinuities introduced in the slope of the scalar fields due to the coupleddiffusion mechanism. Effectively, the flame is broadened, increasing its thicknessand decreasing the gradients. Second, as the turbulent length scale increases atconstant turbulent fluctuation intensity, the peak magnitudes of the SDRs tend inincrease towards the limit of a laminar flame5. Such an effect could be connectedto the increasing Kolmogorov length associated with the integral scale increment.An increase in the Kolmogorov length means the smallest eddies will become toolarge to effectively change the local structure of the flame. In another words, theeddies have more difficulty penetrating the preheat layer and the reaction zone,thus reducing the interaction between the vortices and the flame. With a smallernumber of interactions, the pseudo-turbulent SDR distributions will become moresimilar to the laminar SDR distribution.It is important to emphasise that the current numerical results do not includethe effect of strain. Additional modelling parameters would need to be introducedto the transport equations for such calculations. Previous studies have shownthat regions under highly positive strain are correlated with regions of increasedscalar dissipation, though this is not an exclusive association [90, 178]. For tur-bulent flames exposed to high strain rates, neglecting this effect may lead to non-negligible changes to both the PDF and SDR distributions.The pseudo-invariant conditional SDR distributions and mean SDR values fromthe experimental flames show a more pronouced decrease in magnitude with in-creasing swirl effects. It appears that the intersection of the mixing and shearlayers near z = 30 mm leads to quickly diminishing gradients of c for positionsdownstream of this location. This suggests that additional parameters which quan-tify the dissipative effects of mixing and swirling may be necessary in order toaccurately model the conditional SDR distributions. Moreover, the entrainmentof cold air into the premixed flame caused by the increase in fluid motion may5 For unstrained flows, the maximum conditionally averaged SDR distribution corresponds to thelaminar flame with the same chemistry; any turbulence interactions decrease the magnitude of theSDR distribution with respect to the laminar limit.944.7 Conclusioninduce changes in the temperature based reaction progress variable that the LEMcannot currently simulate. It is, however, curious that the PDF distributions fromboth the LEM and experiments reveal far less susceptibility to changes when ex-posed to physical effects, such as entrainment of air and swirl, in comparison to theconditional SDR distributions.Nevertheless, the quantity that requires closure in combustion models is theunconditional mean SDR, defined by Equation 4.8. In this regard, it appears thatthe LEM can reproduce values of χc well within 50% of the experimental equiva-lent, but the accuracy depends primarily on the swirl number of the flow. For mostpoints, the LEM over-estimates the value of χc, which is understandable consider-ing the conditional SDR distributions from the model indicate higher magnitudescompared to the experiments. Further research could explore the possibility of in-cluding additional modelling parameters into the LEM formulated conditional SDRmodel according to physical effects governing the flow field. In this case, perhapsa post-processing correction can account for the entrainment of air induced by theswirl. This could amount to applying a simple correction factor to the magnitudeof the cSDR, which is calculated as a function of the local anisotropy.4.7 ConclusionThis chapter investigates the effects of variation in swirl and turbulence intensity onthe probability distribution of the reaction progress variable for a series of globallylean, turbulent premixed flames, both experimentally and numerically.The PDF models constructed from the LEM simulations indicate that increas-ing turbulence intensity has a rather modest impact on the distributions. The minorchanges are only observed at medium values (approximately 0.5) of means andnormalised variances. The PDFs from the experimental flames show a bivariatedistribution for all cases, with peaks in the unburned and burned regions. The ther-mal flame thickness is generally larger than the unstrained (and strained) laminarflame thickness. The data therefore suggests that the thermal gradient is smearedby turbulent diffusion at scales on the order of the flame thickness. The LEMPDF model demonstrates good agreement with the experimental results in termsof nominal changes with turbulence. Particularly, it is able to capture the overall954.7 Conclusionshape and the effect of smoothing towards the left and right boundaries of the PDFdistributions.The LEM results suggest that the conditional SDR would decrease towardsthe upper regions of the Borghi diagram and would increase in peak magnitudestowards the unstrained laminar limit with increasing integral scales. The experi-mental conditional SDR for all axial locations of the flame brush shows a greaterdecrease in magnitude than predicted by the numerical model, though this couldbe partially caused by an interaction between shearing and mixing layers withina specific region of the swirl flames, leading to greater dissipation downstream ofthis location. The values for the unconditional mean SDR predicted by LEM aretypically within 50% of the experimental values. Moreover, the accuracy of thenumerical model seems to be dependent on the swirl number of the flow.Overall, the results point towards the idea that the PDF for premixed flames re-mains relatively steady under a variety of turbulent conditions, including changesin the integral length, velocity fluctuation and swirl number. As a consequence,it appears to be practical to use a representative pseudo-turbulent PDF model fora range of turbulent conditions. Future studies could focus on understanding thephysical interpretation behind the relative invariance in the PDF model with varia-tions in both the turbulent fluctuations and integral length scales. Another possibledirection is to determine if this invariance also extends to other fuels in the pre-mixed combustion regime. Swirl appears to have greater influence on the localflame gradients; this, in turn, leads to more substantial changes in the SDR, whichthe model presented here cannot capture in its current form. Perhaps a practicalsolution in the near future would be to implement a correction factor to the LEMformulated SDR models based on the local anisotropy of the flow field.96Chapter 5Conditional Source-termEstimation: Parallel IterativeSolution with DynamicEnsembles5.1 IntroductionTurbulent combustion is an important part of current heat engine designs; exam-ples of these devices include the common automobile engine and gas turbines. Afundamental challenge today is to reduce the pollutant formation of such enginesas emission regulations become increasingly stringent. However, the complex in-teractions between the turbulent flow, chemical kinetics and heat transfer at engine-relevant pressures and temperatures establish a highly non-linear problem spanningorders of magnitude in both length and time scales. The difficulty here is to intro-duce improvements to these highly calibrated heat engines under such dynamicconditions.Computational fluid dynamics (CFD) is indispensable in the development ofcomplex engines due to its low cost and time requirement compared to experi-ments. While direct numerical simulation (DNS) is the most accurate approach975.1 Introductionto understanding the interactions between turbulence and flame, performing DNSfor engine relevant problems will remain impractical for the foreseeable future.With increasing computational capabilities, numerical simulation of these devicesvia the large eddy simulation (LES) paradigm is becoming more accessible to re-searchers. Nevertheless, due to the strong coupling between turbulence and chem-istry in premixed flames, the prediction of chemical reaction source terms continueto be a modelling challenge.A number of simulation strategies has been developed to model the interactionbetween turbulence and flame [84, 86, 132, 134]. For Conditional Moment Closure(CMC) [84], the transport equations for the reactive species are conditionally aver-aged on a conserved scalar. The primary advantages of this strategy are the reducedspatial dependence of the conditional averages and negligible fluctuations aroundthe conditional averages, which considerably simplify the reaction rate closure. Inaddition, detailed chemistry can be included at relatively low computational costwithout the thin flame restriction encountered in flamelet models. However, one ex-tra transport equation for the conditional mass fraction would need to be solved foreach reactive species, which increases the dimensionality of the problem. Further-more, there are unclosed terms associated with these CMC transport equations thatrequire modelling. Bushe and Steiner [26] later proposed a technique based on theCMC known as the Conditional Source-term Estimation (CSE). Rather than solv-ing extra transport equations, the conditional mass fractions are recovered throughan inverse calculation. A more detailed description of CSE is provided in Sec-tion 5.2. This model has been shown to produce promising results for diffusionflames [67, 93, 161, 183]. Further, Salehi et al. has demonstrated that CSE can alsobe applied to premixed flames and produce meaningful results [42, 143–145, 153].Recent effort has extended CSE to lifted turbulent flames [43] and flames operatingin the Moderate and Intense Low Oxygen Dilution (MILD) combustion mode [92].In this chapter, two optimisations to enhance the capability of CSE have beenintroduced. The first optimisation is to increase the flexibility of the CSE modelby dynamically allocating the physical space required for the individual ensem-ble averages. The second optimisation is to integrate a well established iterativesolver, Least-Squares with QR-factorisation (LSQR) [123], to tackle the inverseproblem. These optimisations can be applied to both Reynolds-averaged Navier-985.2 Conditional Source-term EstimationStokes (RANS) and LES paradigms in the premixed, partially-premixed and non-premixed combustion modes. A comparison between the existing CSE methodolo-gies and the new optimisations is made using LES calculations with OpenFOAM2.1.1 [185] simulating the Gu¨lder burner [190].5.2 Conditional Source-term EstimationCombustion models are required for LES or RANS simulations, which, by their na-ture, cannot fully resolve the turbulence-chemistry interactions spatially and tem-porally. First moment closure, whereby the filtered or averaged reaction rates aredirectly computed by the filtered or averaged values of the scalar fields via the Ar-rhenius rate equations, can lead to considerable discrepencies compared to the fullyresolved values, due to the highly non-linear interactions between chemical sourceterms and the fluctuating scalar fields. Conditional Source-term Estimation is onesuch combustion model that provides closure for the mean chemical source termsfor turbulent flames. CSE invokes the first-order Conditional Moment Closure hy-pothesis such that the conditionally-averaged chemical source terms are closed byevaluating the reaction rate expressions with the conditionally-averaged scalars,ω˙α |c∗ ≈ ω˙α(T |c∗,Yα |c∗,ρ|c∗), (5.1)where c∗ is the discretised conditioning variable. T |c∗, Yα |c∗ and ρ|c∗ are theconditional temperature, mass fraction of species α and density, respectively. Theunconditional mean chemical source term is calculated by integrating the condi-tional source term with the probability density function (PDF) of the conditioningvariable, P(c∗;~x, t), ˜˙ωα = ∫ 10ω˙α |c∗P˜(c∗;~x, t)dc∗. (5.2)The Favre-averaged presumed PDF models adopt the functional form:P˜(c∗;~x, t)≈ P˜(c∗; c˜, c˜′2), (5.3)where c˜ and c˜′2 are the Favre-averaged mean and variance of the progress variable.The usual spatial and temporal coordinates are represented by ~x and t. Several995.2 Conditional Source-term EstimationPDF models are tabulated in this format; a PDF model commonly implemented inconjunction with CSE is the modified laminar flamelet PDF [70, 143–145, 153].The central concept behind CSE is that there is more information available ifthe filtered field is considered in its entirety rather than as unrelated data points.Each filtered volume, governed by the same underlying physical processes, thencontributes additional information about the localised reactive flow. Specifically,the relevant conditional scalar averages are obtained by the inversion of the follow-ing integral equation [26],Y˜α(~x, t) =∫ 10P˜(c∗;~x, t)Yα |c∗(~x, t)dc∗, (5.4)where Y˜α(~x, t) is the Favre-averaged mass fraction1. Figure 5.12 illustrates the in-teraction between the CFD module, the CSE module, the chemistry module —which, in this implementation presented here uses a low-dimensional manifoldchemical mechanism reduction scheme called the Flamelet-Generated Manifold(FGM) [121, 142], which reduces the chemical kinetic mechanism down to a sim-ple two-dimensional table lookup based on two reaction progress variables (themass fractions of vapour water and carbon dioxide) — and the PDF. The massfractions, means and variances of the progress variable are calculated by the CFDmodule via transport equations. These means and variances of the progress vari-able are then used to obtain the PDFs, which in conjunction with the mass fractions,serve as the input to the CSE module via Equation 5.4. The computed conditionalmass fractions are subsequently used to determine the mean reaction rates throughthe FGM and the PDF model. In this chapter, the discussion is focused on the op-timisation of two steps within the CSE framework: (1) to dynamically allocate thephysical boundary that governs the cells used as the input to Equation 5.4; and (2)to compute the conditional mass fractions with greater efficiency.In practice, Equation 5.4 can be expressed in a matrix form for each species,b = Ax, (5.5)1Both Yα |c∗(~x, t) and Yα |c∗ represent the conditional mass fraction of species α . The specifica-tion, (~x, t), expresses the notion that the preceding quantity is a function of space and time. Thissimplified notation is applicable to all of the scalar fields and conditional scalar fields in this section.2 This figure is the same as Figure 2.9. It is repeated here for convenience.1005.2 Conditional Source-term EstimationOpenFoam CFDỸ CO2Ỹ H2Oω˙k |c*Y H2O |c*Y CO2 |c*CSẼ˙ωkc̃ c̃ '2PDFFGMFigure 5.1: CSE operational flowchart.where b, A and x respectively represent the mass fraction, the PDF of the condi-tioning variable and the conditional mass fraction. The length of array b corre-sponds to the number of data points within the ensemble and the length of array xcorresponds to the number of bins in which the PDF is discretised. From past CSEexperience [143–145, 153], it appears that a suitable minimum size of b is aroundO(10,000) while the size of x is typically held constant at 50. These two numbersconsequently determine the size of matrix A, which would be in the vicinity ofO(10,000) by 50. As a result, Equation 5.5 represents an over-constrained inver-sion problem; the solution will therefore be sensitive to small perturbations in band regularisation must be imposed for stability and uniqueness of x [63].The Tikhonov regularisation approach has been implemented with CSE exten-sively in the past for both premixed and non-premixed combustion [42, 67, 92, 93,143–145, 153, 161, 183]. For premixed combustion, it was first proposed by Jinet al. [70] to use the solution of an unstrained one-dimensional laminar premixedflame as the regularisation array. The motivation is that this regularisation methodwould implicitly allow the CSE combustion model to switch to the flamelet solu-tion in case of ill-posedness. In this way, one can obtain a robust combustion modelthat is suitable for both flamelet and non-flamelet combustion regimes.As an aside, another feasible a priori solution to implement as the constraintwith the Tikhonov regularisation would be the conditional mass fraction from theprevious timestep. The advantage of this approach is that perhaps the invertedconditional mass fractions would be more reflective of the current values withinthe computation domain; however, the disadvantage is that the errors from thepreviously inverted solutions may slowly accumulate. As a result of this, it ispossible for this regularisation method to decrease stability instead of providing1015.3 Optimisationsstability as the simulation progresses.In general though, there is a tradeoff between stability and accuracy in theinversion problem [142]. One can always obtain a stable solution to the conditionalmass fraction by increasing the regularisation parameter in the system; however,this is not desired for turbulent combustion as the solution will tend towards theregularisation array: an unstrained one-dimensional laminar premixed flame. Thepreferred approach for achieving stability while minimising regularisation is toincrease the size of the system. This translates to ensuring a minimum number ofdata points is included in each instance of inversion of Equation 5.4. As mentionedpreviously, this minimum number of points is roughly O(10,000).5.3 Optimisations5.3.1 Dynamic Ensemble SelectionAs discussed in the preceeding section, the central concept behind CSE is that thereis more information available if the filtered field is considered in its entirety ratherthan as unrelated data points. This leads to the grouping of localised cells referredto as ensembles. It is important to emphasise the distinction between the CFDsubdomains used by OpenFOAM 2.1.1 [185] and the ensembles used by CSE be-fore describing the new algorithm in detail. The number of subdomains is directlylinked to the number of processors available to the user for the simulation. Thesesubdomains are typically divided in such a way as to minimise boundary surfaceareas (hence, communication overhead) between the processors. Moreover, thenumber of cells within each subdomain is made roughly equal for maximum loadbalance. These subdomains are predetermined by the CFD module and remain con-stant throughout the simulation. With the CSE ensembles, the physical boundariesof these partitions must change during the simulation if one wishes to incorporate aroughly equal number of reactive cells in each ensemble. This would ensure properload balance and more importantly, sufficient number of reactive cells in all of theensembles for the inversion of Equation 5.4.The accepted method of dividing up the domain for CSE is to manually splitup the flow field and assign each ensemble to a processor prior to the simulation;1025.3 Optimisationsthis is a close analogue to the partitioning of the computational domain for parallelprocessing in the CFD module. This manual ensemble division ultimately dependson the physical insight of the user for the particular flow field [143]. For exam-ple, a reasonable set of ensembles for a jet flame would be axisymmetric slicesalong the axial direction. The two major advantages with this method are the easeof implementation and good performance for flames with clear symmetries [143–145, 153]. Manual partition of the domain into geometrically equal ensembles,however, has several associated disadvantages. First, this technique may not besuitable for complex geometries or flames with more complicated structures. Insuch cases, the optimal ensemble assignment might not be immediately evident tothe user prior to the initial simulation. Second, the reactive regions (that is, regionsin which chemical reactions are taking place) for turbulent flames will constantlyshift within the simulation field and will be distributed non-uniformly through-out the domain. This could lead to load balance issues where equal amounts ofcomputational resources are assigned to regions of different reactivity. Figure 5.2illustrates an axisymmetric jet flame in a rectangular simulation domain dividedaxially into different number of ensembles; for the illustrated jet flame, the ensem-bles near the base of the flame (z = 0.0 m) will be less “reactive” than the top ofthe flame. In this case, a cell is defined to be reactive if it satisfies the followingcondition,c˜ · (1− c˜)> 10−2. (5.6)This simple condition eliminates the computational cells that are mostly reac-tants and mostly products. Depending on the inlet conditions and defintion of theprogress variable, the value of c˜ in the inflow boundary condition is typically setto 0 or 1; thus the cells within close proximity to the inflow boundary at the baseof the flame would be considered non-reactive according to Equation 5.6. Third,due to the nature of the CSE inverse problem, one must provide a certain numberof data points in order to invert Equation 5.4 appropriately. This issue becomesof concern when a large number of processors is to be used for the simulation.From Figure 5.2, it is evident that increasing the number of processors effectivelydecreases the thickness of each ensemble; when the thickness becomes too small,the number of available reactive cells to perform the inversion would become in-1035.3 Optimisations(a) 10 Ensembles (b) 25 Ensembles (c) 50 EnsemblesFigure 5.2: Computational domain incorporating different number of staticensembles. The volume of each ensemble decreases with increasing numberof processors.sufficient. A practical method to alleviate part of this problem is to introduce anoverlap in the physical spaces between neighbouring ensembles [143–145, 153].This increases the size of the ensembles and includes an added benefit of a moregradual transition of the conditional mass fraction between ensembles as some datais now shared — albeit at the expense of increased communication overhead be-tween the processors sharing data. Clearly, however, this does not provide a longterm solution to the problem.These disadvantages suggest that a more optimal method to decompose theflow field would be to allocate the ensembles and processors in real time accord-ing to the location of the flame for maximum load balance. As a consequence,the physical boundaries of the ensembles will vary and thus likely lead to a morecomplicated implementation, in addition to increased communication of cell in-formation between processors. The main difficulty here is determining a robustalgorithm to construct the ensembles.When constructing dynamic ensembles, the number of processors, N, available1045.3 OptimisationsN⋮21M12⋮(a) N < M⋮21M12⋮N12⋮(b) N = MM12⋮12⋮N⋮(c) N > MFigure 5.3: The three possible scenarios that one would encounter duringthe construction of dynamic ensembles are illustrated: N < M, N = M andN > M, where N and M are respectively the number of processors availableand the number of ensembles required.may be less than, equal to, or greater than the number of ensembles, M. For anygiven simulation, N will be fixed, where as the number of ensembles (M) and thephysical boundaries of each ensemble will continuously change within the domain.The goal is to divide the domain such that O(10,000) reactive cells are incorpo-rated into each ensemble. It is then conceivable for the algorithm to experienceall three scenarios over the course of one simulation. Thus, it is imperative to de-vise strategies to manage each case as it arises. Figure 5.3 illustrates the N and Mdivisions in a rectangular simulation domain for the three possible cases.For the first case, where the number of processors is less than the number ofensembles (N < M), the implementation could be quite straightforward: one couldsimply divide the number of ensembles as evenly as possible between the proces-sors. If N is much smaller than M, the remainders in division can be almost uni-formly distributed and load balancing would not be a primary issue. If N is slightlysmaller than M, then the difference in computational time between the quickest andslowest processors would be approximately one ensemble inversion.The second case, where the number of processors equals the number of ensem-bles (N = M), is ideal: each processor can be assigned to one ensemble. This iswhat often occurs in conventional CSE simulations, where the user would carefullyselect the number of processors appropriate to the specific problem.The third case, where the number of processors is greater than the number ofensembles (N >M), represents a new frontier in CSE. Here, one could conceive oftwo different approaches, depending on the latency associated with message pass-1055.3 Optimisationsing on the particular machine on which the calculation is to be done. The firstwould be to simply have one local master processor collect the information neededfor the CSE calculation for its neighbours, perform the CSE calculation for the en-semble, then distribute the calculated conditionally averaged reaction rates back tothe neighbouring processors. In this way, the conventional approach to CSE wouldstill be applicable. However, computational resources are not efficiently allocatedas a large number of processors will be at an idle state while only the local mas-ters perform the inversions for the ensembles: the resulting code would not scale.The other approach would be to use a different, parallel algorithm to take advan-tage of the parallel computing environment by having the CSE calculation occurin parallel for ensembles spread over a cluster of processors. During practical im-plementation, this would amount to solving each ensemble in parallel to maximiseefficiency and hopefully lead to a code that would scale. As largely-parallel com-putation becomes a standard in LES combustion codes, the N >M case is likely toemerge as the only practical scenario for turbulent flame simulations with CSE inthe future.In this chapter, two dynamic ensemble selection strategies that automaticallyadjust the physical boundaries of the ensembles as required to maintain the correctnumber of reactive cells in each ensemble are devised.Semi-Dynamic Ensemble Selection AlgorithmThe premise of this new algorithm is relatively straightforward. To begin, onemust define what it is for a cell to be “reactive” — and, hence, to be meaningfullycontributing to the CSE inversion process by virtue of having a PDF for the con-ditioning variable that is not merely a delta function at either zero or unity. Forpremixed combustion, a reactive cell is one where the quantity c˜ · (1− c˜) is greaterthan 10−2, as previously defined by Equation 5.6.Prior to the first timestep of the simulation, the domain is decomposed by theOpenFOAM CFD module; the centroid of each subdomain is computed and passedto the master processor. The master then constructs an N × K array, where N isthe number of processors and K is the maximum number of processors with whichto form an ensemble, that indicates the nearest neighbours to each of the proces-1065.3 Optimisations1st 2nd 3rd →Kth1 2 3 4 →2 1 3 4 →3 2 4 1 →↓ ↓ ↓ ↓ →N N-1 N-2 N-3 →Processor Reactive Cells1 3,5412 4,3343 5,8724 3,546↓ ↓N 1,235Nearest NeighboursProcessorFigure 5.4: Semi-Dynamic Ensemble Selection processor map.sors (Figure 5.4). During runtime, an integer containing the number of reactivecells within each subdomain is sent to the master processor and sorted according tothe N × K array. For example, the construction of the ensemble for processor 3 isdemonstrated by Figure 5.4. The algorithm sums the number of reactive cells in or-der of proximity, starting with itself, then processors 2, 4, 1 and so on, until the totalnumber of reactive cells exceeds the minimum predefined number; the data fromthese processors constitutes the ensemble of processor 3. This summation pro-cedure is repeated for each subdomain until all of the ensembles are constructed.Indeed, the number of reactive cells in each subdomain changes at every timestep;the integers shown in the “Reactive Cells” column in Figure 5.4 are for illustrativepurposes only. Unless otherwise specified, all of the inter-processor communica-tions are performed using the standard OpenFOAM class, UPstream [60], which isbased on the Message Passing Interface (MPI) system.Effectively, each CFD subdomain becomes the centre of a CSE ensemble andits nearest neighbours will contribute information to this central core. For par-ticularly reactive regions, perhaps one or two neighbours will contain sufficientreactive cells for a proper inversion; for less reactive regions, cell information isgathered from the surrounding subdomains in increasing radius from the core untila sufficient number of reactive cells is achieved or the maximum number of neigh-bours, K, is exceeded. From past experience, it appears that a sufficient number ofreactive cells for the CSE inversion is on the order of O(10,000). For the currentstudy, this parameter is varied between 5,000 and 15,000 cells depending on thegrid density. This implementation represents a special case of N = M, where thenumber of ensembles is constant in quantity but varies in physical size. This algo-1075.3 Optimisationsrithm is referred to as the Semi-Dynamic Ensemble Selection because the numberof processors always corresponds to the number of ensembles.A possible limitation of this algorithm is that a number of a processors mayrepeat the same CSE calculations in the extreme scenario where only one or twosubdomains contain reactive cells. For this case, it would be possible for O(K)processors to repeat the same inversion. The repeated inversions, however, wouldnot slow down the overall CSE calculation because all of the processors must besynchronised after every timestep; the otherwise idling processors would need towait for the one to two highly reactive subdomains to finish their CSE calculationbefore proceeding. In this regard, the Semi-Dynamic Ensemble Selection algo-rithm simply repeats the calculation on the idling processors without additionalpenalty. This slight inefficiency is one of the motivations to introduce the DynamicEnsemble Selection algorithm, where the load would be equally distributed to allof the processors in this limiting case. Another consideration is that the precedingmethod of ensemble selection is ideally suited for geometries that contain no spa-tial voids and for subdomains that are arranged along the principle flow direction,which is typical of experimental burners. For complex geometries, however, morejudicious decomposition of the domain may be required.Dynamic Ensemble Selection AlgorithmThere are two primary advantages associated with the Dynamic Ensemble Selec-tion algorithm in comparison to the Semi-Dynamic Ensemble Selection algorithmand static ensembles. First, the number of processors used in the simulation can becompletely decoupled from the number of ensembles. This allows for a more flex-ible allocation of resources and reduces complexity for the user during the initialsetup of the simulation. In addition, this allows for a high number of processorsto be used in parallel without introducing exceedingly small ensembles or largeoverlaps between the ensembles. Second, a parallel solver can be implemented tosolve the inverse problem as multiple processors are now assigned to each matrix.The Dynamic Ensemble Selection algorithm consists of six steps with eachstep summarised in Figure 5.5. Similar to the Semi-Dynamic Ensemble Selectionalgorithm, the reactive cells within each CFD subdomain are first determined by the1085.3 Optimisations1 3,5412 4,3343 5,8724 3,5465 6,8956 7,6237 8,1328 7,235N 6,235⋮⋮3,5414,3345,8723,5466,8957,6238,1327,2356,235⋮1357Processor Reactive CellsDetermine EnsemblesLocalMasters⋮2468NLocalSlaves⋮Send reactivecell counts tomasterMaster contacts processorsmasterMPI-LSQRMPI-LSQRMPI-LSQRMPI-LSQRMPI-LSQRDeterminereactivecellsMasterdeterminesensemblesLocal mastersdeterminecommunicationdirectivesLocal masterscontact local slaves1 2 3 4 5 6Figure 5.5: Dynamic Ensemble Selection schematic.criterion set in Equation 5.6. Second, once the number of reactive cells is identifiedon each processor, the N integers are sent to the master processor. Third, the masterallocates the ensembles according to the user defined quantity that is related tothe minimum number of reactive cells per ensemble; as before, this number isof O(10,000). Fourth, the master communicates with all of the processors andassigns their role, local master or local slave, in each ensemble. Fifth, the localmasters compute the communication directives for the local slaves belonging tothe respective ensembles according to the classic round robin algorithm. Sixth, thelocal masters distribute the directives to the appropriate local slaves. This processof ensemble allocation occurs at the beginning of each timestep, prior to any othercalculations.To further consider the layout of the ensembles in this algorithm, the first andlast processors of each ensemble are deliberately overlapped with the previous andnext ensembles. This sharing of reactive cell data between ensembles allows fora smoother transition of the resulting conditional mass fractions from one part ofthe domain to the next. The ensembles illustrated in Figure 5.5 can be examined.For this example, each ensemble consists of data from three processors: Ensemble1 includes data from processors 1–3; ensemble 2 includes data from processors3–5; and ensemble 3 includes data from processors 5–7. This pattern continues1095.3 Optimisationsuntil the N processors are fully divided. To properly invert the matrix in ensemble1, processor 1 must receive the reactive cell data from both processors 2 and 3,while processor 2 must receive the reactive cell data from processors 1 and 3. Onlyprocessors 1 and 2 will be used to compute the inverse of this ensemble in parallel;processor 3 is not directly used in this inversion calculation as it will be used inensemble 2. Similarly, for ensemble 2, processor 3 must receive the reactive celldata from both processors 4 and 5, while processor 4 must receive the reactive celldata from processors 3 and 5. The two processors, 3 and 4, will then compute theinverse of this ensemble in parallel. This pattern continues until the M ensemblesare fully divided. In the event where the last ensemble in the domain contains lessthan the user-defined minimum number of reactive cells, it will be combined withthe nearest ensemble.While this procedure appears to be complex, the additional computational costof the Dynamic Ensemble Selection is relatively low. In the first step, the reac-tive cells must be identified regardless of the ensemble selection method, thus,imposing no additional cost. Second, only one integer is sent to the master fromeach processor to identify the number of reactive cells. Third, the master is onlyrequired to sort N integers, where N is the total number of processors in the simula-tion. Fourth, the master sends three integers back to each processor; these integersidentify the positions of the first, the last and the overlapping processors belong-ing to each ensemble. Fifth, to further reduce computational time, the round robinalgorithm used to determine the communication directives by the local masters ispre-computed and tabulated into lookup tables prior to the simulation. Sixth, thelocal masters send the communication directives to the proper local slaves, witheach ensemble typically incorporating one to three local slaves. These directivesconsist of approximately 2·n integers, where n is the number of processors in eachensemble. Once the ensembles have been assigned, the next step is to compute theinverse solution with the parallel iterative solver described in Section 5.3.2.5.3.2 Matrix Inversion - Iterative SolverInverting the PDF matrix, denoted as A, is the most time-consuming step in a CSEcalculation. As mentioned, the typical inversion matrix isO(10,000) by 50 and the1105.3 Optimisationsaccepted method uses LU-decomposition with a zeroth order Tikhonov regularisa-tion [143–145, 153]. Practical implementation of the LU-decomposition usuallyrequires a mathematical manipulation to first reduce the size of the inversion, in-volving a matrix multiplication of ATA. The benefits of the LU method are that thesolution is exact and the inversion time is deterministic. However, in the contextof CSE, an exact solution is not required for the ensemble average, as there areat least two sources of error implicit to the inversion of Equation 5.4. First, thePDF of the reaction progress variable that is used as the kernel of inversion willinherently contain some modelling error; and second, the zeroth order Tikhonovregularisation using the laminar flamelet solution contributes additional error. Asa consequence, provided the error generated from the iterative solver is relativelysmall compared to the aforementioned sources, it is tenable to obtain a pseudo-inverse of the matrix instead of an exact inverse of the matrix for the conditionalscalars. This relatively modest reduction in accuracy is a small compromise to thesubstantial reduction in computational effort for the matrix inversions.To address this, the LU solver is proposed to be replaced with an establishediterative approach, Least-Squares with QR-factorisation, or LSQR [123]. This it-erative method is designed for the computation of Ax = b, where A and b areknown. The LSQR implements the Lanczos method that leads to a factorization ofthe tridiagonal reduced matrix for ATA. This method is often the method of choicefor over-determined or under-determined systems [140]. Formally, it is equivalentto the following minimisation problems:min || Ax−b ||2, (5.7)min || Ax−b ||2+λ 2d ||x ||2, (5.8)where the matrix A is supposed to be large and sparse3 and λd is the damping co-efficient. The method is based on the bidiagonalization procedure of Golub andKahan [57]. It is analytically equivalent to the standard method of conjugate gra-dients, but possesses more favorable numerical properties [123]. The cost per it-3The PDF matrix is not sparse; however, the LSQR solver has proven to be robust and efficientfor all of the performed tests. These include standalone matrix inversions and CSE inversions in thecontext of LES.1115.3 Optimisationseration for this solver is primarily dependent on the matrix-vector multiplicationsof A·x and AT ·b and typically converges in less than 25 iterations for correctlyregularised one-condition CSE matrices.Another considerable advantage of LSQR over LU-decomposition is that theiterative solver can be parallelised for the computation of large-scale inversionproblems. The implementation of a parallel LSQR solver would allow the userto efficiently solve each CSE inversion problem when the number of processorsexceeds the number of ensembles. A Scalable Parallel LSQR (SPLSQR) algo-rithm has been developed recently by Huang et al. [66] in the field of seismictomography; SPLSQR has demonstrated scalability beyond O(10,000) processorsfor large-scale inversion problems. In this study, a modified version of their scal-able parallel algorithm adapted for CSE problems has been developed; this parallelsolver will be referred to as MPI-LSQR in the subsequent discussions to reduceambiguity. In general, the code structure of the MPI-LSQR and LSQR solvers arevery comparable, though two unique distinctions can be identified. First, each pro-cessor participating in the parallel computation of MPI-LSQR is responsible onlyfor a portion of the solution vector. Second, the partially updated solution vectorswould need to be distributed between all of the participating processors at eachiteration.Figure 5.6 illustrates the steps required by the MPI-LSQR algorithm to obtain aconverged solution for the case of three participating processors. The most currentworking arrays of length x and b, respectively representing Yα |c∗ and Y˜α , are firstgathered and distributed between the processors. These arrays, shown in white,are then checked for convergence according to the user defined error tolerances.If the arrays satisfy the error tolerances, the current solution is accepted. If not,then an iteration of the solver will be performed; each iteration consists of twomatrix-vector multiplications on the order of A·x and AT ·b, followed by two smallnormalisation calculations. Crucially, each processor only needs to compute onethird of the total multiplications for this particular illustration. Each processorthen obtains a local portion of the new x and b arrays. Next, these local arraysare gathered and distributed between the processors. This cycle continues untilconvergence is achieved. In applying the MPI-LSQR solver to the inverse problem,a computational scheme that is parallel in both the CFD and the CSE modules has1125.4 Numerical Methodsx bDistribute Converged?NoA1A2A3b2b3b1=A1 A2 A3 = x2x3x1xbComputeb1x2 b2b3x3YesFinishx1Figure 5.6: MPI-LSQR iteration. Each n participating processor performs1/n of the total matrix-vector multiplications.been obtained.5.4 Numerical Methods5.4.1 Standalone Inversion TestsThe LSQR and MPI-LSQR solvers have been thoroughly tested before implemen-tion with the CSE calculations. Two separate tests have been performed. First,the original formulation of LSQR by Paige and Saunders [123] has been modi-fied to accomodate the input of an initial solution, x0. The tests have shown thatthis improves the convergence rate by 10–15%, reducing the number of iterationsto 20–25 for the average one-condition CSE problem. During practical imple-mentation, x0 is set to the conditional mass fraction of the previous timestep; theusual zeroth order Tikhonov regularisation with the laminar flamelet solution isused [142]. Second, the LSQR test is repeated with the MPI-LSQR solver sub-jected to the same initial condition and Tikhonov regularisation protocols. Duringcombustion simulation with CSE, in the case where the number of processors ex-ceeds that of the ensembles (N > M), it is possible for each of the processors tocarry a different conditional mass fraction from the previous timestep into the en-semble. In this case, the conditional mass fraction used as the input of x0 would1135.4 Numerical Methodshave to be a weighted average of the participating processors.For the standalone LSQR inversion test, five matrices of different dimensionsare examined. The PDF matrix, A, is first randomly generated and multiplied witha known vector, x, to create the sample problem of the form Ax = b. The A andb are then inverted individually by the LU-decomposition (as the control) and theLSQR solver at various settings. In order to better represent the conditions thatthe LSQR solver would encounter during practical CSE implementation, the ini-tial condition input, x0, is generated by adding random +/- 5% fluctuations to theknown x. These random fluctuations emulate the conditional mass fraction changeper timestep. Though this is likely an overestimate of the change, it serves to pro-vide an upper limit of the convergence time. Moreover, the error bounds controllingthe termination condition of the LSQR solver have been adjusted to determine theireffects on accuracy and convergence time. The +/- 5% fluctuations are not neces-sary for the LU-decomposition method as it does not accept an input solution andthe computational times are strictly dependent on the size of the system. The inver-sions are performed on a single i7-2600K processor with a maximum frequency of3.8 GHz.For the standalone MPI-LSQR inversion test, quantifying the scalability of thenew MPI-LSQR algorithm is the primary concern. The inversion problems aregenerated identically as the LSQR inversion test and the parallel solver is sub-jected to the same initial condition and Tikhonov regularisation protocols. Theinversions have been computed using a number of different computational node(1–12) and processor (1–48) combinations. All of the parallel computations in thisstudy are performed on the Mammouth Paralle`le II supercomputer of Universite´ ofSherbrooke. The participating nodes consisted of processors with frequencies of2.2 GHz.5.4.2 LES of Axisymmetric BurnerThe simulation setup of the stoichiometric flame conditions described in a previousstudy [153] has been followed. The CFD code is based on OpenFOAM 2.1.1 [185]and a FGM model is used for the chemistry reduction [121, 142]. Second orderimplicit temporal and second order spatial schemes are employed. The LES model1145.4 Numerical Methodsis the standard oneEqEddy from OpenFOAM, which is equivalent to the k-equationeddy-viscosity model. A cylindrical simulation domain with a radius of 0.025 mand a height of 0.10 m is used along with three grids of different hexahedral cellcounts and processor numbers. The three grids consisting of 52,500, 284,000 and880,000 cells are respectively assigned to 4, 24 and 48 processors. While the griddensity is lower than required for accurate LES simulations, demonstrating thecharacteristics of the new optimisation algorithms takes precedence. The goal hereis to compare the combined performance of the three inversion strategies underidentical simulation conditions: MPI-LSQR solver with Dynamic Ensemble Se-lection, LSQR solver with Semi-Dynamic Ensemble Selection and conventionalLU-decomposition with static (manual) ensemble selection.The inflow condition is isotropic homogeneous turbulence generated from aspectral method using the energy spectrum of Haworth and Poinsot [65]. The tur-bulent kinetic energy (k), peak energy wavenumber and maximum wavenumberare respectively 33.3 m2/s2, 2pi/l0 and 2pi/2∆, where l0 and ∆ are the integral lengthand grid spacing. The integral length is set to 4.8·10−3 m and the grid spacingvaried from 8.96·10−4 m to 2.24·10−4 m. The time correlation for the velocityfluctuation is based on the method developed by Billson [15]. The relevant simu-lation parameters are summarised in Table 5.1. The inlet diameter is 1.12·10−2 mwith stoichiometric premixed methane-air reactants at 300 K and 101.3 kPa. Theseparameters are designed to simulate the experimental conditions as closely as pos-sible [153]. The inlet bulk flow, which includes the spectral generated turbulence,is 15.58 m/s. A laminar, hot pilot surrounds the inlet with an inflow of combustionproducts moving at 16.81 m/s.The simulations are carried out in two stages. The first stage is performed bythe very coarse grid (52,500 cells) to quickly discharge the initial condition andto develop the flame in the domain. This is allowed to carry on for approximatelyeight flow-through times, from t = 0.0 s to t = 0.05 s. The second stage is initiatedby interpolating the result of the very coarse grid to the two finer grids. Data isthen collected from t = 0.05 s to t = 0.25 s for the three grids. Radial temper-ature profiles are collected every 0.001 s at axial positions of z/D = 1–6, whereD is the diameter of the inlet (D = 1.12·10−2 m). The simulations are performedusing the MPI-LSQR solver with Dynamic Ensemble Selection, LSQR solver with1155.5 Results and DiscussionGrid 52,500 cells 284,000 cells 880,000 cellsk 33.3 m2/s2 33.3 m2/s2 33.3 m2/s2Peak energy wavenumber 1.31 · 103 m−1 1.31 · 103 m−1 1.31 · 103 m−1Maximum wavenumber 3.51 · 103 m−1 7.01 · 103 m−1 1.40 · 104 m−1l0 4.80 · 10−3 m 4.80 · 10−3 m 4.80 · 10−3 m∆ 8.96 · 10−4 m 4.48 · 10−4 m 2.24 · 10−4 m∆t 8.00 · 10−6 s 4.00 · 10−6 s 2.50 · 10−6 sτt 9.96 · 10−4 s 9.96 · 10−4 s 9.96 · 10−4 sa 0.992 0.996 0.998b 0.125 0.0884 0.0700Table 5.1: Relevant LES simulation parameters. The values of a and b arecomputed according to the approach of Billson [15]: a = exp(−∆t/τt) andb=√1−a2, where ∆t and τt are respectively the timestep size and the integraltime scale.Semi-Dynamic Ensemble Selection and LU-decomposition with static ensemblesfor the three grids. It is important to emphasise that the only changes between thesimulations are the inversion and ensemble selection algorithms; the grids, CFDmodule, chemistry reduction, timestep sizes and computing facilities are otherwiseidentical. The performance and results are documented in Section 5.5.2.5.5 Results and Discussion5.5.1 Standalone Inversion TestsTable 5.2 compares the computational cost of inverting matrices of various sizesand L2-norms4 of the solutions that are relevant to one and two-condition5 CSEproblems using the LU-decomposition method and the LSQR solver. Figure 5.7graphically compares the computational cost of inversion from Table 5.2. It isapparent that the LSQR solver demonstrates noticeably shorter calculation timesthan LU-decomposition for one-condition CSE and is orders of magnitude faster4L2-norm is defined as: L2 =√n∑i=1(xi,LSQR− xi,exact)2, where n is the number of discrete bins inthe tabulated PDF model, which is typically 50 for one-condition CSE problems.5Two-condition CSE is also referred to as Doubly Conditional Source-term Estimation (DCSE)in the literature [43].1165.5 Results and Discussion10,000 x 5000.010.020.03Inversion Time (s)30,000 x 5000.10.20.32,000 x 1,0000246810x6,000 x 1,0000102030x100,000 x 2,500020406080100x x→ 4 hours + →xFigure 5.7: Comparison of inversion times between LU-decomposition andLSQR at a number of different error tolerance settings (|e|). This is a graphicalrepresentation of the data in Table 5.2. Each of the five plots represents adifferent matrix size. For each plot, the bars represent inversion times of (fromleft to right): LU, LSQR (~0) (|e| = 10−10), LSQR (x0) (|e| = 10−10), LSQR(x0) (|e| = 10−8), LSQR (x0) (|e| = 10−6), and LSQR (x0) (|e| = 10−4). The‘×’ represents the failed cases.for two-condition CSE problems6. These results are anticipated as it has been doc-umented that the LSQR solver is well-suited and exhibits exceptional convergencebehaviour for over-determined inversion problems [41, 140]. Another favourablefeature of LSQR is permitting the user to determine the relative and absolute errorbounds according to the problem. For typical one-condition CSE problems, it isconceivable to increase the error tolerances to produce solutions with L2-norms onthe order of 10−5, further decreasing the convergence time. However, it should bementioned that as the error tolerances are increased, some cases begin to revealconvergence issues as expected. These cases are marked by ‘×’.Table 5.3 compares the computational cost of inverting matrices of varioussizes and L2-norms of the solutions that are relevant to one and two-condition CSEproblems using the MPI-LSQR implementation. Figure 5.8 graphically comparesthe computational cost of inversion from Table 5.3. The inversions in the first6 The focus is on one-condition CSE problems in this study. However, the LSQR solver is inde-pendent of the inversion problem and it is meaningful to illustrate its potential to solve two-conditionCSE matrices for future work.1175.5 Results and DiscussionDimensions 1·104×50 3·104×50 2·103×1·103 6·103×1·103 1·105×2.5·103LU 0.0319 s 0.301 s 9.97 s 30.7 s 4 hours+LSQR (~0)|e| = 10−10 0.0258 s 0.0912 s 1.23 s 1.95 s 98.2 sL2-norm→ 4.5·10−10 4.4·10−10 8.4·10−6 1.3·10−6 1.3·10−6LSQR (x0)|e| = 10−10 0.0253 s 0.0875 s 1.22 s 1.90 s 93.6 sL2-norm→ 1.5·10−11 3.0·10−11 5.0·10−7 4.7·10−8 8.9·10−8LSQR (x0)|e| = 10−8 0.0216 s 0.0753 s 0.975 s 1.54 s 75.3 sL2-norm→ 3.6·10−8 8.9·10−10 4.3·10−5 4.6·10−6 4.2·10−6LSQR (x0)|e| = 10−6 0.0168 s 0.0639 s 0.699 s 1.26 s ×L2-norm→ 8.6·10−6 7.4·10−8 3.9·10−3 3.7·10−4LSQR (x0)|e| = 10−4 0.0142 s 0.0491 s × × ×L2-norm→ 1.4·10−4 1.6·10−5Table 5.2: Comparison of inversion times between LU-decomposition andLSQR at a number of different settings. The |e|,~0 and x0 respectively denotethe error tolerances in the solver, a zero initial solution and a non-zero initialsolution. The associated L2-norms are printed below the convergence times.Typical matrix sizes for one-condition CSE are represented by the first twocolumns while typical matrix sizes for two-condition CSE are represented bythe last three columns. Failed cases are denoted by ‘×’, where lower errortolerances are required to provide stable solutions. From the smallest to thelargest matrix, 1,000, 500, 100, 20 and 2 samples are used to obtain the aver-age inversion times.1185.5 Results and Discussion10,000 x 5000.10.2Inversion Time (s)30,000 x 5000.10.20.30.40.50.60.72,000 x 1,0000246810126,000 x 1,0000102030405060100,000 x 2,5000100200300400xFigure 5.8: Comparison of inversion times with the MPI-LSQR implementa-tion using a number of different processors counts. This is a graphical repre-sentation of the data in Table 5.3. Each of the five plots represents a differentmatrix size. For each plot, the bars represent inversion times of (from left toright): LU, LSQR (~0), LSQR (x0), (1 × 2), (2 × 1), (1 × 3), (3 × 1), (1 × 4),(4 × 1), (1 × 8), (8 × 1), (1 × 12), (12 × 1), (1 × 24), and (1 × 48). Thenumbers (x × y) represent the number of nodes (x) and the number of proces-sors per node (y) used in the MPI-LSQR computation. The inversion for thematrix of size 100,000 × 2,500 is not performed by the LU-decompositiondue to an exceedingly exhaustive computational requirement; it is marked by‘×’.three rows are performed with one processor via the LU-decomposition and LSQRsolver; these are considered the reference states. The inversions from rows 4–15are performed by the MPI-LSQR solver with different numbers of nodes and pro-cessors. Moreover, the error tolerence (|e|) is held constant at 1·10−10 for all of theLSQR and MPI-LSQR cases for consistency. The L2-norms produced by this valueof error tolerence would be much lower than required by CSE; as such, the listedtimes illustrate an upper bound for the inversion calculations. During practicalimplementation, one would increase the |e|-tolerence to the vicinity of 1·10−8–1·10−6 to increase the convergence rate by approximately 30%. The method usedto generate the sample problems and initial condition inputs for rows 3–15 are aspreviously described in Section 5.4.1.In general, it can be observed from Table 5.3 and Figure 5.8 that the inversion1195.5 Results and DiscussionDimensions 1·104×50 3·104×50 2·103×1·103 6·103×1·103 1·105×2.5·103LU 0.228 s 0.694 s 12.3 s 59.0 s N/ALSQR (~0) 0.0877 s 0.249 s 2.03 s 5.39 s 313 sL2-norm→ 4.1·10−10 4.3·10−10 8.6·10−6 1.3·10−6 1.3·10−6LSQR (x0) 0.0881 s 0.337 s 2.80 s 5.41 s 361 sL2-norm→ 1.6·10−11 2.9·10−11 4.9·10−7 5.5·10−8 1.6·10−81 × 2 0.0717 s 0.205 s 1.85 s 2.98 s 126 sL2-norm→ 1.3·10−11 3.1·10−11 5.0·10−7 5.7·10−8 1.6·10−82 × 1 0.0739 s 0.197 s 1.88 s 2.72 s 137 sL2-norm→ 1.3·10−11 2.9·10−11 5.0·10−7 5.6·10−8 1.5·10−81 × 3 0.0532 s 0.156 s 1.01 s 1.60 s 86.1 sL2-norm→ 1.3·10−11 2.9·10−11 4.9·10−7 6.1·10−8 1.5·10−83 × 1 0.0573 s 0.168 s 0.978 s 1.50 s 86.5 sL2-norm→ 1.3·10−11 2.9·10−11 4.9·10−7 5.0·10−8 1.6·10−81 × 4 0.0503 s 0.137 s 0.574 s 1.08 s 67.0 sL2-norm→ 1.2·10−11 3.0·10−11 4.9·10−7 5.8·10−8 1.8·10−84 × 1 0.0460 s 0.134 s 0.574 s 0.991 s 63.7 sL2-norm→ 1.1·10−11 3.0·10−11 4.9·10−7 5.7·10−8 1.6·10−81 × 8 0.0386 s 0.107 s 0.242 s 0.573 s 47.2 sL2-norm→ 1.4·10−11 3.0·10−11 4.8·10−7 5.6·10−8 1.7·10−88 × 1 0.0398 s 0.0998 s 0.229 s 0.571 s 50.5 sL2-norm→ 1.2·10−11 2.9·10−11 4.9·10−7 5.1·10−8 1.6·10−81 × 12 0.0347 s 0.106 s 0.192 s 0.512 s 44.3 sL2-norm→ 1.1·10−11 2.9·10−11 4.8·10−7 6.0·10−8 1.7·10−812 × 1 0.0337 s 0.0978 s 0.202 s 0.559 s 38.3 sL2-norm→ 1.7·10−11 2.9·10−11 4.9·10−7 5.4·10−8 1.5·10−81 × 24 0.0311 s 0.0927 s 0.161 s 0.392 s 24.3 sL2-norm→ 1.3·10−11 2.9·10−11 5.0·10−7 6.1·10−8 1.7·10−81 × 48 0.0295 s 0.0983 s 0.152 s 0.390 s 23.6 sL2-norm→ 1.9·10−11 3.0·10−11 4.9·10−7 5.8·10−8 1.7·10−8Table 5.3: Comparison of inversion times for single ensembles of varyingsizes with the MPI-LSQR implementation using a number of different pro-cessors. The numbers x × y in the first column represent the number of nodes(x) and the number of processors per node (y) used in the computation. Theerror tolerances are kept constant at 10−10 between all of the cases for consis-tency.1205.5 Results and Discussiontimes decrease steadily, to a certain limit, with increasing number of processors.The limit appears to be around 24–48 processors, where gains tend to be mini-mal with additional processors. For some cases, the inversion time of LSQR (x0)is slightly longer than LSQR (~0). However, this is simply a consequence of themethod by which the terminating condition is calculated from the relative errortolerance in the LSQR solver. An inspection of Table 5.3 reveals that the L2-normfrom LSQR (x0) is typically an order of magnitude smaller than the L2-norm ofLSQR (~0), leading to longer inversion times. With the inclusion of an initial con-dition (x0), it is possible to decrease the relative error tolerance in the solver toachieve the same L2-norm as cases with a zero initial condition (~0). Nonetheless,for consistency between cases, the relative and absolute error tolerances are held atthe same values.Larger systems of equations (two-condition CSE matrices) tend to demonstratemore benefit from increasing the processor count, though this is anticipated ascommunication overhead differentially affects smaller systems. A rather surpris-ing result is that the internodal and intranodal communication times appear to becomparable for the machine used. For instance, examining the 1× 12 and the 12×1 inversion times, where the first number represents the number of nodes while thesecond number represents the number of processors used per node, the values aretypically within 15% of one another. More importantly, inversions performed byprocessors centralised on one node do not consistently yield greater performancethan processors scattered over a number of nodes. This implies that the ensem-ble assignment could be very flexible for future implementation of this solver withCSE, as the penalty for computing a parallel inversion with processors from differ-ent nodes appears to be small.5.5.2 LES of Axisymmetric BurnerThe temperature and species mass fraction profiles at six axial locations for thenine simulations are depicted by Figure 5.9. For each grid density, the results aver-aged over 200 ms for the LU-decomposition with static ensembles and LSQR withSemi-Dynamic Ensemble Selection are nearly indistinguishable from one another;this implies that the conditional mass fractions produced by both of the methods1215.5 Results and Discussionare very comparable. The slight deviations could be caused by a difference inthe ensemble sizes between the static ensemble and Semi-Dynamic Ensemble Se-lection methods; additionally, the spectral turbulent inlet condition could play aminor role. It is interesting, however, to observe the transition exhibited by theCO mass fraction profiles with increasing grid resolution; particularly, the CO ap-pears to increase around r/D = 0.6 farther downstream from the burner nozzle.This behaviour of the CO mass fraction could be associated with the currentlyimplemented chemistry reduction technique, where the FGM may not fully cap-ture the formation of species with slow time scales. If that is the case, it wouldnot be possible for the algorithm, in its present form, to accurately represent theCO profiles. As the main purpose of this study is to compare the three ensem-ble selection algorithms and the three inversion strategies, the LES results fromthe LU-decomposition with static ensembles and the LSQR with Semi-DynamicEnsemble Selection are remarkably similar.The temperature and species mass fraction profiles produced by the MPI-LSQRand Dynamic Ensemble Selection are slightly different in comparison with theother two algorithms. This can be attributed to at least two possible factors. First,the total number of ensembles no longer corresponds to the number of processorsin the MPI-LSQR cases, as explained in Section 5.3.1. The number of ensemblesfor the 284,000 grid simulation is typically around 10–15, instead of 24 for the LUand LSQR simulations. Similarly, the number of ensembles for the 880,000 gridsimulation is typically around 20–25, instead of 48 for the LU and LSQR simula-tions. This change in the number of ensembles in the domain would have an impacton the conditional mass fractions, which would then alter the overall reaction ratesin each CFD subdomain, leading to different temperature profiles. Second, theconditional mass fraction of the previous timestep that is used as the initial condi-tion input is now a weighted average of the participating CFD subdomains in eachensemble instead of one CFD subdomain. While this change in the initial conditionwill not lead to a noticeable difference in the solution over one timestep, the differ-ence would accumulate over thousands of timesteps. Despite these two relativelyimportant modifications, the MPI-LSQR results are very comparable to the LU andLSQR results. In general, it appears that the number of ensembles and the exactlayout of the ensembles do not significantly alter the converged solution. Similar1225.5 Results and DiscussionFigure 5.9: The radial temperature and species mass fraction profiles at sixaxial locations downstream of the inlet for the three grids are shown. The pro-files at z/D = 1 are shown in black, then each successive jet diameter down-stream is shown in a lighter shade of gray, up to z/D = 6. The jet diameter,D, is 1.12·10−2 m. The symbols represent: LU–manual: ‘—’; LSQR–semi-dynamic: ‘4’; PLSQR–dynamic: ‘•’.1235.5 Results and Discussionto the LU and LSQR results, it is observed that the CO mass fraction exhibits alocal maximum near r/D = 0.6 with increasing grid resolution farther downstreamfrom the burner nozzle.The computational time for the nine simulations are tabulated in Table 5.4.The execution times indicate that the LSQR solver in conjunction with the Semi-Dynamic Ensemble Selection method can reduce the total computational cost byup to 50% in comparison with the conventional LU-decomposition method. More-over, the MPI-LSQR solver with Dynamic Ensemble Selection can further reducethe computational times by another 2–22% in comparison with LSQR7. While theperformance increase is relatively small by changing the implementation from theLSQR solver with the Semi-Dynamic Ensemble Selection to the MPI-LSQR solverwith Dynamic Ensemble Selection, the primary advantage is to enable the CSEcombustion module to continuously form more suitable ensembles as the flameevolves in shape and position within the domain. From these simulations and thestandalone tests, it is conceivable that a parallelised version of LSQR along withthe Dynamic Ensemble Selection algorithm could serve as a good candidate to re-place the LU-decomposition with static ensembles for future versions of CSE forall modes of combustion.7 The execution times for the 880,000 cell cases are averaged over two runs to eliminate possiblenetwork instability errors.1245.6 ConclusionGrid 52,500 52,500 284,000 880,000Sim. time (s) 0.0→ 0.05 0.05→ 0.25 0.05→ 0.25 0.05→ 0.25∆t (s) 8·10−6 8·10−6 4·10−6 2.5·10−6LU exec. (hrs) 1.7 11.9 88.8 302.0Processors 4×3.8 GHz 4×2.1 GHz 24×2.1 GHz 48×2.2 GHzLSQR exec. (hrs) 1.2 8.6 37.5 177.2Processors 4×3.8 GHz 4×2.1 GHz 24×2.1 GHz 48×2.2 GHzMPI-LSQR exec. (hrs) 0.98 7.4 35.2 173.6Processors 4×3.8 GHz 4×2.1 GHz 24×2.1 GHz 48×2.2 GHzTable 5.4: Summary of simulation and execution times. The execution timesfor LSQR with Semi-Dynamic Ensemble Selection are typically 50% of theLU-decomposition method for larger systems. The executions times for MPI-LSQR with Dynamic Ensemble Selection are further decreased by 2–22% incomparison with the LSQR procedure. The ensemble sizes were decreased forthe coarse grid (52,500) case in order to maintain locality for cell information,leading to more comparable execution times.5.6 ConclusionA first study of Conditional Source-term Estimation using a modified version ofthe LSQR iterative solver along with Dynamic Ensemble Selection for a premixed,turbulent methane-air Bunsen burner has been performed. The LSQR solver hasbeen modified to accept an initial condition as the input and parallelised for scala-bility. The standalone tests and LES results indicate that the iterative solver is fullycapable of reproducing similar solutions to the LU-decomposition method whilesubstantially reducing computational time. The data sharing algorithm betweenensembles has also proven to be robust in achieving the proper communicationdirectives for grids of three sizes, three different processor counts and differentcomputational hardware. The integration of these two algorithms allows one toestablish the framework for a highly autonomous and parallel version of CSE thatcan adapt to the changing flame distribution in complex geometries.125Chapter 6LES of Turbulent PremixedFlames with Optimised CSE andLEM SubmodelsThe work in this chapter serves to unify the submodels and optimisations developedin Chapters 4 and 5. A turbulent, premixed, stoichiometric methane-air flame of alaboratory-scale Bunsen burner operating in the thin reaction zone has been sim-ulated via large eddy simulation (LES). The optimised and parallelised version ofConditional Source-term Estimation (CSE) developed in Chapter 5 has been imple-mented; furthermore, the new pseudo-turbulent probability density function (PDF)and scalar dissipation rate (SDR) models formulated from the Linear-Eddy Model(LEM), as discussed in Chapter 4, are used to complete the CSE inversion and toclose the transport equation for the variance of the progress variable, respectively.6.1 IntroductionAs described in the previous chapter, a number of simulation strategies has been de-veloped to model the interaction between turbulence and flame [84, 86, 132, 134].For Conditional Moment Closure (CMC) [84], the transport equations for the reac-tive species are conditionally averaged on a conserved scalar. The primary advan-tages of this strategy are the reduced spatial dependence of the conditional averages1266.1 Introductionand negligible fluctuations around the conditional averages, which considerablysimplifies the reaction rate closure. In addition, detailed chemistry can be includedat relatively low computational cost without the thin flame restriction encounteredin flamelet models. However, one extra transport equation for the conditional massfraction would need to be solved for each reactive species, increasing the dimen-sionality of the problem. Furthermore, there are unclosed terms associated withthe CMC transport equations that require modelling. Bushe and Steiner [26] laterproposed a technique based on the CMC known as the Conditional Source-termEstimation (CSE). Rather than solving extra transport equations, the conditionalmass fractions are recovered through an inverse calculation. A detailed descriptionof the CSE combustion model has been provided in Section 5.2. This model hasbeen shown to produce promising results for diffusion flames [67, 93, 161, 183].Further, Salehi et al. has demonstrated that CSE can also be applied to premixedflames and produce meaningful results [42, 143–145, 153]. Recent effort has ex-tended CSE to lifted turbulent flames [43] and flames operating in the Moderateand Intense Low Oxygen Dilution (MILD) combustion mode [92].In this chapter, the performance of an optimised version of the CSE combustionmodel designed for parallel computations over a cluster of processors via a largeeddy simulation of a turbulent, premixed methane-air laboratory burner operatingin the thin reaction zone [173] is evaluated. In addition, recently developed proba-bility density function (PDF) and scalar dissipation rate (SDR) submodels requiredto complete the CSE model and the variance transport equation for the progress ofreaction are implemented concurrently in a large eddy simulation combustion codefor the first time. The pseudo-turbulent PDF and SDR submodels are formulatedthrough the Linear-Eddy Model (LEM) [148, 157] and benefit from the inclusion ofturbulence and correct chemistry in the framework compared to previous submod-els; as a result, the LEM-based PDF model has demonstrated improved accuracyover previous modelling approaches [175]. In the following section, a detailed de-scription of the underlying theory behind each numerical model used in this workwill be provided. Next, the simulation setup for both the non-reactive and reactiveflows will be described. This will be followed by the LES result and concludingremarks.1276.2 Theory6.2 Theory6.2.1 Governing EquationsThe governing equations for the gaseous flow field is based on the standard Open-FOAM 2.1.1 package [185], which would be similar to the ones outlined in Sec-tion 2.2. The energy equation, however, is solved in the form of the total enthalpy.To integrate the CFD code with the CSE combustion module, two additional trans-port equations for the mean and variance of the reaction progress variable are re-quired. The transport equation of the normalised mean of the progress variabletakes the following form [132],∂ (ρ¯ c˜)∂ t+∂ (ρ¯ u˜ic˜)∂xi=∂∂xi(ρ¯νTSct∂ c˜∂xi)+ ω˙c, (6.1)where Sct is the turbulent Schmidt number1. The mass fraction of CO2 is used asthe progress variable in this study and the normalisation is calculated in accordancewith Equation 2.73. The transport equation of the variance of the progress variablewith some simplification reads [132],∂ (ρ¯ c˜′2)∂ t+∂ (ρ¯ u˜ic˜′2)∂xi=∂∂xi(ρ¯νTSct∂ c˜′2∂xi)+2ρ¯νTSct∂ c˜∂xi∂ c˜∂xi−2ρDc ∂c∂xi∂c∂xi+2c′ω˙c,(6.2)where ρDc ∂c∂xi∂c∂xi ≈ ρ¯ χ˜c 2; χ˜c is the scalar dissipation rate of the progress vari-able. An implicit, second order temporal scheme from the standard OpenFOAMpackage [185], designated as backward, is implemented in the calculation. TheLES subgrid model is known as the oneEqEddy in the framework of OpenFOAM,which is equivalent to the k-equation eddy-viscosity model.6.2.2 Chemistry ReductionTo maintain reasonable computational times for combustion simulations, severalnotable techniques have been developed in the past decades to simplify the reac-1Sct is set to 0.7 unless otherwise specified.2In some literature, the factor of 2 could be included in the definition of the scalar dissipationrate.1286.2 Theorytion kinetics. These methods typically involve the use of quasi steady-state andpartial-equilibrium assumptions [32, 125] to reduce the complexity of the sys-tems. However, low-dimensional manifold approaches, first proposed by Maasand Pope [104], are particularly suited for complex simulations as they retain thefeatures of detailed mechanisms while significantly reducing the overall computa-tional cost. For this study, we have used a manifold approach referred to as theFlamelet-Generated Manifold (FGM) [121]. The particular manifold we are usingis tabulated as functions of two progress variables, the normalised CO2 and H2Omass fractions; it is thus referred to as a two-dimensional FGM, or 2D-FGM. The2D-FGM approach has been successfully coupled with CSE in the past [153].For premixed combustion involving short hydrocarbon chains, it is generallyacceptable to use a one-dimensional, unstrained laminar flame to construct themodel manifold. Each one-dimensional flame is provided with a realisable unburntfuel-air mixture; the solution to the flame is tabulated in the manifold with the endpoint represented by the combustion products. The FGM in this study is tabulatedas functions of two progress variables, the normalised CO2 and H2O mass frac-tions. During the LES calculation, the mass fractions and production rates of allmajor species can then be retrieved from the FGM if the two progress variables areprovided. It should be mentioned that FGM would typically encounter difficultiesin predicting slow chemical species, such as NOx and, to a lesser degree, CO. Thisis caused by the observation that the true, instantaneous reaction rates of these slowspecies are not always close to the tabulated manifold. As the primary objectivein this study is to determine the overall performance and differences between thenew ensemble selection and inversion algorithms, the FGM model is suitable forthe task.6.2.3 Presumed PDF and SDR modelsA number of different presumed PDF models have been previously investigated forpremixed combustion. The often used β -PDF does recover the extreme propertiesexpected of the true PDF, such as δ -functions at the zero and unity extremes ofreaction progress for maximal variance and single δ -functions at the mean for zerovariance. However, it fails to reproduce the shape of the true PDF in more general1296.2 Theorycases of premixed flames [70]. The issue is related to the fact that the shape ofthe true PDF appears to be a function of how the chemical reaction rates vary as afunction of the progress variable; hence, different chemical kinetics lead to differ-ent shapes of the PDF of progress variable. To account for the effects of chemistryon the shape of the PDF, Bray et al. [22] proposed using a premixed laminar flameto model the functional dependence of the PDF on progress variable. Their orig-inal formulation only provided coverage for flames with very high variance. Jinet al. [70] later proposed a modification of the Bray PDF that extends the originalformulation to cover all possible mean and variance combinations. This is accom-plished by truncating the PDF shape function as needed to match the mean andvariance parameters. It was found to significantly improve the fit of the PDF tothat extracted from DNS results. This PDF model by Jin et al., referred to as themodified laminar flamelet PDF (MLF-PDF), is commonly implemented in con-junction with the CSE combustion model [70, 143–145, 153]. A shortcoming ofthis method is that, at the point of truncation, the model PDF has a sharp drop tozero, whilst the true PDF tends to be more rounded.To address this issue, a one-dimensional turbulent method was proposed to takethe place of the typical laminar flame calculation to tabulate pseudo-turbulent PDFmodels for RANS and LES closures [175]. The Linear-Eddy Model [148, 157], aninexpensive one dimensional stochastic mixing model, has demonstrated the abilityto capture important effects from the interaction between chemistry and turbulenceon the PDF distributions sufficiently well [175]. This PDF model, formulated bythe pseudo-turbulent LEM calculation, is implemented in place of the conventionalmodified laminar flamelet PDF for the LES calculations in the current chapter.A related quantity is the scalar dissipation rate, χc, which is an important pa-rameter to both premixed and non-premixed combustion modelling [45, 181]. Forpremixed combustion, the SDR, χc, of a species-based reaction progress variableis,χc(T,φ) = Dc(T,φ)|∇c| · |∇c|, (6.3)where Dc(T,φ) and ∇c are the diffusivity of the selected species at the local tem-perature and equivalence ratio and the gradient of the progress variable, respec-tively. The SDR appears as an unclosed term in the transport equation for the vari-1306.2 Theoryance of progress variable, where χc directly measures the decay rate of fluctuationsthrough turbulent micromixing [181]. Since the burning rate of many combustionprocesses depends on the contact area and local gradient between the reactants, itis reasonable for most combustion models to assume that the mean burning rate ofthe flames either explicitly or implicitly depends on the scalar dissipation rate. Forexample, Conditional Moment Closure uses the scalar dissipation rate conditionedby the progress variable to calculate micromixing [181]; not surprisingly, mod-elling the conditional scalar dissipation term conditioned on the local and globalprogress of reaction, χc|c∗, emerges as one of the main difficulties in applyingCMC to turbulent premixed flames [4].In common practice, the SDR can be estimated through additional algebraicmodelling [40] or by including an extra transport equation related to the flame sur-face density as proposed in the Coherent Flame Model (CFM) [132]. However, afundamental relationship can be established between the conditional scalar dissi-pation rate (χc|c∗) and the PDF of the reaction progress variable (P(c∗)) [4],χc =∫ 10(χc|c∗)P(c∗) dc∗, (6.4)where χc is the unconditional mean scalar dissipation rate, or simply the meanscalar dissipation rate. Through the Linear-Eddy Model, it is possible to obtainpseudo-turbulent flame data that can directly compute both the χc|c∗ and the P(c∗),leading to the quantity required to close the variance transport equation, χc. Thissingle method can create both the PDF and SDR submodels, which are consistentwith each other. Furthermore, this formulation of the SDR model has the addedbenefit of pre-tabulation, which would reduce the computational time for the LESor RANS calculations. A detailed description of the methodology used in con-structing the PDF and SDR submodels can be found in [175]. This LEM SDRmodel is implemented in place of the conventional algebraic and transport equa-tion methods in the current chapter.1316.3 Numerical Methods / Simulation Setup6.3 Numerical Methods / Simulation Setup6.3.1 Pre-ProcessThe 2D-FGM, PDF and SDR submodels require pre-tabulation prior to the largeeddy simulation of the reactive flow. The temperature, pressure and equivalenceratio corresponding to the experimental condition for the methane-air burner arerespectively, 300 K, 101.325 kPa and 1.0. All of the submodels are specific tothe chemistry and thermodynamics conditions; thus, the temperature, pressure andequivalence ratio adhere to the quantities used in the experiments. The 2D-FGM isconstructed with the open source chemical kinetics and thermodynamics software,Cantera [58], using the GRI-MECH 3.0 mechanism [156]. A detailed descriptionof the 2D-FGM chemistry reduction technique can be found here [142]. Approxi-mately five to ten hours are required to construct each 2D-FGM table with a singlei7-2600K processor; however, the procedure can be readily parallelised if shorterconvergence time is desired as the unstrained laminar flame solutions with differentrealisable initial conditions are independent of one another.In addition to the thermo-chemical composition, the current PDF and SDRsubmodels are dependent on the turbulent intensity. As such, the overall turbulentfluctuation is approximately matched to the experimental condition (1.89 m/s) inthe LEM calculations. Following the procedure from previous work [175], approx-imately 10,000 pseudo-turbulent flame profiles are collected at 1/100 large eddyturnover time intervals at the specified equivalence ratio. These profiles are subse-quently used to construct the converged PDF and SDR submodels implemented asthe input to the CSE combustion model and the variance transport equation. Thecomputational time required to accumulate the 10,000 flame profiles and construct-ing the PDF and SDR submodels is approximately 48 hours on a single i7-2600Kprocessor using a six-step global mechanism.6.3.2 Non-Reactive SimulationA premixed, stoichiometric methane-air flame in the thin reaction zone has beensimulated. Similar to the experiments, the turbulent fluctuation (u′) is first de-termined via the cold flow. This is possible as it has been found that the total1326.3 Numerical Methods / Simulation Setupturbulence intensity value along the axial direction on the centerline of this burnerremains relatively constant between reacting and non-reacting conditions [33].For the experimental burner, a turbulence grid is fitted a distance, hg, upstreamof the nozzle. The value of hg for the current flame is 44.5 mm. Figure 6.1 illus-trates the schematic of the axisymmetric burner nozzle, while a detailed descrip-tion of the burner assembly can be found in [191]. It should be mentioned thatthe pilot has been turned off during the experimental cold flow measurements; thesimulations also adhere to this implemention for consistency. The relevant param-eters characterising the experimental conditions in this study are summarised inTable 6.1.11.2 mm0.889 mm1.397 mmhgPremixedGasPilot PilotSyntheticturbulence inletFigure 6.1: Burner cross-section.UB (m/s) u′ (m/s) ΛL (mm) Ret Ka Da21.0 1.89 2.64 316 2.6 9.6Table 6.1: Summary of experimental conditions. Symbols: UB = bulk flowvelocity in the central jet; u′ = total turbulence intensity; ΛL = longitudi-nal integral length scale; Ret = turbulent Reynolds number; Ka = turbulentKarlovitz number; Da = turbulent Damko¨hler number. The values are re-trieved from [173].To simulate the grid generated turbulence achieved in the experiments, a spec-tral inlet condition assuming isotropic homogeneous turbulence using the energy1336.3 Numerical Methods / Simulation Setupspectrum of Haworth and Poinsot [65] has been implemented. The time correlationfor the velocity fluctuation is based on the method developed by Billson [15]. Thetime correlation parameters, denoted a and b, represent the fraction of the currentvelocity field, a, to be added to the newly generated synthetic fluctuations, b. Theyare computed by the following relations: a= exp(−∆t/τt) and b=√1−a2, where∆t and τt are respectively the timestep size and the integral time scale. The param-eters for the spectral inlet condition required to generate the turbulent intensity forthe cold flow observed in the experiment are determined on the basis of trial anderror. The values are summarised in Table 6.2.Following the procedure of the experiments, the reported value of u′ (1.89 m/s)is averaged within the cylindrical volume of−0.25< r/D< 0.25 and 0.2< h/D<0.5, where r and h are respectively the radius and the axial distance downstream ofthe burner nozzle. The overall simulation domain length, diameter and hexahedralcell count for the cold flow are 100 mm, 60 mm and 2.25 million cells, respec-tively. The lower boundary of the simulation domain is aligned with location ofthe turbulence grid.hg (mm) K (m2/s2) Max. κ (m−1) Λ (mm) a b ∆t (s)44.5 570 2.10 · 104 2.64 0.985 0.171 2 · 10−6Table 6.2: Relevant synthetic turbulence parameters. Symbols: hg = gridlocation; K = turbulent kinetic energy; Λ = integral length; a and b = temporalcorrelation parameters; ∆t = timestep size.6.3.3 Reactive SimulationWith the appropriate turbulent parameters determined from the cold flow simula-tion, the reactive flow is initialised with premixed reactants in the central jet andstoichiometric methane-air combustion products in the pilot stream. The central jetbulk flow velocity is as described in Table 6.1, while the pilot bulk flow velocityis prescribed a magnitude of 1.0 m/s at the lower boundary of the concentric pilotpipe (Figure 6.1). The overall simulation domain length, diameter and hexahedralcell count for the reactive flow are 178.9 mm, 60 mm and 3.16 million cells, re-spectively. The simulation is performed with 96 2.2 GHz processors. Instantaneous1346.4 Results and Discussiontemperature profiles are collected after ten characteristic flame height flow-throughtimes from initialisation for 30 ms at 0.2 ms intervals.6.4 Results and DiscussionThe temperature field for the current flame is measured via Mie scattering3, wherebythe progress of reaction from the experiments is computed via binarised imagesdepending on light intensity [173]. The light intensity binarisation threshold dueto the evaporation of the oil particles corresponds to approximately 700 K; theprogress of reaction is assigned a value of 0 for temperatures below 700 K and 1otherwise. These binarised images are subsequently averaged and the iso-contourof 0.5 is used to determine the characteristic flame height; this approach has alsobeen adopted for the analysis of the LES temperature field.From Figure 6.2, it can be observed that the overall flame height and temper-ature distribution are well matched between the simulation and the experiment.The slight difference in the width of the flame could be caused by the syntheticturbulence inlet condition at the location of the turbulence grid. Improvementsto the central jet velocity profile at the synthetic turbulence inlet can be made, asthe velocity profile for the experimental burner exhibits a parabolic shape which isnot fully captured by the current implementation of the boundary condition. Theparabolic shape leads to a slightly slower transition of the progress variable fromzero to one along the centerline of the jet and faster transition of the progress vari-able from zero to one along the outer rim of the jet, observed in the right half ofFigure 6.2. The difference in overall behavior of the flame, however, may not besignifcant as the turbulence intensity has been carefully calibrated from the coldflow simulation.The accurate prediction of the overall flame distribution and height can be at-tributed to a number of improvements made in the current LES calculation in com-parison with previous CSE implementations for this laboratory burner, which op-erated in slightly different turbulent conditions [153]. In terms of the CFD module,the synthetic turbulent boundary condition is applied at the base of the premixed3The experimental data is provided by Mr. P. Tamadonfar and Dr. O¨. Gu¨lder of the University ofToronto [173].1356.4 Results and DiscussionFigure 6.2: Average binarised temperature-based progress variable profileand characteristic flame height from the LES result (left) and experimentalmeasurements (right), where D = 11.2 mm. The solid line represents thehalf-burning surface, at which the average progress variable takes on a valueof 0.5.gas pipe, which is located 44.5 mm from the burner nozzle (Figure 6.1); the in-clusion of the premixed gas pipe in the simulation allows the turbulent flow fieldto develop under more realistic conditions before interacting with the premixedflame front. Additionally, the synthetic turbulent intensity is carefully calibratedin the non-reactive simulation to match the experimental measurement prior to thereactive calculation.1366.4 Results and DiscussionFor the CSE module, we have implemented a recently developed dynamic en-semble selection algorithm [177] and a parallel iterative solver in a detailed LEScalculation for the first time. The dynamic ensemble selection algorithm increasesthe flexibility of the CSE module as the physical spaces in which the conditionalaverages are computed are automatically adjusted to the changing flame positionsthroughout the simulation. Such a feature ensures that an optimal number of re-active cells is provided for the inversion of Equation 5.4 for all of the ensembles,leading to increased load balance and more accurately inverted conditional massfractions. Furthermore, this automatic procedure reduces the complexity of im-plementation from the perspective of the user in two notable ways: i) substantialfore-knowledge of the flame distribution and flow is no longer required as the en-sembles can adapt to these changes as the simulation progresses; and ii) the numberof processors used in the LES computation can now be optimised only accordingto the CFD module, without having to account for the manual division of physicalspaces associated with the ensembles in the CSE module. In addition, the imple-mentation of the parallel iterative solver, parallel-LSQR, contributes to a reductionin computational time by approximately 50%, while maintaining accuracy in theinverse solution as compared to conventional CSE approaches [177].The LEM-formulated PDF model, which has been rigorously compared to ex-perimentally measured PDFs in a rudimentary sense [176], is implemented forthe first time in a LES of a turbulent premixed flame. This LEM PDF model hasdemonstrated improved accuracy compared to its predecessor, the MLF-PDF [175],which could account for part of the improvement in the prediction of the flame dis-tribution and height observed in this study. Similarly, the LEM-formulated SDRmodel is implemented for the first time in a LES calculation; it has been shownthat the LEM SDR matches experimental results satisfactorily for non-swirlingflames [176]. The SDR term is embedded in the transport equation of the vari-ance of the progress variable; improvements to this parameter can lead to a moreaccurate representation of the PDF model in each control volume.The computational time required to achieve the current result is documentedin Table 6.3. However, the times are expected to vary depending on the machineand communication characteristics between nodes. In general, one can anticipate aconverged solution for similar Bunsen-type burners using the simulation strategies1376.5 Concluding Remarksoutlined in this work within 200–300 hours of computation using 96 2.2 GHz pro-cessors. Indeed, increasing the number of processors would also increase the rateof convergence, though scalability has not been thoroughly studied.Processors Cell count Simulation time (s) Run time (hr)96 × 2.2-GHz 3.16 · 106 0.030 225Table 6.3: Summary of simulation times.6.5 Concluding RemarksA turbulent, premixed, stoichiometric methane-air laboratory Bunsen flame oper-ating in the thin reaction zone has been simulated with a recently developed versionof the Conditional Source-term Estimation combustion model designed for parallelcomputational environments. This version of CSE, in conjunction with newly for-mulated pseudo-turbulent PDF and SDR models, are implemented concurrently ina LES calculation for the first time. Comparison with experimental data suggeststhis combination of models can accurately predict the flame height and overall tem-perature distribution for our test case. Moreover, this implementation of CSE canproduce converged solutions with rather modest computational effort.Each of the implemented models, the CSE, 2D-FGM, LEM PDF and LEMSDR, is self-contained; in principle, it is possible to simulate any other turbulentpremixed flame with the same thermo-chemical composition and similar turbulentintensity without the need to readjust any model parameters. The approach outlinedin this work allows for progression towards a fully dynamic, predictive CFD solverfor turbulent premixed combustion.138Chapter 7Conclusion and Future Work7.1 ConclusionsThe objective of this thesis is to advance numerical simulation strategies for tur-bulent premixed combustion in terms of accuracy, autonomy and efficiency. Thefirst half of this work is focused on the development of new presumed PDF andSDR models formulated by the Linear-Eddy Model for the purposes of subsequentLES and RANS simulations. The LEM-based PDF and SDR models, tabulated asa function of the mean and variance of the reaction progress variable are able toaccount for the basic features of turbulence-chemistry interactions, a characteris-tic which all of the previous models lack. These new models are then comparedto experiments in detail, through a range of turbulent intensities and swirl condi-tions for validation. The second half of this work is focused on the optimisationand parallelisation of a combustion model for LES and RANS, the ConditionalSource-term Estimation. The methodologies developed herein can be straightfor-wardly extended to partially-premixed and non-premixed combustion. As the finalstudy of this thesis, the new LEM-based PDF and SDR models are combined withthe redesigned CSE module in a large eddy simulation of a premixed, methane-airlaboratory flame operating in the thin reaction zone.This work began with the observation that a single PDF model, the β -PDF,can be applied to non-premixed combustion model based on conditional momentapproaches; however, there does not appear to be a strong consensus on the most1397.1 Conclusionssuitable PDF model for turbulent premixed combustion simulation strategies. Theβ -PDF, Bray PDF and modified laminar flamelet PDF have been previously studiedin the context of premixed, methane-air flames; each of these models has demon-strated some form of deficiency in capturing the correct PDF shape. For the pur-pose of CSE calculations, it is critical that the PDF distribution is modelled asaccurately as computationally allowable, since the PDF is implemented twice foreach reaction rate calculation — once to obtain the conditional mass fractions andonce to obtain the source terms — at every timestep.Presented with this challenge, it would be appropriate to develop a compu-tationally inexpensive methodology for the construction of a more accurate PDFmodel for premixed combustion. It has been found that the Linear-Eddy Model, aone-dimensional stochastic mixing model, can reproduce turbulent statistics withrelatively low computational cost. As a result, this model is selected as a suitablecandidate to construct a PDF model that includes turbulence in its formulation.The one-dimensional premixed flame simulation code, including a suitable Linear-Eddy Model stochastic mixing module, is developed in this work. The resultantLEM-based PDFs show greater resemblance to the DNS and experimental mea-surements in comparison with existing presumed PDF models. Concurrently, theLEM data used to construct the PDF model can be applied to formulate a pseudo-turbulent conditional scalar dissipation rate model. The combination of these twomodels would allow for the implementation of consistent LEM-based PDF andSDR models in the subsequent LES calculations. It is shown that the LEM SDR isable to predict the decrease in the magnitude of these distributions caused by flamethickening due to turbulent effects in premixed combustion, albeit over-predictingthe experimental results in peak magnitudes.The PDF and SDR models formulated in this work are suitable for a numberof simulation strategies for turbulent premixed combustion, including PresumedConditional Moment, Conditional Moment Closure and Conditional Source-termEstimation. The LEM-based SDR model is particularly useful for the CMC ap-proach as the lack of accurate SDR models is limiting the general applicability ofCMC in this combustion mode. The successful development of a one-dimensionalpseudo-turbulent flame code, LEM PDF model and LEM SDR model motivatedthe subsequent topic of research, which is to extend the current study from a purely1407.1 Conclusionspractical basis — that is, to simply provide better submodels for advanced combus-tion models — to a theoretical investigation, where perhaps it would be possibleto understand the underlying characteristics of the PDF and SDR distributions as afunction of turbulent fluctuations and swirling conditions.The Cambridge swirl burner, an experimental bluff body stabilised burner, isselected as the basis of validation of the new LEM-based PDF and SDR models. Abroad range of turbulent conditions, including changes in intensities and swirl num-bers, are simulated and compared with the experimental results. The algorithms forthe construction of the LEM PDF and LEM SDR models are also parallelised andoptimised at this stage.The LEM PDF model demonstrates good agreement with the experimental re-sults in terms of nominal changes with turbulence. Particularly, it is able to capturethe overall shape and the effect of smoothing towards the left and right bound-aries of the PDF distributions — a feature that cannot be reproduced by any ofthe previous PDF models. Further, the LEM results suggest that the conditionalSDR would decrease towards the upper regions of the Borghi diagram and wouldincrease in peak magnitudes towards the unstrained laminar limit with increasingintegral scales. The experimental conditional SDRs for all axial locations of theflame brush show a greater decrease in magnitude than predicted by the numericalmodel, though this could be partially caused by an interaction between shearing andmixing layers within a specific region of the swirl flames, leading to greater dissi-pation downstream of this location. The values for the unconditional mean SDRpredicted by LEM are typically within 50% of the experimental values. Moreover,the accuracy of the numerical model seem to be dependent on the swirl number ofthe flow.Overall, the results suggest that the PDF tabulated as a function of the meanand variance of the reaction progress variable for premixed flames remains rela-tively steady under a variety of turbulent conditions, including changes in the in-tegral length, velocity fluctuation and swirl number. From a practical perspective,this finding is remarkable as it implies that one can use a representative pseudo-turbulent PDF model for a range of turbulent conditions in subsequent RANS orLES calculations, which greatly simplifies the complexity of implementation withthe corresponding combustion models. The validation with the Cambridge swirl1417.1 Conclusionsburner data firmly establishes that LEM can be used to formulate a high-fidelityPDF model for turbulent premixed combustion.Swirl appears to have greater influence on the local flame gradients; this, inturn, leads to more substantial changes in the SDR, which the LEM cannot fullycapture in its current form. However, the magnitude of the SDR as calculated byLEM is, in the most extreme cases, over-predicted by a factor of two in comparisonwith the experimental flames in this study. Moreover, the shape of the conditionalSDR distribution is well-captured with the LEM methodology. These character-istics imply: (i) that it could be possible to apply an empirical correction for theconditional SDR as a simple function of the swirl number and turbulent fluctuation;and (ii) the LEM-formulated SDR is a competitive candidate amongst the availableSDR models for turbulent premixed combustion modelling.As a result of this study, an accurate, pseudo-turbulent PDF model suitable forsubsequent RANS and LES calculations has been formulated. Further, a compan-ion SDR model that is consistent with the LEM PDF has been also been estab-lished. It is hoped that these new models would not only serve to demonstrate thetheoretical findings described above, but would also be useful to a number of thesimulation strategies for turbulent premixed combustion. For the reason that thesemodels are pre-tabulated prior to the detailed flame simulation, they are ideallysuited for large-scale, complex problems by focusing the computational resourcesto solving the fluid dynamics calculation rather than the submodel parameters.The next part of this thesis work is dedicated to the parallelisation and opti-misation of the Conditional Source-term Estimation combustion model for com-putations in large clusters. Prior to this research, the CSE model had severalimplementational obstacles that limited its performance in largely parallel com-putational environments; these include the pre-allocation of static ensembles andsingle-processor inversion calculations. In reference to the former, it would bedifficult for the user to achieve load balance during the simulation as the flame dis-tribution evolves dynamically within the domain. Further, the exact division of theensembles is entirely dependent on the foreknowledge of the flame morphologyby the user; this leads to the question of completeness of the combustion modelas different ensemble divisions could potentially lead to different solutions. An-other problem is that the number of reactive cells may not be fully satisfied if the1427.1 Conclusionsensembles are not carefully divided according to the available processors.To address these shortcomings, two optimisations to enhance the capabilityof CSE are introduced. The first optimisation is to increase the flexibility of theCSE model by dynamically allocating the physical space required for the individ-ual ensemble averages. The second optimisation is to integrate a well-establishediterative solver, Least-Squares with QR-factorisation, to tackle the CSE inversionproblem. As the purpose of this work is to prepare the CSE module for cluster op-eration, the LSQR solver is also parallelised to ensure scalability of the combustioncode. These optimisations can be applied to both RANS and LES paradigms in thepremixed, partially-premixed and non-premixed combustion modes.The new implementations are subsequently validated with several generic LEScalculations of a turbulent, premixed methane-air Bunsen burner. The standaloneinversion tests and LES results indicate that the parallelised LSQR solver is fullycapable of reproducing similar solutions to the conventional CSE methods whilesubstantially reducing computational time. The data sharing algorithm betweenensembles also proves to be robust in achieving the proper communication direc-tives for three sets of grids, three different processor counts and three different setsof computational hardware. The integration of the dynamic ensemble selectionand parallel LSQR algorithms allows one to establish the framework for a highlyautonomous and parallel version of CSE that can adapt to the changing flame dis-tribution in complex geometries.As combustion modelling for turbulent premixed flames is currently a subjectof significant research because of the intricate interactions between turbulence andchemistry, the newly optimised version of CSE demonstrates itself to be a highlypractical approach in several notable ways. First, the computational times havebeen substantially reduced in comparison with previous implementations, by ap-proximately 50% in most LES cases. Second, this variant can be considered morecomplete as the need for manual adjustment of the ensembles has been eliminated.Third, the CSE combustion model is more accessible to the average user as the en-sembles are automatically constructed to best suit the flame position and reactivitywithin the simulation domain at any given time step — substantial foreknowledgeof the flame morphology is no longer required for the division of the ensembles.Another benefit to the optimisations described herein is that the CSE model is1437.1 Conclusionsdesigned for RANS or LES calculations under all ranges of air-fuel ratios and tur-bulent intensities; this implies that the algorithms developed is applicable to all ofthe combustion modes within the framework of both RANS and LES. In conse-quence, the current variant of CSE is not only more suitable for academic purposesof validation, but is also more suitable for industrial problems from a practicalperspective.In the final study of this thesis work, a detailed LES calculation of a pre-mixed methane-air laboratory-scale burner operating in the thin reaction zone isperformed; this simulation neatly connects the LEM-formulated PDF and SDRmodels with the redesigned CSE combustion module. Comparison with exper-imental data, measured by Mie scattering, suggests this combination of modelscan accurately predict the flame height and overall distribution. Moreover, thisimplementation of CSE can produce converged solutions with rather modest com-putational effort. Each of the implemented models, the CSE, FGM, LEM PDFand LEM SDR, is self-contained; in principle, it is possible to simulate any otherturbulent premixed flame with the same thermo-chemical composition and similarturbulent intensity without the need to adjust any modelling parameters.This LES of the laboratory-scale burner is an appropriate terminus for the cur-rent thesis as it serves to demonstrate the accomplishments of all of the previouschapters. The pseudo-turbulent PDF and SDR models presented in Chapters 3and 4 are fully integrated into the optimised CSE combustion model of Chapter 5for a detailed simulation. Each chapter represents an important advancement in themodelling of turbulent premixed combustion, in terms of accuracy and efficiency.The summation of approaches outlined herein allows for progression towards afully dynamic, predictive CFD solver for turbulent premixed combustion devoid ofany tunable parameters. With continued research, it is hoped that these methodscan eventually be implemented in industrially relevant LES calculations to accu-rately predict the formation of minor species and the thermodynamic characteris-tics of the next generation of combustion devices. Futhermore, perhaps the outlinedapproaches can be refined in the future to a level where they can provide a morefundamental understanding of the nature of the turbulence-chemistry interactionwith high-fidelity large eddy simulations.1447.2 Future Work7.2 Future WorkThe first half of this research is concerned with the development of more accuratePDF and SDR models for turbulent premixed combustion. While it appears that theLEM PDF can replicate the measured PDF to a high degree of precision for a broadrange of turbulence conditions, the same cannot be said of the LEM SDR model. Itwould be of importance to understand the underlying mechanisms that decrease themagnitude of the SDR with increasing swirl intensity. Further, it may be possibleto formulate an empirical correlation between the magnitude of the SDR and thevarious turbulent characteristics, such as swirl and turbulent fluctuations.The second half of this research pertains to further developing the ConditionalSource-term Estimation combustion model for turbulent premixed flames. Havingoptimised a version of CSE that is capable of predicting basic flame attributesfrom a Bunsen-type laboratory burner with very modest computational resources,it would be interesting to examine other burners, fuels and combustion modes.A particularly attractive set of experiments for simulation purposes is producedby the Cambridge swirl burner [167]. The temperature and species mass fractionsare measured accurately and well-documented [168, 170–172, 194]. Further, theburner geometry and boundary conditions are well-characterised and suitable forsimulations. It would be interesting to use the new formulation of CSE to simulatethis burner and to study the swirling effects. The ability to simulate swirl burn-ers is fundamental to combustion models designed for industrial applications asmodelling gas turbine combustors is one of their most important functions.Another direction of research would be to test different fuels with the CSEcombustion model. This study has focused on methane-air combustion as it is thesimplest of hydrocarbon fuels. More complex hydrocarbons and alternative fuelsources, such as biofuels, could be interesting candidates due to their increasingpopularity as fuels for a new generation of engines and slower chemical time scales.The extended chemical times allows for greater turbulent mixing within the reactiveflow, leading to more substantial chemistry-turbulence interactions. This extendedperiod allows for more thorough validation of the combustion model, provideddetailed experimental data is available. In addition, perhaps it would be feasible toimplement different chemistry reduction techniques in conjunction with CSE, such1457.2 Future Workas skeletal mechanisms or other reduced manifold approaches for such complexfuels.In reality, combustion rarely occurs in a purely premixed or non-premixedmode; it is, thus, imperative for combustion models to perform equally well inthe partially-premixed mode. The primary challenge for CSE here is that the inver-sion calculations are substantially more expensive. However, the newly developeddynamic ensemble selection algorithms and the newly integrated parallel iterativesolver from this thesis work should alleviate a considerable part of the computa-tional load. Thus, another possible research direction is the optimisation of thePLSQR and ensemble selection algorithms for the partially-premixed combustionmode in the context of CSE.146BibliographyBibliography[1] Large eddy simulation of combustion instabilities in turbulent premixedburners, 1997. SAND86-8841. → pages 39[2] Key World Energy Statistics 2015. https://www.iea.org/publications/freepublications/publication/key-world-energy-statistics-2015.html, 2015.Accessed: 2016-01-14. → pages 1[3] R. Abdel-Gayed, D. Bradley, and F.-K. Lung. Combustion regimes and thestraining of turbulent premixed flames. Combustion and Flame, 76(2):213–218, 1989. → pages 23[4] 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. → pages 41, 58, 68, 86, 131[5] G. K. Bachelor. The Theory of Homogeneous Turbulence. CambridgeUniversity Press, 1982. → pages 19[6] B. S. Baldwin and H. Lomax. Thin-layer approximation and algebraicmodel for separated turbulent flows. In Paper 91-0610, AIAA, 1978. →pages 30[7] J. Bardina, J. Ferziger, and W. Reynolds. chapter Improved subgrid-scalemodels for large eddy simulation. Fluid Dynamics and Co-locatedConferences. American Institute of Aeronautics and Astronautics, Jul1980. → pages 33[8] R. Barlow, G.-H. Wang, P. Anselmo-Filho, M. Sweeney, and S. Hochgreb.Application of Raman/Rayleigh/LIF diagnostics in turbulent stratifiedflames. Proceedings of the Combustion Institute, 32(1):945–953, 2009. →pages 59147Bibliography[9] F. Bastin, P. Lafon, and S. Candel. Computation of jet mixing noise due tocoherent structures: the plane jet case. Journal of Fluid Mechanics, 335:261–304, 3 1997. → pages 30[10] A. Belcadi, M. Assou, E. Affad, and E. Chatri. Construction of a reducedmechanism for modelling premixed combustion of methane-air.Combustion Theory and Modelling, 11(4):603–613, 2007. → pages 14[11] A. Belcadi, M. Assou, E. Affad, and E. Chatri. CH4/NOx reducedmechanisms used for modeling premixed combustion. Energy and PowerEngineering, 4(4):264–273, 2012. → pages 14[12] R. Bilger, S. Pope, K. Bray, and J. Driscoll. Paradigms in turbulentcombustion research. Proceedings of the Combustion Institute, 30(1):21–42, 2005. → pages 42[13] R. W. Bilger. Conditional Moment Closure for turbulent reacting flow.Physics of Fluids A, 5(2):436–444, 1993. → pages 41[14] R. W. Bilger, S. B. Pope, K. N. C. Bray, and J. F. Driscoll. Paradigms inturbulent combustion research. Proceedings of the Combustion Institute,30:21–42, 2005. → pages 38[15] M. Billson, L.-E. Eriksson, and L. Davidson. Jet noise prediction usingstochastic turbulence modeling. In AIAA paper 2003-3282, 9thAIAA/CEAS Aeroacoustics Conference, 2003. → pages x, 115, 116, 134[16] R. B. Bird, W. E. Stewart, and E. N. Lightfoot. Transport phenomena.AIChE Journal, 7(2):5J–6J, 1961. → pages 9[17] M. Boger, D. Veynante, H. Boughanem, and A. Trouve´. Direct numericalsimulation analysis of flame surface density concept for large eddysimulation of turbulent premixed combustion. Symposium (International)on Combustion, 27(1):917–925, 1998. → pages 37[18] M. Boileau, G. Staffelbach, B. Cuenot, T. Poinsot, and C. Be´rat. LES of anignition sequence in a gas turbine engine. Combustion and Flame, 154(12):2–22, 2008. → pages 39[19] R. Borghi. Recent Advances in the Aerospace Sciences: In Honor of LuigiCrocco on His Seventy-fifth Birthday, chapter On the Structure andMorphology of Turbulent Premixed Flames, pages 117–138. Springer US,Boston, MA, 1985. ISBN 978-1-4684-4298-4. → pages 23148Bibliography[20] K. Bray and J. Moss. A unified statistical model of the premixed turbulentflame. Acta Astronautica, 4(34):291–319, 1977. → pages 35, 36[21] K. Bray, P. A. Libby, and J. Moss. Unified modeling approach for premixedturbulent combustion–part I: General formulation. Combustion and Flame,61(1):87–102, 1985. → pages 35[22] K. Bray, M. Champion, P. Libby, and N. Swaminathan. Finite ratechemistry and presumed PDF models for premixed turbulent combustion.Combustion and Flame, 146(4):665–673, 2006. → pages 45, 46, 66, 130[23] K. N. C. Bray. Turbulent Reacting Flows, chapter Turbulent flows withpremixed reactants, pages 115–183. Springer Berlin Heidelberg, Berlin,Heidelberg, 1980. ISBN 978-3-540-38273-7. → pages 23[24] K. N. C. Bray, M. Champion, and P. A. Libby. Turbulent Reactive Flows,chapter The Interaction Between Turbulence and Chemistry in PremixedTurbulent Flames, pages 541–563. Springer US, New York, NY, 1989.ISBN 978-1-4613-9631-4. → pages 35[25] A. Burcat. Thermochemical data for combustion calculations. InCombustion Chemistry, pages 455–473. Springer US, 1984. → pages 11[26] W. K. Bushe and H. Steiner. Conditional Moment Closure for large eddysimulation of non-premixed turbulent reacting flows. Physics of Fluids, 11(7):1896–1906, 1999. → pages 42, 43, 98, 100, 127[27] T. Butler and P. O’Rourke. A numerical method for two dimensionalunsteady reacting flows. Symposium (International) on Combustion, 16(1):1503–1515, 1977. → pages 38[28] S. Cannon, B. Brewster, and L. Smoot. PDF modeling of lean premixedcombustion using in situ tabulated chemistry. Combustion and Flame, 119(3):233–252, 1999. → pages 41[29] W. Chang and J. Chen. Reduced mechanisms for premixed andnon-premixed combustion, 1999. URLhttp://firebrand.me.berkeley.edu/griredu.html. Accessed: 2016-01-14. →pages 14, 51, 71[30] F. Charlette, C. Meneveau, and D. Veynante. A power-law flame wrinklingmodel for LES of premixed turbulent combustion part I: non-dynamicformulation and initial tests. Combustion and Flame, 131(12):159–180,2002. → pages 37149Bibliography[31] F. Charlette, C. Meneveau, and D. Veynante. A power-law flame wrinklingmodel for LES of premixed turbulent combustion part II: dynamicformulation. Combustion and Flame, 131(12):181–197, 2002. → pages 37[32] J. Y. Chen. A general procedure for constructing reduced reactionmechanisms with given independent relations. Combustion Science andTechnology, 57(1-3):89–94, 1988. → pages 13, 129[33] Y.-C. Chen, N. Peters, G. Schneemann, N. Wruck, U. Renz, and M. S.Mansour. The detailed flame structure of highly stretched turbulentpremixed methane-air flames. Combustion and Flame, 107(3):223–244,1996. → pages 133[34] T. J. Chung. Computational Fluid Dynamics. Cambridge University Press,2002. ISBN 0521594162. → pages 27[35] O. Colin, F. Ducros, D. Veynante, and T. Poinsot. A thickened flame modelfor large eddy simulations of turbulent premixed combustion. Physics ofFluids, 12(7):1843–1863, 2000. → pages 39[36] G. Damko¨hler. Der Einfluss der Turbulenz auf dieFlammengeschwindigkeit in Gasgemischen. Z. Electrochem, 46:601–652,1940. → pages 22, 23[37] J. W. Deardorff. A numerical study of three-dimensional turbulent channelflow at large Reynolds numbers. Journal of Fluid Mechanics, 41:453–480,4 1970. → pages 31[38] J. W. Deardorff. The use of subgrid transport equations in athree-dimensional model of atmospheric turbulence. ASME Transactions,Journal of Fluids Engineering, 95:429–438, 1973. → pages 33[39] P. Domingo, L. Vervisch, S. Payet, and R. Hauguel. DNS of a premixedturbulent V flame and LES of a ducted flame using a FSD-PDF subgridscale closure with FPI-tabulated chemistry. Combustion and Flame, 143(4):566–586, 2005. → pages 39, 44, 65, 74[40] P. Domingo, L. Vervisch, and D. Veynante. Large eddy simulation of alifted methane jet flame in a vitiated coflow. Combustion and Flame, 152(3):415–432, 2008. → pages 40, 131[41] J. J. Dongarra, I. S. Duff, D. C. Sorensen, and H. A. van der Vorst.Numerical Linear Algebra on High-Performance Computers. Society forIndustrial and Applied Mathematics, 1998. → pages 117150Bibliography[42] D. Dovizio, M. M. Salehi, and C. B. Devaud. RANS simulation of aturbulent premixed bluff body flame using Conditional Source-termEstimation. Combustion Theory and Modelling, 17(5):935–959, 2013. →pages 42, 98, 101, 127[43] D. Dovizio, J. W. Labahn, and C. B. Devaud. Doubly ConditionalSource-term Estimation (DCSE) applied to a series of lifted turbulent jetflames in cold air. Combustion and Flame, 162(5):1976–1986, 2015. →pages 42, 98, 116, 127[44] J. F. Driscoll. Turbulent premixed combustion: Flamelet structure and itseffect on turbulent burning velocities. Progress in Energy and CombustionScience, 34(1):91–134, 2008. → pages 38[45] J. Duclos, D. Veynante, and T. Poinsot. A comparison of flamelet modelsfor premixed turbulent combustion. Combustion and Flame, 95(1 - 2):101–117, 1993. → pages 58, 68, 69, 130[46] M. Du¨sing, A. Sadiki, and J. Janicka. Towards a classification of modelsfor the numerical simulation of premixed combustion based on ageneralized regime diagram. Combustion Theory and Modelling, 10(1):105–132, 2006. → pages 23[47] T. Echekki and E. Mastorakos. Turbulent Combustion Modeling:Advances, New Trends and Perspectives, volume 95. Springer Science,2011. → pages 37, 41[48] C. L. Fefferman. Existence and smoothness of the Navier-Stokes equation.URL http://www.claymath.org/millennium-problems/navier%E2%80%93stokes-equation. Accessed: 2016-01-15. → pages 18, 27[49] C. Fenimore. Formation of nitric oxide in premixed hydrocarbon flames.Symposium (International) on Combustion, 13(1):373–380, 1971. → pages66[50] C. FUREBY. Large eddy simulation of combustion instabilities in a jetengine afterburner model. Combustion Science and Technology, 161(1):213–243, 2000. → pages 37[51] C. Fureby. A fractal flame-wrinkling large eddy simulation model forpremixed turbulent combustion. Proceedings of the Combustion Institute,30(1):593–601, 2005. → pages 37151Bibliography[52] J. Galpin, A. Naudin, L. Vervisch, C. Angelberger, O. Colin, andP. Domingo. Large eddy simulation of a fuel-lean premixed turbulentswirl-burner. Combustion and Flame, 155(12):247–266, 2008. → pages 40[53] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamicsubgridscale eddy viscosity model. Physics of Fluids A, 3(7):1760–1765,1991. → pages 33[54] S. Ghosal and P. Moin. The basic equations for the large eddy simulationof turbulent flows in complex geometry. Journal of Computational Physics,118(1):24–37, 1995. → pages 32, 33[55] O. Gicquel, N. Darabiha, and D. The´venin. Laminar premixedhydrogen/air counterflow flame simulations using flame prolongation ofILDM with differential diffusion. Proceedings of the Combustion Institute,28(2):1901–1908, 2000. → pages 40, 44, 65, 74[56] P. Givi. Model-free simulations of turbulent reactive flows. Progress inEnergy and Combustion Science, 15(1):1–107, 1989. → pages 41[57] G. Golub and W. Kahan. Calculating the singular values andpseudo-inverse of a matrix. Journal of the Society for Industrial andApplied Mathematics Series B Numerical Analysis, 2(2):205–224, 1965. →pages 111[58] D. G. Goodwin, H. K. Moffat, and R. L. Speth. Cantera: Anobject-oriented software toolkit for chemical kinetics, thermodynamics,and transport processes. http://www.cantera.org, 2015. Version 2.2.0. →pages xiv, xv, 17, 51, 52, 58, 72, 132[59] S. Gordon and B. J. McBride. Computer program for calculation ofcomplex chemical equilibrium compositions and applications: I. analysis.Technical report, NASA Reference Publication 1311, 1994. → pages 10[60] C. J. Greenshields. Upstream Class Reference, 2011. URLhttp://foam.sourceforge.net/docs/cpp/a02750.html. Accessed: 2016-01-20.→ pages 107[61] F. Grinstein, L. Margolin, and W. Rider. Implicit Large Eddy Simulation:Computing Turbulent Fluid Dynamics. Cambridge University Press, 2007.ISBN 9780521869829. → pages 31[62] R. W. Grout. An age extended progress variable for conditioning reactionrates. Physics of Fluids, 19(10):105–107, 2007. → pages 53152Bibliography[63] P. C. Hansen. Numerical tools for analysis and solution of Fredholmintegral equations of the first kind. Inverse Problems, 8(6):849, 1992. →pages 101[64] E. Hawkes and R. Cant. A flame surface density approach to large eddysimulation of premixed turbulent combustion. Proceedings of theCombustion Institute, 28(1):51–58, 2000. → pages 36, 37[65] D. C. Haworth and T. J. Poinsot. Numerical simulations of Lewis numbereffects in turbulent premixed flames. Journal of Fluid Mechanics, 244:405–436, 11 1992. → pages 115, 134[66] H. Huang, J. M. Dennis, L. Wang, and P. Chen. A scalable parallel LSQRalgorithm for solving large-scale linear system for tomographic problems:A case study in seismic tomography. Procedia Computer Science, 18(0):581–590, 2013. → pages 112[67] J. Huang and W. K. Bushe. Simulation of transient turbulent methane jetignition and combustion under engine-relevant conditions usingConditional Source-term Estimation with detailed chemistry. CombustionTheory and Modelling, 11(6):977–1008, 2007. → pages 42, 98, 101, 127[68] Y. Huang, H.-G. Sung, S.-Y. Hsieh, and V. Yang. Large eddy simulation ofcombustion dynamics of lean-premixed swirl-stabilized combustor.Journal of Propulsion and Power, 19(5):782–794, Sep 2003. → pages 38[69] H. G. Im, T. S. Lund, and J. H. Ferziger. Large eddy simulation of turbulentfront propagation with dynamic subgrid models. Physics of Fluids, 9(12):3826–3833, 1997. → pages 38[70] B. Jin, R. Grout, and W. Bushe. Conditional Source-term Estimation as amethod for chemical closure in premixed turbulent reacting flow. Flow,Turbulence and Combustion, 81:563–582, 2008. → pages xiii, xiv, 40, 43,44, 45, 47, 49, 62, 65, 66, 68, 74, 79, 100, 101, 130[71] M. M. Kamal, R. Zhou, S. Balusamy, and S. Hochgreb. Favre- andReynolds-averaged velocity measurements: Interpreting PIV and LDAmeasurements in combustion. Proceedings of the Combustion Institute, 35(3):3803–3811, 2015. → pages x, 75[72] V. Katta and W. M. Roquemore. C/H atom ratio inrecirculation-zone-supported premixed and non-premixed flames.Proceedings of the Combustion Institute, 34(1):1101–1108, 2013. → pages69153Bibliography[73] R. J. Kee, G. Dixon-lewis, J. Warnatz, M. E. Coltrin, and J. A. Miller. Afortran computer code package for the evaluation of gas-phase,multicomponent transport properties. Technical report, 1986. → pages 8,10[74] R. J. Kee, F. M. Rupley, and J. A. Miller. CHEMKIN-II: A FORTRANchemical kinetics package for the analysis of Gas-Phase chemical kinetics.1989. → pages 8, 13, 51, 71[75] R. J. Kee, F. M. Rupley, and J. A. Miller. The CHEMKIN ThermodynamicData Base, 1991. → pages 10, 11, 71[76] R. J. Kee, F. M. Rupley, J. A. Miller, M. E. Coltrin, J. F. Grcar, E. Meeks,H. K. Moffat, A. E. Lutz, G. Dixon-Lewis, M. D. Smooke, J. Warnatz,G. H. Evans, R. S. Larson, R. E. Mitchell, L. R. Petzold, W. C. Reynolds,M. Caracotsios, W. E. Stewart, P. Glarborg, C. Wang, and O. Adigun.PREMIX Program, CHEMKIN Collection, Release 3.6, 2000. → pages 49[77] A. R. Kerstein. Linear-eddy modeling of turbulent transport. Part 2:Application to shear layer mixing. Combustion and Flame, 75(3-4):397–413, 1989. → pages xiv, 48, 49, 50, 51, 69[78] A. R. Kerstein. Linear-eddy modeling of turbulent transport. Part 3:Mixing and differential molecular diffusion in round jets. Journal of FluidMechanics, 216:411–435, 1990. → pages[79] A. R. Kerstein. Linear-eddy modeling of turbulent transport. Part 6:Microstructure of diffusive scalar mixing fields. Journal of FluidMechanics, 231:361–394, 1991. → pages 70[80] A. R. Kerstein. Linear-eddy modeling of turbulent transport. Part 4:Structure of diffusion flames. Combustion Science and Technology, 81(1-3):75–96, 1992. → pages[81] A. R. Kerstein. Linear-eddy modeling of turbulent transport. Part 7:Finite-rate chemistry and multi-stream mixing. Journal of FluidMechanics, 240:289–313, 1992. → pages 49, 69[82] S. H. Kim and H. Pitsch. Conditional filtering method for large eddysimulation of turbulent non-premixed combustion. Physics of Fluids, 17(10), 2005. → pages 41154Bibliography[83] W.-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. → pages 38[84] A. Klimenko and R. Bilger. Conditional Moment Closure for turbulentcombustion. Progress in Energy and Combustion Science, 25(6):595–687,1999. → pages 41, 98, 126[85] A. Y. Klimenko. Multicomponent diffusion of various admixtures inturbulent flow. Fluid Dynamics, 25(3):327–334, 1990. → pages 41[86] A. Y. Klimenko and S. B. Pope. The modeling of turbulent reactive flowsbased on multiple mapping conditioning. Physics of Fluids (1994-present),15(7):1907–1925, 2003. → pages 98, 126[87] R. Knikker, D. Veynante, and C. Meneveau. A priori testing of a similaritymodel for large eddy simulations of turbulent premixed combustion.Proceedings of the Combustion Institute, 29:2105–2111, 2002. → pages 37[88] E. Knudsen and H. Pitsch. A dynamic model for the turbulent burningvelocity for large eddy simulation of premixed combustion. Combustionand Flame, 154(4):740–760, 2008. → pages 38[89] A. N. Kolmogorov. The local structure of turbulence in incompressibleviscous fluid for very large reynolds numbers. Dokl. Akad. Nauk SSSR, 30:299–303, 1921. → pages 20, 31[90] P. S. Kothnur and N. T. Clemens. Effects of unsteady strain rate on scalardissipation structures in turbulent planar jets. Physics of Fluids, 17(12):–,2005. → pages 94[91] K. K. Kuo. Principles of Combustion. John Wiley & Sons, 1986. → pages14[92] J. Labahn, D. Dovizio, and C. Devaud. Numerical simulation of thedelft-jet-in-hot-coflow (DJHC) flame using Conditional Source-termEstimation. Proceedings of the Combustion Institute, 35(3):3547–3555,2015. → pages 42, 98, 101, 127[93] J. W. Labahn and C. B. Devaud. Investigation of Conditional Source-termEstimation applied to a non-premixed turbulent flame. Combustion Theoryand Modelling, 17(5):960–982, 2013. → pages 42, 98, 101, 127155Bibliography[94] B. E. Launder and D. B. Spalding. Mathematical Models of Turbulence.Academic Press, 2002. → pages 30[95] B. E. Launder, G. J. Reece, and W. Rodi. Progress in the development of aReynolds-stress turbulence closure. Journal of Fluid Mechanics, 68:537–566, 4 1975. → pages 30[96] C. K. Law. Combustion Physics. Cambridge University Press, 2006. →pages 14[97] A. Leonard. Energy cascade in large eddy simulations of turbulent fluidflows. In F. Frenkiel and R. Munn, editors, Turbulent Diffusion inEnvironmental Pollution Proceedings of a Symposium held atCharlottesville, volume 18, Part A of Advances in Geophysics, pages237–248. Elsevier, 1975. → pages 32[98] P. A. Libby and K. N. C. Bray. Variable density effects in premixedturbulent flames. AIAA Journal, 15(8):1186–1193, Aug 1977. → pages 35[99] P. A. Libby and F. A. Williams. Turbulent Reacting Flows. London :Academic, 1994. ISBN 0124479456. → pages 28[100] D. K. Lilly. A proposed modification of the Germano subgrid-scale closuremethod. Physics of Fluids A, 4(3), 1992. → pages 33[101] R. Lindstedt and E. Vaos. Transported PDF modeling ofhigh-Reynolds-number premixed turbulent flames. Combustion and Flame,145(3):495–511, 2006. → pages 41[102] T. Lu and C. K. Law. Toward accommodating realistic fuel chemistry inlarge-scale computations. Progress in Energy and Combustion Science, 35(2):192–215, 2009. → pages 13[103] T. Lund and E. A. Novikov. Parameterization of subgrid-scale stress by thevelocity gradient tensor. Center for Turbulence Research (StanfordUniversity and NASA), 1992. → pages 33[104] U. Maas and S. Pope. Simplifying chemical kinetics: Intrinsiclow-dimensional manifolds in composition space. Combustion and Flame,88(34):239–264, 1992. → pages 14, 129[105] B. Magnussen and B. Hjertager. On mathematical modeling of turbulentcombustion with special emphasis on soot formation and combustion.Symposium (International) on Combustion, 16(1):719–729, 1977. → pages35156Bibliography[106] T. Mantel and R. Borghi. A new model of premixed wrinkled flamepropagation based on a scalar dissipation equation. Combustion andFlame, 96(4):443–457, 1994. → pages 36[107] G. H. Markstein. Non-steady Flame Propagation. Pergamon Press, 1964.→ pages 38[108] S. Mathur, P. K. Tondon, and S. C. Saxena. Thermal conductivity of binary,ternary and quaternary mixtures of rare gases. Molecular Physics, 12:569–579, 1967. → pages 9[109] S. McAllister, J. Y. Chen, and A. C. Fernandez-Pello. Fundamentals ofCombustion Processes. Springer, 2011. → pages 1[110] B. J. McBride and S. Gordon. Computer program for calculation ofcomplex chemical equilibrium compositions and applications: II. usersmanual and program description. Technical report, NASA ReferencePublication 1311, 1996. → pages 10[111] C. Meneveau, T. S. Lund, and W. H. Cabot. A Lagrangian dynamicsubgrid-scale model of turbulence. Journal of Fluid Mechanics, 319:353–385, 7 1996. → pages 33[112] P. Moin and K. Mahesh. Direct Numerical Simulation: A tool in turbulenceresearch. Annual Review of Fluid Mechanics, 30(1):539–578, 1998. →pages 28[113] C. Montgomery, G. Kosa´ly, and J. Riley. Direct numerical solution ofturbulent non-premixed combustion with multistep hydrogen-oxygenkinetics. Combustion and Flame, 109(1-2):113–144, 1997. → pages 88[114] A. Msaad, A. Belcadi, M. Mahdaoui, E. Aaffad, and M. Mouqallid.Reduced detailed mechanism for methane combustion. Energy and PowerEngineering, 4(1):28–33, 2012. → pages 14[115] J. R. Nanduri, D. R. Parsons, S. L. Yilmaz, I. B. Celik, and P. A. Strakey.Assessment of RANS-based turbulent combustion models for prediction ofemissions from lean premixed combustion of methane. CombustionScience and Technology, 182(7):794–821, 2010. → pages 41[116] S. Navarro-Martinez, A. Kronenburg, and F. D. Mare. Conditional MomentClosure for large eddy simulations. Flow, Turbulence and Combustion, 75(1):245–274, 2005. → pages 41157Bibliography[117] K.-J. Nogenmyr, C. Fureby, X. Bai, P. Petersson, R. Collin, and M. Linne.Large eddy simulation and laser diagnostic studies on a low swirl stratifiedpremixed flame. Combustion and Flame, 156(1):25–36, 2009. → pages 38[118] I. V. Novosselov and P. C. Malte. Development and application of aneight-step global mechanism for CFD and CRN simulations oflean-premixed combustors. Journal of Engineering for Gas Turbines andPower, 130(2):021502–021502–9, 2008. → pages 14[119] M. Oberlack, H. Wenzel, and N. Peters. On symmetries and averaging ofthe G-equation for premixed combustion. Combustion Theory andModelling, 5(3):363–383, 2001. → pages 38[120] M. Oevermann, H. Schmidt, and A. Kerstein. Investigation of autoignitionunder thermal stratification using linear eddy modeling. Combustion andFlame, 155(3):370–379, 2008. → pages 71[121] V. Oijen and D. Goey. Modelling of premixed laminar flames usingflamelet-generated manifolds. Combustion Science and Technology, 161(1):113–137, 2000. → pages 14, 44, 65, 74, 100, 114, 129[122] E. S. Oran and J. P. Boris. Numerical Simulation of Reactive Flow.Elsevier Science Publishing Co., Inc., 1987. → pages 5[123] C. C. Paige and M. A. Saunders. LSQR: An algorithm for sparse linearequations and sparse least squares. ACM Trans. Math. Softw., 8(1):43–71,Mar. 1982. → pages 98, 111, 113[124] G. D. 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. → pages 41[125] N. Peters. Numerical and asymptotic analysis of systematically reducedreaction schemes for hydrocarbon flames. In R. Glowinski, B. Larrouturou,and R. Temam, editors, Numerical Simulation of Combustion Phenomena,volume 241 of Lecture Notes in Physics, pages 90–109. Springer BerlinHeidelberg, 1985. ISBN 978-3-540-16073-1. → pages 13, 129[126] N. Peters. Laminar flamelet concepts in turbulent combustion. Symposium(International) on Combustion, 21(1):1231–1250, 1988. → pages 23[127] N. Peters. The turbulent burning velocity for large-scale and small-scaleturbulence. Journal of Fluid Mechanics, 384:107–132, 1999. → pages 23158Bibliography[128] N. Peters. Turbulent Combustion. Cambridge University Press, 2000. →pages 15, 16, 25, 38[129] H. Pitsch. A consistent level set formulation for large eddy simulation ofpremixed turbulent combustion. Combustion and Flame, 143(4):587–598,2005. → pages 38[130] H. Pitsch. Large eddy simulation of turbulent combustion. Annual Reviewof Fluid Mechanics, 38(1):453–482, 2006. → pages 37[131] H. Pitsch and L. D. de Lageneste. Large eddy simulation of premixedturbulent combustion using a level-set approach. Proceedings of theCombustion Institute, 29(2):2001–2008, 2002. → pages 23, 38[132] T. Poinsot and D. Veynanate. Theoretical and Numerical Combustion.Edwards, Philadelphia, 2nd edition, 2005. → pages 5, 6, 26, 37, 38, 41, 98,126, 128, 131[133] T. Poinsot, D. Veynante, and S. Candel. Quenching processes andpremixed turbulent combustion diagrams. Journal of Fluid Mechanics,228:561–606, 7 1991. → pages 23[134] S. Pope. PDF methods for turbulent reactive flows. Progress in Energy andCombustion Science, 11(2):119–192, 1985. → pages 40, 98, 126[135] S. Pope. Computations of turbulent combustion: Progress and challenges.Symposium (International) on Combustion, 23(1):591–612, 1991.Twenty-Third Symposium (International) on Combustion. → pages 41[136] S. B. Pope. Turbulent Flows. Cambridge University Press, 2000. → pages5, 18, 19, 21, 28, 32[137] F. Proch and A. M. Kempf. Numerical analysis of the cambridge stratifiedflame series using artificial thickened flame LES with tabulated premixedflame chemistry. Combustion and Flame, 161(10):2627–2646, 2014. →pages 69[138] 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. → pages 18[139] L. F. Richardson. Weather Prediction by Numerical Process. CambridgeUniversity Press, 1922. → pages 20159Bibliography[140] Y. Saad and H. A. van der Vorst. Iterative solution of linear systems in the20th century. Journal of Computational and Applied Mathematics, 123(12):1–33, 2000. → pages 111, 117[141] P. Sagaut. Large Eddy Simulation for Incompressible Flows.Springer-Verlag Berlin Heidelberg, 2006. ISBN 3540263446. → pages 27,30, 31, 33[142] M. Salehi. Numerical Simulation of Turbulent Premixed Flames withConditional Source-term Estimation. PhD thesis, University of BritishColumbia, 2012. → pages 100, 102, 113, 114, 132[143] M. M. Salehi and W. K. Bushe. Presumed PDF modeling for RANSsimulation of turbulent premixed flames. Combustion Theory andModelling, 14(3):381–403, 2010. → pages 40, 42, 43, 47, 65, 98, 100, 101,103, 104, 111, 127, 130[144] 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, 2012.→ pages 48[145] M. M. Salehi, W. K. Bushe, N. Shahbazian, and C. P. Groth. Modifiedlaminar flamelet presumed probability density function for LES ofpremixed turbulent combustion. Proceedings of the Combustion Institute,34(1):1203–1211, 2013. → pages 40, 42, 43, 65, 98, 100, 101, 103, 104,111, 127, 130[146] A. Sa´nchez, A. Le´pinette, M. Bollig, A. Lin¯a´n, and B. La´zaro. The reducedkinetic description of lean premixed combustion. Combustion and Flame,123(4):436–464, 2000. → pages 14[147] R. Sankaran, E. R. Hawkes, J. H. Chen, T. Lu, and C. K. Law. Structure ofa spatially developing turbulent lean methane-air bunsen flame.Proceedings of the Combustion Institute, 31(1):1291–1298, 2007. → pages63, 93[148] V. Sankaran and S. Menon. Structure of premixed turbulent flames in thethin-reaction-zones regime. Proceedings of the Combustion Institute, 28(1):203–209, 2000. → pages 67, 70, 127, 130[149] U. Schumann. Subgrid-scale model for finite difference simulations ofturbulent flows in plane channels and annuli. Journal of ComputationalPhysics, 18(4):376–404, 1975. → pages 33160Bibliography[150] A. Scotti, C. Meneveau, and D. K. Lilly. Generalized Smagorinsky modelfor anisotropic grids. Physics of Fluids A, 5(9):2306–2308, 1993. → pages31[151] A. Scotti, C. Meneveau, and M. Fatica. Dynamic Smagorinsky model onanisotropic grids. Physics of Fluids, 9(6):1856–1858, 1997. → pages 31[152] L. Selle, G. Lartigue, T. Poinsot, R. Koch, K.-U. Schildmacher, W. Krebs,B. Prade, P. Kaufmann, and D. Veynante. Compressible large eddysimulation of turbulent combustion in complex geometry on unstructuredmeshes. Combustion and Flame, 137(4):489–505, 2004. → pages 39[153] N. Shahbazian, M. M. Salehi, C. P. Groth, O¨. L. Gu¨lder, and W. K. Bushe.Performance of Conditional Source-term Estimation model for LES ofturbulent premixed flames in thin reaction zones regime. Proceedings ofthe Combustion Institute, 35(2):1367–1375, 2015. → pages 40, 42, 43, 98,100, 101, 103, 104, 111, 114, 115, 127, 129, 130, 135[154] M. Sheikhi, T. Drozda, P. Givi, F. Jaberi, and S. Pope. Large eddysimulation of a turbulent non-premixed piloted methane jet flame (SandiaFlame D). Proceedings of the Combustion Institute, 30(1):549 – 556, 2005.→ pages 41[155] J. Smagorinsky. General Circulation Experiments with the PrimitiveEquations. Monthly Weather Review, 91:99, 1963. → pages 33[156] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer,M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, W. C. Gardiner,V. V. Lissianski, and Z. Qin. GRI-Mech 3.0.http://www.me.berkeley.edu/gri mech/, 1997. Accessed: 2016-01-14. →pages 12, 17, 71, 132[157] T. M. Smith and S. Menon. One-dimensional simulations of freelypropagating turbulent premixed flames. Combustion Science andTechnology, 128(1-6):99–130, 1997. → pages 67, 127, 130[158] Y. Sommerer, D. Galley, T. Poinsot, S. Ducruix, F. Lacas, and D. Veynante.Large eddy simulation and experimental study of flashback and blow-off ina lean partially premixed swirled burner. Journal of Turbulence, 5:N37,2004. → pages 39[159] P. Spalart and S. Allmaras. A one-equation turbulence model foraerodynamic flows. In Paper 92-0439, AIAA, 1992. → pages 30161Bibliography[160] D. Spalding. Mixing and chemical reaction in steady confined turbulentflames. Symposium (International) on Combustion, 13(1):649–657, 1971.→ pages 34[161] H. Steiner and W. K. Bushe. Large eddy simulation of a turbulent reactingjet with Conditional Source-term Estimation. Physics of Fluids, 13(3):754–769, 2001. → pages 42, 98, 101, 127[162] N. H. Stern. Stern Review: on the Economics of Climate Change. LondonHM Treasury, 2006. → pages 1[163] M. Sto¨llinger and S. Heinz. PDF modeling and simulation of premixedturbulent combustion. Monte Carlo Methods and Applications, 14(4):343–377, 2008. → pages 41[164] N. Swaminathan and R. Bilger. Assessment of combustion submodels forturbulent non-premixed hydrocarbon flames. Combustion and Flame, 116(4):519–545, 1999. → pages 88[165] N. Swaminathan and R. W. Bilger. Scalar dissipation, diffusion anddilatation in turbulent H2-air premixed flames with complex chemistry.Combustion Theory and Modelling, 5(3):429–446, 2001. → pages 88[166] N. Swaminathan and R. W. Bilger. Analyses of Conditional MomentClosure for turbulent premixed flames. Combustion Theory and Modelling,5(2):241–260, 2001. → pages 41[167] M. Sweeney, S. Hochgreb, and R. Barlow. Cambridge stratified slot burnerdata. 2009. URL http://www.dspace.cam.ac.uk/handle/1810/226470. →pages xiv, 61, 63, 145[168] M. Sweeney, S. Hochgreb, and R. Barlow. The structure of premixed andstratified low turbulence flames. Combustion and Flame, 158(5):935–948,2011. → pages 60, 68, 145[169] M. Sweeney, S. Hochgreb, M. Dunn, and R. Barlow. A comparativeanalysis of flame surface density metrics in premixed and stratified flames.Proceedings of the Combustion Institute, 33(1):1419–1427, 2011. → pages69[170] M. S. Sweeney, S. Hochgreb, M. J. Dunn, and R. S. Barlow. The structureof turbulent stratified and premixed methane/air flames I: Non-swirlingflows. Combustion and Flame, 159(9):2896–2911, 2012. → pages 68, 75,76, 145162Bibliography[171] M. S. Sweeney, S. Hochgreb, M. J. Dunn, and R. S. Barlow. The structureof turbulent stratified and premixed methane/air flames II: Swirling flows.Combustion and Flame, 159(9):2912–2929, 2012. → pages 69[172] M. S. Sweeney, S. Hochgreb, M. J. Dunn, and R. S. Barlow. Multiplyconditioned analyses of stratification in highly swirling methane/air flames.Combustion and Flame, 160(2):322–334, 2013. → pages 75, 76, 145[173] P. Tamadonfar and O¨. L. Gu¨lder. Effects of mixture composition andturbulence intensity on flame front structure and burning velocities ofpremixed turbulent hydrocarbon/air bunsen flames. Combustion andFlame, pages –, 2015. → pages xi, 127, 133, 135[174] B. Thornber, R. Bilger, A. Masri, and E. Hawkes. An algorithm for LES ofpremixed compressible flows using the Conditional Moment Closuremodel. Journal of Computational Physics, 230(20):7687–7705, 2011. →pages 41[175] H. Tsui and W. Bushe. Linear-eddy model formulated probability densityfunction and scalar dissipation rate models for premixed combustion. Flow,Turbulence and Combustion, 93(3):487–503, 2014. → pages 67, 74, 81,127, 130, 131, 132, 137[176] H. Tsui, M. Kamal, S. Hochgreb, and W. Bushe. Direct comparison of PDFand scalar dissipation rates between LEM simulations and experiments forturbulent, premixed methane air flames. Combustion and Flame, 165:208–222, 2016. → pages 137[177] H. P. Tsui and W. K. Bushe. Conditional Source-term Estimation usingdynamic ensemble selection and parallel iterative solution. CombustionTheory and Modelling, 2016. (in press). → pages 137[178] M. Tsurikov. Experimental investigation of the fine-scale structure inturbulent gas-phase jet flows. PhD thesis, University of Texas at Austin,2002. → pages 94[179] D. Vandromme and H. Haminh. Turbulence and Coherent Structures:Selected Papers from “Turbulence 89: Organized Structures andTurbulence in Fluid Mechanics”, Grenoble, 18–21 September 1989, pages507–523. Springer Netherlands, Dordrecht, 1991. ISBN978-94-015-7904-9. → pages 30163Bibliography[180] 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:N4,2004. → pages 39[181] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress inEnergy and Combustion Science, 28(3):193–266, 2002. → pages 5, 34, 35,58, 68, 69, 90, 130, 131[182] W. Vicente, M. Salinas, E. Barrios, and C. Dopazo. PDF modeling of COand NO formation in lean premixed methane flames. Combustion Scienceand Technology, 176(4):585–601, 2004. → pages 41[183] M. Wang, J. Huang, and W. Bushe. Simulation of a turbulent non-premixedflame using Conditional Source-term Estimation with trajectory generatedlow-dimensional manifold. Proceedings of the Combustion Institute, 31(2):1701–1709, 2007. → pages 42, 98, 101, 127[184] H. Weller, G. Tabor, A. Gosman, and C. Fureby. Application of aflame-wrinkling LES combustion model to a turbulent mixing layer.Symposium (International) on Combustion, 27(1):899–907, 1998. → pages37[185] H. G. Weller. OpenFOAM 2.1.1. http://www.openfoam.org/version2.1.1/,2012. Accessed: 2016-01-14. → pages 99, 102, 114, 128[186] F. M. White. Fluid Mechanics. Boston: McGraw-Hill Book Company, 5edition, 2003. → pages 5, 18[187] D. C. Wilcox. Reassessment of the scale-determining equation foradvanced turbulence models. AIAA Journal, 26(11):1299–1310, 1988. →pages 30[188] C. R. Wilke. A viscosity equation for gas mixtures. The Journal ofChemical Physics, 18(4):517–519, 1950. → pages 9[189] S. L. Yilmaz, M. B. Nik, P. Givi, and P. A. Strakey. Scalar filtered densityfunction for large eddy simulation of a bunsen burner. Journal ofPropulsion and Power, 26(1):84–93, Jan 2010. → pages 41[190] F. T. Yuen and O¨. L. Gu¨lder. Premixed turbulent flame front structureinvestigation by Rayleigh scattering in the thin reaction zone regime.Proceedings of the Combustion Institute, 32(2):1747–1754, 2009. → pages99164Bibliography[191] F. T. C. Yuen. Experimental Investigation of the Dynamics and Structure ofLean-Premixed Turbulent Combustion. PhD thesis, University of Toronto,2009. → pages 133[192] K. K. yun Kuo and R. Acharya. Fundamentals of Turbulent andMulti-Phase Combustion. Wiley, 1st edition, 2012. → pages 6[193] J. Zhang, F. Gao, G. Jin, and G. He. Conditionally statistical description ofturbulent scalar mixing at subgrid-scales. Flow, Turbulence andCombustion, 93(1):125–140, 2014. → pages 68[194] R. Zhou, S. Balusamy, M. S. Sweeney, R. S. Barlow, and S. Hochgreb.Flow field measurements of a series of turbulent premixed and stratifiedmethane/air flames. Combustion and Flame, 160(10):2017–2028, 2013. →pages x, 69, 75, 145165
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Turbulent premixed combustion simulation with Conditional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Turbulent premixed combustion simulation with Conditional Source-term Estimation and Linear-Eddy Model… Tsui, Hong P. 2016
pdf
Page Metadata
Item Metadata
Title | Turbulent premixed combustion simulation with Conditional Source-term Estimation and Linear-Eddy Model formulated PDF and SDR models |
Creator |
Tsui, Hong P. |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | Computational fluid dynamics (CFD) is indispensable in the development of complex engines due to its low cost and time requirement compared to experiments. Nevertheless, because of the strong coupling between turbulence and chemistry in premixed flames, the prediction of chemical reaction source terms continues to be a modelling challenge. This work focuses on the improvement of turbulent premixed combustion simulation strategies requiring the use of presumed probability density function (PDF) models. The study begins with the development of a new PDF model that includes the effect of turbulence, achieved by the implementation of the Linear-Eddy Model (LEM). Comparison with experimental burners reveals that the LEM PDF can capture the general PDF shapes for methane-air combustion under atmospheric conditions with greater accuracy than other presumed PDF models. The LEM is additionally used to formulate a new, pseudo-turbulent scalar dissipation rate (SDR) model. Conditional Source-term Estimation (CSE) is implemented in the Large Eddy Simulation (LES) of the Gülder burner as the closure model for the chemistry-turbulence interactions. To accommodate the increasingly parallel computational environments in clusters, the CSE combustion module has been parallelised and optimised. The CSE ensembles can now dynamically adapt to the changing flame distributions by shifting their spatial boundaries and are no longer confined to pre-allocated regions in the simulation domain. Further, the inversion calculation is now computed in parallel using a modified version of an established iterative solver, the Least-Square QR-factorisation (LSQR). The revised version of CSE demonstrates a significant reduction in computational requirement — a reduction of approximately 50% — while producing similar solutions as previous implementations. The LEM formulated PDF and SDR models are subsequently implemented in conjunction with the optimised version of CSE for the LES of a premixed methane-air flame operating in the thin reaction zone. Comparison with experimental measurements of temperature reveals that the LES results are very comparable in terms of the flame height and distribution. This outcome is encouraging as it appears that this work represents a significant step towards the correct direction in developing a complete combustion simulation strategy that can accurately predict flame characteristics in the absence of ad hoc parameters. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-01-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0340685 |
URI | http://hdl.handle.net/2429/60295 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_february_tsui_hong.pdf [ 2.82MB ]
- Metadata
- JSON: 24-1.0340685.json
- JSON-LD: 24-1.0340685-ld.json
- RDF/XML (Pretty): 24-1.0340685-rdf.xml
- RDF/JSON: 24-1.0340685-rdf.json
- Turtle: 24-1.0340685-turtle.txt
- N-Triples: 24-1.0340685-rdf-ntriples.txt
- Original Record: 24-1.0340685-source.json
- Full Text
- 24-1.0340685-fulltext.txt
- Citation
- 24-1.0340685.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0340685/manifest