Improving combustion simulations through a novel principal component analysis-based reduction technique and a new pressure projection algorithm

Turbulent combustion modeling is a complex computational problem. Several factors including the large number of unknowns and equations, stiffness in the chemical source terms, and turbulence-chemistry interaction combine to make simulation of turbulent combustion a grand-challenge problem. Direct Numerical Simulation (DNS) is the most accurate approach to simulate turbulent combustion processes. This approach solves for the full set of chemical variables in the system and is fully resolved in space and time; therefore, it is computationally expensive. There are several models trying to increase the efficiency of turbulent combustion modeling with reducing the number of unknowns, reducing the stiffness of the problem, or decreasing the resolution with the least error possible. In this research, two novel models are introduced to increase the efficiency of turbulent combustion modeling in the Large Eddy Simulation (LES) context. Each method tries to make the modeling more efficient in a different aspect. The first one is a method to reduce the number of species equations that must be solved, via application of Principal Component Analysis (PCA). This technique provides a robust methodology to reduce the number of species equations by identifying correlations in state-space and defining new variables that are linear combinations of the original variables. Here we first present results from \emph{a priori} studies to show the strengths and weaknesses of such a modeling approach. Results suggest that the PCA-based model can identify manifolds that exist in state space which are insensitive to filtering, suggesting that the model is directly applicable for use in Large Eddy Simulation. Second, we explore the invariance of the manifolds identified by PCA with respect to the problem's parameters. In order to simulate a turbulent process using a PCA-based model, the PCA mapping should be trained using an empirical dataset. This \emph{a priori} study clarifies the important factors for choosing a training dataset. Results indicate that, for given reactant compositions and temperatures, over modest ranges of Reynolds number where the combustion regime does not change dramatically, PCA-derived manifolds are invariant with respect to Reynolds number. It also further confirms PCA manifolds invariance to the filter width, which is an interesting result that suggests the applicability of the model in LES. Finally, an \emph{a posteriori} study of PCA is presented as a combustion model applied to a nonpremixed CO/H$_2$ temporally evolving jet flame with extinction and reignition. As a basis for comparison, results from detailed chemistry calculations are compared with the PCA-transport results to verify the model and evaluate its performance. Invariance of the model's error to the Reynolds number, the number of retained PCs, the PCA scaling factor, and the training dataset is evaluated in this research. The second proposed method is a new explicit variable-density pressure projection method with a focus on transient low-Mach-number reacting flows in order to avoid implicitness and iterative schemes. The method is based on solving the pressure Poisson equation and is suitable for implementation in fully explicit codes. It has been verified against novel closed-form analytical solutions as well as manufactured solutions for time-varying, variable-density test cases. These cases range from predominantly diffusive to purely convective conditions, and are suitable for use in verification of transient, variable-density flow codes such as those employed in low-Mach-number turbulent combustion simulations. Finally, the algorithm has been used to simulate an annular nonpremixed, nonreacting, variable-density jet flow to qualitatively demonstrate its performance on a practical case.

IMPROVING COMBUSTION SIMULATIONS THROUGH A NOVEL PRINCIPAL COMPONENT ANALYSIS-BASED REDUCTION TECHNIQUE AND A NEW PRESSURE PROJECTION ALGORITHM by Amir Biglari A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2015 Copyright c Amir Biglari 2015 All Rights Reserved The Univers i ty of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Amir Biglari has been approved by the following supervisory committee members: James Sutherland , Chair 4/28/2015 Date Approved Philip Smith , Member 5/19/2015 Date Approved Jeremy Thornock , Member 4/28/2015 Date Approved John McLennan , Member 4/28/2015 Date Approved 􀀭􂴀􀁄􄐀􀁐􅀀􀁈􄠀􀁖􅘀􀀃􀍒R􀀑􁄠 Stoll , Member 4/28/2015 Date Approved and by Milind Deo , Chair/Dean of the Department/College/School of Chemical Engineering and by David B. Kieda, Dean of The Graduate School. ABSTRACT Turbulent combustion modeling is a complex computational problem. Several factors including the large number of unknowns and equations, sti ness in the chemical source terms, and turbulence-chemistry interaction combine to make simulation of turbulent com-bustion a grand-challenge problem. Direct Numerical Simulation (DNS) is the most ac-curate approach to simulate turbulent combustion processes. This approach solves for the full set of chemical variables in the system and is fully resolved in space and time; therefore, it is computationally expensive. There are several models trying to increase the e ciency of turbulent combustion modeling with reducing the number of unknowns, reducing the sti ness of the problem, or decreasing the resolution with the least error possible. In this research, two novel models are introduced to increase the e ciency of turbulent combustion modeling in the Large Eddy Simulation (LES) context. Each method tries to make the modeling more e cient in a di erent aspect. The first one is a method to reduce the number of species equations that must be solved, via application of Principal Component Analysis (PCA). This technique provides a robust methodology to reduce the number of species equations by identifying correlations in state-space and defining new variables that are linear combinations of the original variables. Here we first present results from a priori studies to show the strengths and weaknesses of such a modeling approach. Results suggest that the PCA-based model can identify manifolds that exist in state space which are insensitive to filtering, suggesting that the model is directly applicable for use in Large Eddy Simulation. Second, we explore the invariance of the manifolds identified by PCA with respect to the problem's parameters. In order to simulate a turbulent process using a PCA-based model, the PCA mapping should be trained using an empirical dataset. This a priori study clarifies the important factors for choosing a training dataset. Results indicate that, for given reactant compositions and temperatures, over modest ranges of Reynolds number where the combustion regime does not change dramatically, PCA-derived manifolds are invariant with respect to Reynolds number. It also further confirms PCA manifolds invariance to the filter width, which is an interesting result that suggests the applicability of the model in LES. Finally, an a posteriori study of PCA is presented as a combustion model applied to a nonpremixed CO/H2 temporally evolving jet flame with extinction and reignition. As a basis for comparison, results from detailed chemistry calculations are compared with the PCA-transport results to verify the model and evaluate its performance. Invariance of the model's error to the Reynolds number, the number of retained PCs, the PCA scaling factor, and the training dataset is evaluated in this research. The second proposed method is a new explicit variable-density pressure projection method with a focus on transient low-Mach-number reacting flows in order to avoid implic-itness and iterative schemes. The method is based on solving the pressure Poisson equation and is suitable for implementation in fully explicit codes. It has been verified against novel closed-form analytical solutions as well as manufactured solutions for time-varying, variable-density test cases. These cases range from predominantly di usive to purely convective conditions, and are suitable for use in verification of transient, variable-density flow codes such as those employed in low-Mach-number turbulent combustion simulations. Finally, the algorithm has been used to simulate an annular nonpremixed, nonreacting, variable-density jet flow to qualitatively demonstrate its performance on a practical case. iv To my dearest parents and lovely wife. CONTENTS ABSTRACT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : iii LIST OF FIGURES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : ix LIST OF TABLES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xii ACKNOWLEDGMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xiv CHAPTERS 1. INTRODUCTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 1.1 Principal Component Analysis-Based Model for Turbulent Combustion Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 An Explicit Pressure Projection Model for Low-Mach Variable Density Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2. A FILTER-INDEPENDENT MODEL IDENTIFICATION TECHNIQUE FOR TURBULENT COMBUSTION MODELING : : : : : : : : : : : : : : : : : : : : 9 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Introduction and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Parameterization Using Principal Component Analysis . . . . . . . . . . . . . . . 12 2.4.1 Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.2 Multivariate adaptive nonlinear regression . . . . . . . . . . . . . . . . . . . . 17 2.5 Filtering and Turbulent Closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5.1 PCA sensitivity to filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5.2 Turbulent closure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5.3 Comparison with SLFM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.6 Considerations for Model Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.9 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3. IDENTIFICATION OF INVARIANT MANIFOLDS IN TURBULENT COMBUSTION SYSTEMS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 38 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3 Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.1 Dependence on chemical composition . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4.2 Invariance to Reynolds number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.4.3 Invariance to Reynolds number and filter width . . . . . . . . . . . . . . . . . 45 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4. AN A-POSTERIORI EVALUATION OF PRINCIPAL COMPONENT ANALYSIS-BASED MODELS FOR TURBULENT COMBUSTION SIMULATIONS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 51 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3.1 Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3.2 Formulation of PC transport equations . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.3 One-dimensional turbulence simulation . . . . . . . . . . . . . . . . . . . . . . . 57 4.4 Computational Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.5.1 E ect of the Reynolds number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.5.2 E ect of the number of retained PCs . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.5.3 E ect of scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.5.4 E ect of training data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.7 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.8 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5. AN EXPLICIT PRESSURE PROJECTION METHOD FOR LOW-MACH VARIABLE-DENSITY FLOWS : : : : : : : : : : : : : : : : : : : : : : : : : : : : 73 5.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Estimation of the Pressure Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5.1 Divergence formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.5.2 Density formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.5.3 Unified formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.6 Benchmark Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.6.1 Convection of a nonreacting mixture . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.6.2 One-dimensional manufactured solution . . . . . . . . . . . . . . . . . . . . . . 82 5.6.3 Two-dimensional manufactured solution . . . . . . . . . . . . . . . . . . . . . . 83 5.7 Determining 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 vii 5.8.1 Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.8.2 One-dimensional manufactured solution . . . . . . . . . . . . . . . . . . . . . . 88 5.8.3 Two-dimensional manufactured solution . . . . . . . . . . . . . . . . . . . . . . 89 5.9 An Annular Jet-flow Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.11 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.12 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6. CONCLUSION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 96 6.1 PCA-Based Reduction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.2 Explicit Pressure Projection Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3 Recommendations for Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3.1 PCA-based reduction model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3.2 Explicit pressure projection method . . . . . . . . . . . . . . . . . . . . . . . . . 99 APPENDICES A. PC SOURCE TERMS PARAMETERIZATION IMPROVEMENTS : : : : : : 100 B. SUMMARY OF THE PRESSURE PROJECTION ALGORITHM : : : : : : : 103 C. SOURCE TERMS OF THE ONE-DIMENSIONAL MANUFACTURED SOLUTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 104 D. SOURCE TERMS OF THE TWO-DIMENSIONAL MANUFACTURED SOLUTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 106 viii LIST OF FIGURES 1.1 Conventional flow regimes and their associated Mach numbers. . . . . . . . . . . 4 2.1 Illustration of the principal components of a hypothetical 2D dataset where we retain a single principal component. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Eigenvalue magnitude (left axis, bars) and percent variance captured (right axis, line) by retaining the given number of components. . . . . . . . . . . . . . . . 13 2.3 Comparison of PCA and MARS reconstructions for OH mass fraction for a two-dimensional model based on principal components 1 and 2. VAST scaling was used. (a) PCA (linear) reconstruction of YOH in ( 1; 2)-space. (b) MARS (nonlinear) reconstruction of YOH in ( 1; 2)-space. . . . . . . . . . . . 19 2.4 E ect of filtering on temperature profile for a specific time and realization from the temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Instantaneous profile of K􀀀¯K ¯K = K0 ¯K indicating the magnitude of the unresolved kinetic energy at = x = 4 and 16. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6 Changes in largest (most important) components of the first three eigenvec-tors with respect to changes in normalized filter width in temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 R2 value changes with respect to the changes in normalized filter width (nor-malized with grid spacing length) for several variables in temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.8 Original (solid) and reconstructed (dashed) profiles for CO2 and OH for no filtering and a filter width of = x = 16 and n = 2. . . . . . . . . . . . . . . . . . . . 26 2.9 S¯ 1 (¯ 1; ¯ 2) obtained directly from the data (points) as well as the prediction based on PCA/MARS (surface) for various filter widths using the temporal CO/H2 dataset. Figure 2.10 shows the corresponding R2 values. . . . . . . . . . . 27 2.10 R2 value changes with respect to the changes in normalized filter width, = x for the source term of the first and the second PCs, in temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.11 Original (solid) and reconstructed (dashed) profiles for S 1 and S 2 for no filtering and a filter width of = x = 16 and n = 2. . . . . . . . . . . . . . . . . . . . 28 2.12 Weights for the first three eigenvectors (which define the first three PCs) for Pareto scaling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.13 Parity plots for temperature (left) and OH (right) reconstructions using PCA/- MARS (top row) and SLFM (bottom row). . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.14 Reconstruction of temperature (left) and OH mass fraction (right) for a 2D PCA/MARS model (top) and the SLFM model (bottom) for the temporal ODT dataset [14]. The R2 value of these reconstructions are reported on each plot as well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.15 R2 value changes with respect to the changes in normalized filter width ( = x) for several state variables comparing PCA and MARS results with SLFM and SLFM- -PDF results for the temporal CO/H2 dataset. See Table 2.6 for more information. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Structure of the first four PCs for the JHC data with 3% and 9% O2 in the co-flow, showing the variability of the PC structure. . . . . . . . . . . . . . . . . . . . 42 3.2 Largest five contributions to the first three PCs in the PCA for the three Re cases at various filter widths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.1 Contour plots of temperature, YH2O, YO2 , and YOH for di erent Re numbers (case L (2510), M (4478) and H (9079)). Each contour plot shows the full-chemistry ODT solution on its left half-plane and PC transported ODT solution on the right half-plane. Here, N = 2 and the model is trained on case M. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2 Mean values of temperature, YH2O, YO2 , and YOH conditioned on the mixture fraction (Z) for di erent Re numbers. Each plot compares full-chemistry solution and PC transported solution at several points in time. Here is the nondimensional time scale of the problem, t=t j. See also Figure 4.1. . . . . . . . 61 4.3 Contour plots for case M using N = 1 (top row) and N = 3 (bottom row). These are comparable to the middle column (case M) of Figure 4.1. The left half-plane in each contour plot shows the full-chemistry while the PC transported solution is on the right side. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.4 Mean temperature, YH2O, YO2 , and YOH conditioned on the mixture fraction (Z) at three di erent key points in time (at the beginning, = 11, and the end, = 21, of extinction and after reignition is completed, = 46, where is the nondimensional time scale of the problem, t=t j). . . . . . . . . . . . . . . . . . . 63 4.5 Evolution of the PC manifolds from 1 retained PC to 3 retained PCs. . . . . . . 64 4.6 Mean temperature, YH2O, YO2 , and YOH conditioned on the mixture fraction. Each plot shows the statistical means of these variable solutions using VAST and Pareto scaling in the PC transported simulation and the full-chemistry ODT simulation at a specific time (at the beginning, = 11, and the end, = 21, of extinction and after reignition is completed, = 46). . . . . . . . . . . 65 4.7 Profiles of T, YH2O, and YOH at the onset of extinction ( = 11) and reignition ( = 21). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.8 The manifold coverage of the two datasets considered in this work. Note that the full-chemistry ODT dataset size is reduced 25 times and SLFM dataset size is reduced 5 times for better visibility. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 x 4.9 Conditional means with respect to the mixture fraction. Each plot shows the statistical means of the variable solutions using the SLFM dataset and the full-chemistry ODT dataset section 4.4 from the PC transported simulation and compares them against the full-chemistry ODT simulation at three dif-ferent times (at the beginning, = 11, and the end, = 21, of the extinction and after reignition is completed, = 46). . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.1 Time-averaged L2 norm of the x-velocity error versus 0. The plot is shown for the one-dimensional MMS described in 5.6.2 . . . . . . . . . . . . . . . . . . . . . 85 5.2 Comparison between the numerical and analytical solutions for the convec-tion problem defined in section 5.6.1. Results are shown at ¯t = 0; 1:6; 3:2; and 5. The first two columns from the left correspond to the mixture fraction and its absolute error. The last two columns correspond to the velocity and its absolute error, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3 Temporal and spatial order of accuracy for the convection benchmark re-ported in section 5.6.1. Results were obtained using the mixture fraction. . . . 87 5.4 Temporal (left) and spatial (right) order of accuracy for the 1D MMS bench-mark reported in section 5.6.2 for the scalar field (i.e., mixture fraction). . . . . 88 5.5 Comparison between the numerical and analytical solutions for the 1D MMS problem defined in section 5.6.2. Results are shown at ¯t = 0; 4; 8; and 12:5. The first two columns from the left correspond to the mixture fraction and its absolute error. The last two columns correspond to the velocity and its absolute error, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.6 Temporal and spatial order of accuracy for the 2D MMS benchmark reported in section 5.6.3. Results were produced using the mixture fraction. . . . . . . . . 90 5.7 Evolution of the dimensionless density (first row), mixture fraction (second row), and axial velocity (third row) for the two-dimensional MMS (section 5.6.3). Results are reported at the dimensionless times ¯t = 0; 0:125; 0:25; 0:375; and 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.8 Temporal and spatial evolution of the absolute error for the density (first row), mixture fraction (second row), and axial velocity (third row). All reported quantities are in dimensionless form. Columns correspond to the dimensionless times ¯t = 0; 0:125; 0:25; 0:375; and 1. . . . . . . . . . . . . . . . . . 92 5.9 Volume rendered and filled contour plots (on the jet mid-plane) of the mix-ture fraction for the annular jet-flow simulation. Results are shown at t = 0.5 s (top), 2, and 4 s (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 xi LIST OF TABLES 2.1 Brief descriptions of the scaling options considered here. . . . . . . . . . . . . . . . 14 2.2 R2 values for PCA projection of state variables with di erent scaling meth-ods using 2 PCs on the temporal CO/H2 dataset. . . . . . . . . . . . . . . . . . . . . . . 15 2.3 R2 values for MARS regression of state variables on principal components with di erent scaling methods in PCA using 2 PCs on the temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4 R2 values for MARS regression of source terms on principal components with di erent scaling methods in PCA using 3 PCs on the temporal CO/H2 dataset [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 R2 values for di erent variables in flame D showing that no closure is re-quired to reconstruct the original variables, from the principal components, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Summary of models compared in Figure 2.15. . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 R2 values for the JHC 3% and 9% O2 datasets. R2 3% and R2 9% indicate recon-struction using the PCA obtained from the 3% and 9% O2 cases, respectively. 43 3.2 Reconstruction accuracy, R2, for the flames C, D, E and F datasets by the PCA reduction. Note that R2C and R2 F refer to the accuracy by which vari-ables are reconstructed (via 2.1) using the PCA obtained for flame C and F, respectively . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 R2 values across a range of filter widths and Reynolds numbers for the CO/H2 datasets. R2d irect indicates the R2 values for reconstructions from PCA obtained directly on the data at the given eRe and . R2 H;1 indicates reconstruction using a PCA obtained on case H (high Re) with = 1. R2 L;32 indicates a reconstruction using a PCA obtained on case L (low Re) with = 32. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1 R2 values for MARS regression of state-variables for N = 1 to N = 3 with the Pareto scaling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 R2 values for MARS regression of state-variables with VAST and Pareto scaling for N = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.3 R2 values for MARS regression of PC source terms with VAST and Pareto scaling for N = 1 and N = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Parameter values for the two-dimensional MMS problem. . . . . . . . . . . . . . . 89 A.1 R2 values for MARS regression of PC source terms with VAST scaling for N = 1 to N = 3 comparing the classic PCA against hybrid PCA. . . . . . . . . 102 A.2 R2 values for MARS regression of PC source terms with VAST scaling for N = 2 and N = 3 using the coupled PC source term method. . . . . . . . . . . . 102 xiii ACKNOWLEDGMENTS First of all, I want to thank my parents and wife for their unconditional love, support, encouragement, and attention throughout my life and PhD program. I would like to express my special appreciation and thanks to my adviser, Dr. James Sutherland, who has been a tremendous mentor for me. I would like to thank him for his invaluable guidance, all the support and patience, his mentorship, and for sharing his knowledge with me to enlighten my path throughout my PhD research. I take this opportunity to express gratitude to all of my committee members, Dr. Philip Smith, Dr. Jeremy Thornock, Dr. Rob Stoll, and Dr. John McLennan, for their support and helpful comments. I would also like to thank Dr. Tony Saad, who was involved in my variable density pressure projection research. I appreciate his helpful advice and unconditional assistance in a co-operative atmosphere in this research. This work has been supported by financial funds from the Department of Energy un-der award number DE-NT0005015, the National Science Foundation PetaApps project 0904631, and the National Nuclear Security Administration through DOE Research Grant DE-NA0000740 and award DE-NA0002375. CHAPTER 1 INTRODUCTION Turbulent combustion problems have several complexities such as large number of unknowns and equations, chemical sti ness, coupling of di erent properties, and turbulent mixing that make their modeling and simulation challenging. There are also numerical limitations in turbulent combustion simulations that a turbulent combustion model should deal with. Small time steps, fine grid resolutions, and iterative schemes are examples of these limitations. Several methods have been developed to reduce these complexities and limitations. For instance, mechanism reduction techniques [1-3] reduce the number of chemical reactions and therefore the chemical sti ness in the simulation, while Reynolds-averaged NavierStokes (RANS) [4] and Large Eddy Simulation (LES) [5-7] methods target the temporal and spatial limitations in turbulent simulations and reduce them to increase the e ciency. One of the proposed methods in this research targets the issue of large number of unknowns and equations and tries to reduce them in a turbulent combustion system. The other proposed method's objective is to remove a limitation on pressure projection schemes in turbulent combustion processes by introducing a novel explicit scheme in order to avoid implicitness and iterative methods. 1.1 Principal Component Analysis-Based Model for Turbulent Combustion Modeling Modeling turbulent combustion processes require solution of a large number of equa-tions due to the large number of reacting species present, as well as the computational cost of resolving turbulent flows scales as Re3. Reducing the range of scales that must be resolved as well as the number of equations to be solved is, therefore, of utmost importance to achieve simulations of practical combustion systems. Classical turbulence theory indicates that resolution requirements scale with the Reynolds 2 number as Re3 for isotropic, homogeneous turbulent flow [8]. Species with large Schmidt numbers further increase the range of scales. In addition to the separation of length scales due to turbulence, the large number of species involved in combustion and the sti chem-istry associated with the reactions further increase the cost of direct simulation so that it is prohibitively expensive for all but the simplest of systems. Typically, time averaging (RANS) or spatial filtering (LES) is used to reduce the res-olution requirements. To reduce the number of thermochemical degrees of freedom, there are two broad approaches: Mechanism reduction, where the chemical mechanism is modified to reduce the number of species and the sti ness, and State-space parameterization, where the state of the system is assumed to be param-eterized by a small number of variables (smaller manifold), which are evolved and transported in the simulation calculations. The techniques proposed here fall into the second category, where they seek to obtain a set of variables that parameterizes the thermochemical state, and these variables are then evolved in the CFD calculation. Turbulent combustion is characterized by a spectrum of length and time scales for both the chemistry and fluid dynamics. At large Dahmkoler number (Da), when these scales are segregated and mixing timescales are much slower than chemical timescales, chemical equilibrium prevails. One of the major challenges in turbulent combustion lies in modeling the situation where the mixing and chemical timescales overlap. As the degree of overlap increases (Da ! 1), more timescales become coupled and finite-rate chemistry e ects become increasingly important. Most combustion modeling approaches begin at the large Da limit and then attempt to incorporate some coupling between mixing and reaction. The steady laminar flamelet model is such an example that introduces the scalar dissipation rate as a mixing timescale that perturbs the state of the system away from chemical equilibrium [9]. Many variations on this have been proposed to model regimes where there is increased overlap in timescales. These models frequently add a chemical progress variable to account for this additional coupling [10-14]. 3 The definition of progress variables is a major challenge. One would like each addi-tional progress variable (parameterizing variable) to be "orthogonal" to the previous ones so that it captures information not already represented. Additionally, progress variables should be chosen so that they represent as much of the variation in the system (departure from equilibrium) as possible. Principal Component Analysis (PCA) has been shown as a viable technique to identify progress variables for use in identifying manifolds in combustion [15-17]. In this research, we review PCA as a technique to obtain a reduced parameter set, discuss how PCA can be formulated as a predictive model, and introduce adaptive regression to enable parameteriza-tion of nonlinear functions of the principal components. Then we examine the model in the context of turbulent closure and show that the model is closed, i.e., it requires no explicit closure model for the thermochemistry. Afterwards, we investigate the structure of the PCA (i.e., the definition of progress variables) to see if it is invariant with respect to system parameters such as Da, Re, and filter width. Then the proposed model will be applied to simulate an ODT problem for validation. The model is trained on simulations with detailed chemistry involving 11 species. Multivariate Adaptive Regression Spline (MARS) is used generally for regressing source terms and regenerating state-variables from PCs to capture the nonlinearity e ects. Results of a posteriori studies on a turbulent nonpremixed simulation involving extinction and reignition are presented, together with a study on the e ect of parameters such as the Reynolds number, the scaling method used in training PCA, the number of retained PCs, and the training data itself. 1.2 An Explicit Pressure Projection Model for Low-Mach Variable Density Flows The pressure projection method is one of the most versatile and widely used techniques to resolve the pressure-velocity coupling in the numerical solution of the Navier-Stokes equations. Originally developed by Chorin [18] for incompressible constant density flows, it was later extended to low-Mach-number variable-density flows [19-22]. The importance of this extension is that one is no longer bound to use the notoriously expensive compress-ible algorithms to simulate this category of flow fields. This cost stems from the scale separation between the convective and acoustic speeds thus requiring the use a small time step size to resolve the fast acoustic scales. 4 From a qualitative perspective, the pressure projection method filters out compressible sti ness and all acoustic waves via instantaneous pressure equilibration causing the speed of sound to become infinite. In other words, the scale separation between the convective and acoustic motions disappears and the convective (and di usive) scales become dominant. This approximation, however, applies to a specific class of flow problems with a low Mach number (M < 0:3) (Figure 1.1). Equivalently, this class corresponds to configurations where pressure variations are small or loosely dependent on the density. In this case, one can e ectively exclude the pressure from the equation of state (EOS) and set , (p). This model is known as the low-Mach-number approximation and is the predominant tool used in the majority of combustion applications and industrial reacting processes. Perhaps one of the earliest attempts at solving the low-Mach-number equations was ac-complished by Majda and Sethian [19] who derived the governing equations for low-Mach combustion for both inviscid and viscous conditions. Their algorithm consisted of solving a highly specialized Poisson's equation in conjunction with a nonlinear ordinary di erential equation for the mean pressure. Significant progress was made by Bell and Marcus circa 1992 who developed a second order projection method for variable density flows [20] based on the method developed by Bell et al. [23] for constant density incompressible flows. Although their model did not allow for any dilatation e ects (r u = 0), it was a significant improvement over the predominance of the Boussinesq approximation. Najm [21] devised a conservative predictor-corrector projection scheme for low-Mach reacting flows. Najm et al. [22] and Knio et al. [24] later developed a semi-implicit method for solving reacting flow problems with chemical reactions. Their physical model was based on a zero-Mach-number formulation of the compressible conservation equa- 􀀁􀄀􀀂􀈀􀀃􀌀􀀄􀐀􀀅􀔀􀀆􀘀􀀇􀜀􀀈􀠀􀀉􀤀􀀊􀨀􀀋􀬀 􀀌􀰀 􀀌􀰀􀀎􀸀􀀏􀼀 􀀍􀴀 􀀍􀴀􀀎􀸀􀀐􁀀 􀀑􁄀􀀒􁈀􀀈􀠀􀀓􁌀􀀋􀬀􀀊􀨀􀀔􁐀􀀔􁐀􀀕􁔀􀀉􀤀􀀖􁘀􀀊􀨀 􀀗􁜀􀀇􀜀􀀓􁌀􀀊􀨀􀀋􀬀􀀔􁐀􀀒􁈀􀀆􀘀􀀕􁔀􀀃􀌀 􀀘􁠀􀀙􁤀􀀓􁌀􀀊􀨀􀀋􀬀􀀔􁐀􀀒􁈀􀀆􀘀􀀕􁔀􀀃􀌀 􀀘􁠀􀀕􁔀􀀚􁨀􀀄􀐀􀀛􁬀􀀄􀐀􀀙􁤀􀀓􁌀􀀊􀨀􀀋􀬀􀀔􁐀􀀒􁈀􀀆􀘀􀀕􁔀􀀃􀌀 􀀌􀰀􀀎􀸀􀀜􁰀 􀀖􁘀􀀒􁈀􀀝􁴀􀀅􀔀􀀁􀄀􀀂􀈀􀀃􀌀􀀄􀐀 􀀞􁸀􀀋􀬀􀀂􀈀􀀆􀘀􀀔􁐀􀀒􁈀􀀆􀘀􀀕􁔀􀀃􀌀 􀀟􁼀 􀀍􀴀􀀌􀰀 􀀠􂀀􀀆􀘀􀀃􀌀􀀒􁈀􀀈􀠀􀀓􁌀􀀋􀬀􀀊􀨀􀀔􁐀􀀔􁐀􀀕􁔀􀀉􀤀􀀖􁘀􀀊􀨀 􀀡􂄀􀀒􁈀􀀇􀜀􀀔􁐀􀀔􁐀􀀕􁔀􀀆􀘀􀀊􀨀􀀔􁐀􀀢􂈀 􀀣􂌀􀀆􀘀􀀊􀨀􀀖􁘀􀀂􀈀􀀔􁐀􀀤􂐀􀀕􁔀􀀃􀌀 􀀗􁜀􀀤􂐀􀀊􀨀􀀖􁘀􀀖􁘀􀀂􀈀􀀋􀬀􀀅􀔀􀀘􁠀􀀙􁤀􀀥􂔀􀀋􀬀􀀒􁈀􀀥􂔀􀀙􁤀􀀆􀘀􀀂􀈀􀀈􀠀􀀕􁔀􀀃􀌀􀀔􁐀 􀀑􁄀􀀒􁈀􀀈􀠀􀀉􀤀􀀇􀜀􀀔􁐀􀀤􂐀􀀕􁔀􀀒􁈀􀀆􀘀 􀀫􂬀􀀬􂰀􀀪􂨀􀀯􂼀􀀭􂴀􀀑􁄀􀀞􁸀􀀠􂀀􀀪􂨀􀀰􃀀 􀀦􂘀 􀀧􂜀 􀀨􂠀 􀀨􂠀 􀀩􂤀 􀀅􀔀 􀀑􁄀 􀀪􂨀 􀀁􀄀 􀀫􂬀 􀀬􂰀 􀀭􂴀 􀀗􁜀 􀀗􁜀 􀀠􂀀 􀀡􂄀 􀀨􂠀 􀀭􂴀 􀀅􀔀 􀀣􂌀 􀀨􂠀 􀀮􂸀 􀀪􂨀 􀀬􂰀 􀀠􂀀 􀀞􁸀 􀀘􁠀 􀀁􀄀 􀀗􁝆 Figure 1.1: Conventional flow regimes and their associated Mach numbers. 5 tions. During the same period, Almgren et al. [25] crafted an adaptive projection method for variable density, incompressible flows in which they used a modified fractional step scheme, called subcycling in time. In addition, they used cell-centered variables for all quantities except for the pressure that was located at the cell corners and staggered in time. Advection-di usion equations were called upon to predict intermediate velocities. Those were subsequently projected onto a space of approximately divergence-free vector fields. Pierce and Moin [11] proposed another semi-implicit scheme to solve variable-density flows. Their technique consisted of solving the pressure Poisson equation for an additive correction of pressure and momentum fields by time-splitting the momentum equation. Furthermore, they used a staggered grid for the velocity in both space and time. Their approach is based on a mixture of ideal gases in thermodynamic equilibrium and chemical nonequilibrium. Later, Shunn et al. [26] presented a spatially-collocated, unstructured version of the reacting flow algorithm of Pierce and Moin. Their variables are advanced in time using a semi-implicit fractional-step method similar to [11]. Pressure and density are decoupled by defining the density through an EOS expressed in terms of transported scalars, which can be defined by an analytical expression, or it may be precomputed and tabulated. Then, a nonlinear solver with Picard iterations is applied at each time step to converge the nonlinear system of equations. This literature survey reveals that the existing approaches for solving low-Mach vari-able density flow make exclusive use of fully- or semi-implicit schemes to advance the transported variables which use iterative algorithms to solve the discretized system of nonlinear equations. Iterative algorithms do not completely satisfy the discrete equations and contain inevitable residual errors. This is because in these approaches, the solution is iteratively refined from an initial guess until it approximately satisfies all the equations at given a point in time. It can be significantly expensive to converge all the solvers to machine precision; therefore, they only conduct a limited number of iterations that results in residual errors. Furthermore, the surveyed methods use fractional time steps along with time staggered calculations leading to increased coding and formulation complexity, especially in the context of large-scale codes. Another limitation in most of the surveyed methods is that they are based on a specific form of the EOS, which needs further e ort to generalize them for any EOS. There are few approaches [27, 28] where they di erentiate 6 a general form of EOS to calculate a term for velocity divergence field or [26] where they use a tabulated EOS in order to generalize their formulation for any type of EOS. This research proposes a general explicit pressure projection method for low-Mach number variable density flows. The method is intended to accommodate an arbitrary equation of state and provide an e cient procedure of calculating the pressure field. The way that the EOS is utilized here is a modified version of [26] where we fully converge on the EOS to ensure consistency between the density and scalars of the system. The method is then verified by two valuable verification tools that have been first provided in literature by the author. 1.3 References [1] K. Kuo, Principles of Combustion. Wiley, 1986. [2] W. Jones and S. Rigopoulos, "Rate-controlled constrained equilibrium: formulation and application to nonpremixed laminar flames," Combustion and Flame, vol. 142, no. 3, pp. 223-234, 2005. [3] S. Rigopoulos, "The rate-controlled constrained equilibrium RCCE method for reduc-ing chemical kinetics in systems with time-scale separation," International Journal for Multiscale Computational Engineering, vol. 5, no. 1, pp. 11-18, 2007. [4] O. Reynolds, "On the dynamical theory of incompressible viscous fluids and the determination of the criterion," Philosophical Transactions of the Royal Society of London. A, vol. 186, pp. 123-164, 1895. [5] J. Smagorinsky, "General circulation experiments with the primitive equations: I.the basic experiment," Monthly weather review, vol. 91, no. 3, pp. 99-164, 1963. [6] J. Deardor , "A numerical study of three-dimensional turbulent channel flow at large reynolds numbers," Journal of Fluid Mechanics, vol. 41, no. 2, pp. 453-480, 1970. [7] H. Pitsch, "Large-eddy simulation of turbulent combustion," Annual Review of Fluid Mechanics, vol. 38, pp. 453-482, 2006. [8] H. Tennekes and J. Lumley, A First Course in Turbulence. MIT Press, 1972. [9] N. Peters, "Laminar di usion flamelet models in non-premixed turbulent combus-tion," Progress in Energy and Combustion Science, vol. 10, pp. 319-339, 1984. [10] J. van Oijen and L. de Goey, "Modelling of premixed laminar flames using flamelet-generated manifolds," Combustion Science and Technology, vol. 161, pp. 113-137, 2000. 7 [11] C. Pierce, "Progress-variable approach for large-eddy simulation of turbulent com-bustion," Ph.D. dissertation, Stanford University, 2001. [12] B. Fiorina, R. Baron, O. Gicquel, D. Thevenin, S. Carpentier, and N. Darabiha, "Mod-elling non-adiabatic partially premixed flames using flame-prolongation of ILDM," Combustion Theory and Modelling, vol. 7, pp. 449-470, 2003. [13] J. Sutherland, P. Smith, and J. Chen, "A quantitative method for a priori evaluation of combustion reaction models," Combustion Theory and Modelling, vol. 11, no. 2, pp. 287-303, 2007. [14] B. Cuenot, F. Egolfopoulos, and T. Poinsot, "An unsteady flamelet model for non-premixed combustion," Combustion Theory and Modelling, vol. 4, pp. 77-97, 2000. [15] A. Parente, J. Sutherland, L. Tognotti, and P. Smith, "Identification of low-dimensional manifolds in turbulent flames," Proceedings of the Combustion Institute, vol. 32, pp. 1579-1586, 2009. [16] J. Sutherland and A. Parente, "Combustion modeling using principal component analysis," Proceedings of the Combustion Institute, vol. 32, pp. 1563-1570, 2009. [17] A. Parente, J. Sutherland, B. Dally, L. Tognotti, and P. Smith, "Investigation of the MILD combustion regime via principal component analysis," Proceedings of the Combustion Institute, vol. 33, no. 2, pp. 3333-3341, 2011. [Online]. Available: http://dx.doi.org/10.1016/j.proci.2010.05.108 [18] A. Chorin, "Numerical solution of the Navier-Stokes equations," Mathematics of Computation, vol. 22, no. 104, pp. 745-762, 1968. [19] A. Majda and J. Sethian, "The derivation and numerical solution of the equations for zero mach number combustion," Combustion Science and Technology, vol. 42, no. 3-4, pp. 185-205, 1985. [20] J. Bell and D. Marcus, "A second-order projection method for variable-density flows," Journal of Computational Physics, vol. 101, no. 2, pp. 334-348, 1992. [21] H.N.Najm, "A conservative low mach number projection method for reacting flow modelling," Sandia National Labs., Livermore, CA (United States), Tech. Rep., 1995. [22] H. Najm, P. Wycko , and O. Knio, "A semi-implicit numerical scheme for reacting flow: I. sti chemistry," Journal of Computational Physics, vol. 143, pp. 381-402, 1998. [23] J. Bell, P. Colella, and H. Glaz, "A second-order projection method for the incom-pressible Navier-Stokes equations," Journal of Computational Physics, vol. 85, no. 2, pp. 257-283, 1989. [24] O. Knio, H. Najm, and P. Wycko , "A semi-implicit numerical scheme for reacting flow: Ii. sti , operator-split formulation," Journal of Computational Physics, vol. 154, 8 no. 2, pp. 428-467, 1999. [25] A. Almgren, J. B. Bell, P. Colella, L. Howell, and M. Welcome, "A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations," Journal of Computational Physics, vol. 142, pp. 1-46, 1998. [26] L. Shunn, F. Ham, and P. Moin, "Verification of variable-density flow solvers using manufactured solutions," Journal of Computational Physics, vol. 231, pp. 3801-3827, 2012. [27] P. McMurtry, W. Jou, J. Riley, and R. Metcalfe, "Direct numerical simulations of a reacting mixing layer with chemical heat release," AIAA Journal, vol. 24, pp. 962- 970, 1985. [28] J. Bell, M. Day, C. Rendleman, S. Woosley, and M. A. Zingale, "Adaptive low mach number simulations of nuclear flame microphysics," Journal of Computational Physics, vol. 195, pp. 677-694, 2004. CHAPTER 2 A FILTER-INDEPENDENT MODEL IDENTIFICATION TECHNIQUE FOR TURBULENT COMBUSTION MODELING2 2.1 Abstract In this paper, we address a method to reduce the number of species equations that must be solved via application of Principal Component Analysis (PCA). This technique provides a robust methodology to reduce the number of species equations by identifying correlations in state-space and defining new variables that are linear combinations of the original variables. We show that applying this technique in the context of Large Eddy Simulation allows for a mapping between the reduced variables and the full set of variables that is insensitive to the size of filter used. This is notable since it provides a model to map state variables to progress variables that is a closed model. As a linear transformation, PCA allows us to derive transport equations for the principal components, which have source terms. These source terms must be parameterized by the reduced set of principal components themselves. We present results from a priori studies to show the strengths and weaknesses of such a modeling approach. Results suggest that the PCA-based model can identify manifolds that exist in state space which are insensitive to filtering, suggesting that the model is directly applicable for use in Large Eddy Simulation. However, the resulting source terms are not parameterized with an accuracy as high as the state variables. 2The material presented in this chapter has been accepted and available in Combustion and Flame 159 (2012) 19601970. Reprinted and adapted with permission from Elsevier. 10 2.2 Introduction and Background Modeling turbulent combustion processes requires solution of a large number of equa-tions due to the large number of reacting species present. Furthermore, the computational cost of resolving turbulent flows scales as Re3. Reducing the range of scales that must be resolved as well as the number of equations to be solved is, therefore, of utmost importance to achieve simulations of practical combustion systems. Classical turbulence theory indicates that resolution requirements scale with the Reynolds number as Re3 for isotropic, homogeneous turbulent flow [1]. Species with large Schmidt numbers further increase the range of scales. In addition to the separation of length scales due to turbulence, the large number of species involved in combustion and the sti chem-istry associated with the reactions further increase the cost of direct simulation so that it is prohibitively expensive for all but the simplest of systems. Typically, time averaging (RANS) or spatial filtering (LES) is used to reduce the res-olution requirements. To reduce the number of thermochemical degrees of freedom, there are two broad approaches: mechanism reduction, where the chemical mechanism is modified to reduce the number of species and the sti ness, and state-space parameterization, where the state of the system is assumed to be parame-terized by a small number of variables which are evolved in the CFD. The techniques proposed in this paper fall into the second category: they seek to obtain a set of variables that parameterizes the thermochemical state, and these variables are then evolved in the CFD calculation. There have been numerous e orts to reduce the dimensionality of a combustion process (see, e.g. [2-12] for a few). Flamelet models such as Steady Laminar Flamelet Method (SLFM) [2-4], flamelet-generated manifold (FGM) [5-7, 9] or flamelet-prolongation of ILDM model (FPI) [11-13] are examples of state-space parameterization model. The remainder of this chapter is organized as follows: we first identify the datasets that will be used to evaluate the proposed model in section 2.3. We then review Principal Component Analysis (PCA) as a technique to obtain a reduced parameter set (section 2.4), discuss how PCA can be formulated as a predictive model (section 2.4.1.2), and introduce 11 adaptive regression to enable parameterization of nonlinear functions of the principal com-ponents (section 2.4.2). Section 2.5 then examines the model in the context of turbulent closure and shows that the model is closed, i.e. it requires no explicit closure model for the thermochemistry. Finally, conclusions are presented in section 2.7. 2.3 Datasets In the discussions below we will consider two datasets: 1. A dataset from a One-Dimensional Turbulence (ODT) simulation which has been done on a temporally evolving nonpremixed CO/H2-air jet with extinction and reig-nition [14, 15]. This was shown to be a statistically accurate representation of a corre-sponding high-fidelity DNS dataset [14]. The calculations include detailed chemical kinetics, thermodynamics, and transport and exhibit significant local extinction and reignition and the dataset is, therefore, a modeling challenge. The state variables are: temperature and species mass fractions for H2, O2, O, OH, H2O, H, HO2, CO, CO2, HCO and N2. 2. Sandia TNF CH4/air Flame D [16]. This flame does not exhibit significant amounts of extinction or reignition, and is a standard modeling target flame. The state variables are temperature and mass fractions for O2, N2, H2, H2O, CH4, CO, CO2, OH and NO. The Flame D dataset is "incomplete" in that it does not contain species reaction rates or a complete set of species. The ODT/DNS dataset, on the other hand, is "complete" in that it has the full set of species, reaction rates, etc. resolved in space and time, but relies on simulation to obtain the data, and is only as accurate as the thermodynamic, kinetic and transport properties that were used in the simulation. When comparing against the datasets, we report R2 values to measure the accuracy with which the model represents the original data, R2 = 1 􀀀 PN i=1( i 􀀀 i )2 PN i=1( i 􀀀 h i)2 ; (2.1) where i is the observed value, i is the predicted value, and h i is the mean value of : For the PCA analysis, we consider data sampled from all space and time in the ODT dataset, 12 and the full dataset for the TNF data. In other words, the PCA does not vary in ~x or t since we sample all ~x and t to obtain the PCA. 2.4 Parameterization Using Principal Component Analysis 2.4.1 Principal Component Analysis Principal Component Analysis (PCA) provides a robust methodology to reduce the number of species equations by identifying correlations in state-space and defining new variables (principal components) that are linear combinations of the original variables (state variables) [17-20]. Details of the formulation have been published elsewhere [19-22], and here we only review the concepts behind the PCA analysis. The basic process of a PCA reduction is 1. Identify a new basis in the multidimensional dataset that is a rotation of the original basis. We call this new basis and the original data : The new basis is obtained via an eigenvalue decomposition of the correlation matrix for many observations of for a system. At this point, we have only performed a rotation, and no information loss has occurred. 2. Truncate the new basis and project the data onto the new basis to obtain an approxi-mation (compression) of the data on the new basis. 3. Given an observation in the truncated basis, we can approximate the value of the original data. This "reconstruction" is a linear reconstruction and is thus very e - cient. Steps 1 and 2 are illustrated conceptually in Figure 2.1. The PCA modeling approach thus requires "training" data which should (ideally) be observations of the a system at conditions close to where we wish to apply the model. Once i is known, the original state variables (e.g. T, yj) can be easily obtained. Furthermore, the accuracy of the parameterization is obtained a priori, and can be adjusted to obtain arbitrary accuracy by increasing the number of retained PCs. This is illustrated in Figure 2.2, where the eigenvalues (relative importance of a given PC in representing the data) as well as the percent of the variance in the original variables are shown as a function of the 13 η1 η2 !2 !1 !2 !1 η1 truncation Figure 2.1: Illustration of the principal components of a hypothetical 2D dataset where we retain a single principal component. Figure 2.2: Eigenvalue magnitude (left axis, bars) and percent variance captured (right axis, line) by retaining the given number of components. number of retained eigenvalues. The eigenvalues can assist in determining how many PCs should be retained to maintain a desired level of accuracy in the resulting model. 2.4.1.1 E ects of scaling Prior to applying PCA, the original data should be centered and scaled [19, 20, 23-25]. There are many di erent scaling options, some of which are enumerated in Table 2.1. Further details regarding scaling may be found in the aforementioned sources. For the purposes of this paper, it is su cient to recognize that the choice of scaling a ects the accuracy of the resulting PCA parameterization. To illustrate this point, consider the results shown in Table 2.2, where R2 values are shown for parameterization of the original state variables by three PCs for various choices of scaling. The e ects of scaling will become even more pronounced when source terms are considered in section 2.4.2.1. However, from Table 2.2, it is evident that scaling can have an appreciable impact in the accuracy of a PCA reconstruction. Unless explicitly stated otherwise, the results presented in this paper were obtained with 14 Table 2.1: Brief descriptions of the scaling options considered here. Method Factor used to scale each state variable STD Standard deviation VAST Ratio of the standard deviation to the mean. Range Maximum-minimum Level Mean Max Maximum Pareto Square-root of the standard deviation 15 Table 2.2: R2 values for PCA projection of state variables with di erent scaling methods using 2 PCs on the temporal CO/H2 dataset. Scaling T H2 O2 O OH H2O H HO2 CO CO2 HCO Average VAST 0.999 0.952 1.000 0.777 0.853 0.954 0.822 0.075 0.996 0.985 0.893 0.846 STD 0.978 0.939 0.995 0.862 0.920 0.927 0.841 0.056 0.988 0.984 0.911 0.855 Level 0.944 0.947 0.978 0.906 0.951 0.877 0.824 0.037 0.976 0.965 0.933 0.849 Range 0.987 0.946 0.999 0.817 0.881 0.952 0.862 0.059 0.996 0.985 0.876 0.851 Max 0.984 0.945 0.999 0.820 0.884 0.950 0.866 0.057 0.996 0.984 0.875 0.851 Pareto 1.000 0.948 0.998 0.746 0.832 0.956 0.809 0.085 1.000 0.981 0.864 0.838 16 VAST scaling. 2.4.1.2 Transport equations for PCs The governing equations can be written as @ @t = 􀀀r ~u 􀀀 r ~J + S ; (2.2) where = 1; ~u; yi; T (or any suitable energy variable in place of T), ~u is the mass-averaged velocity, ~J is the di usive flux of , and S is the volumetric rate of production of . Due to the large number of species present in combustion, Eq. (2.2) represents a large number of strongly coupled partial di erential equations that must be solved. The thermochemical state variables (T, p and ns􀀀1 species mass fractions Yi) define an (ns+1)-dimensional state space which is widely recognized to have lower-dimensional attractive manifolds [26]. Since PCA is a linear transformation, we may apply it directly to the subset of Eq. (2.2) associated with T and yi to derive the transport equations for the PCs. The full derivation has been presented elsewhere [21], and results in @ i @t = 􀀀r i~u 􀀀 r ~J i + S i : (2.3) The source term for the PCs, S i ; is a linear combination of the original (scaled) species and temperature source terms, and must be parameterized in terms of to close the model. It is important to note that (2.3) requires that the PCA definition is independent of space and time so that commutativity with di erential operators is maintained. This can be achieved by using data from all space and time in constructing the PCA reduction, and all analyses presented herein adhere to this principle. In previous work where this approach was originally proposed [21], preliminary results were shown where PCA was performed locally in mixture fraction space (i.e. conditioned on mixture fraction). Here we consider unconditional PCA, and extend the analysis to examine: 1) the e ects of scaling (see section 2.4.1.1) on the source term parameterization, 2) the e ects of filtering on the accuracy of the source term parameterization, and 3) multivariate regression, which will be discussed in section 2.4.2. Just as the species source terms are highly nonlinear functions of yi and T; the PC source terms (S i ) are highly nonlinear functions of the PCs. The original state variables are well-parameterized by the PCs (given a su cient number of retained PCs) because this 17 is the objective of PCA: to identify correlations in the original variables. However, the PCA transformation does not necessarily identify the ideal basis for representing source terms. Furthermore, although the original state variables can be well-characterized by linear functions of , the same is not necessarily true for S i : Thus, several questions remain to be addressed relative to a model based on PCA: 1. Can the truncated basis (see section 2.4.1) adequately represent the PC source terms? 2. Given that the relationship between i and S j is highly nonlinear, can an adaptive regression technique be employed to obtain the functions S j = Fj( 1; 2; : : : n ) (2.4) for the j = 1 : : : n retained PCs? 3. Are the functions represented by Eq. (2.4) sensitive to filtering? In other words, is Fj a function of the filter width, ? This paper aims to address these questions using a priori analysis of high-fidelity combustion data. We next turn our attention to question 2 and outline a methodology to obtain Fj. 2.4.2 Multivariate adaptive nonlinear regression Because we have no physical insight into the appropriate basis functions to form Fj in Eq. (2.4), we need an adaptive method. Multivariate Adaptive Regression Splines (MARS) [27-29] is a technique that allows adaptive selection of basis functions to obtain nonlinear functions such as Fj. At each iteration of the MARS algorithm, a basis function is selected that results in the largest reduction in the regression error. The iterative procedure is re-peated until convergence is achieved. To avoid over-fitting the data, we choose lower-order basis functions (typically quadratic or cubic at most) and subdivide the high-dimensional space into only a few sub-spaces to fit the data (5 sub-spaces were used for the results presented here). Table 2.3 shows the results of applying MARS to map the state variables onto the PCs. Comparing Table 2.3 to Table 2.2, where the state variables were mapped onto the PCs directly via the (linear) PCA transformation, we note an increase in the accuracy of all state 18 Table 2.3: R2 values for MARS regression of state variables on principal components with di erent scaling methods in PCA using 2 PCs on the temporal CO/H2 dataset [14]. Scaling T H2 O2 O OH H2O H HO2 CO CO2 HCO Average VAST 1.000 0.996 1.000 0.996 0.994 0.997 0.989 0.937 0.999 0.998 0.998 0.991 STD 0.995 0.996 0.999 0.986 0.995 0.994 0.997 0.938 0.999 0.991 0.997 0.990 Level 0.990 0.996 0.997 0.982 0.991 0.991 0.994 0.938 0.997 0.982 0.997 0.987 Range 0.997 0.996 1.000 0.991 0.996 0.995 0.993 0.925 0.999 0.993 0.996 0.989 Max 0.996 0.996 1.000 0.989 0.995 0.995 0.994 0.924 0.999 0.992 0.996 0.989 Pareto 1.000 0.993 1.000 0.987 0.984 0.997 0.985 0.921 1.000 0.997 0.953 0.983 19 variables (but particularly minor species and most notably HO2), indicating a nonlinear relationship between the state variables and PCs. This nonlinear relationship has also been observed elsewhere [19, 20, 22], but the MARS approach allows us to capture the nonlinearity between and quite well. Figure 2.3 shows the OH mass fraction, YOH projected into the two-dimensional space defined by the first two principal components, ( 1; 2) : Also shown is a reconstruction of YOH using the (linear) PCA reconstruction (Figure 2.3a) and the nonlinear MARS recon-struction (Figure 2.3b). This clearly illustrates the advantages of the nonlinear reconstruc-tion. 2.4.2.1 MARS for parameterizing PC source terms In contrast to the state variables themselves, where the PCA defines a linear relationship with the PCs, the PC source terms have no linear relationship to the PCs, and adaptive regression is the only plausible method to obtain Fj in Eq. (2.4). Table 2.4 shows the R2 values for the regression of the source terms for various scaling approaches. Notably, that there is a much more significant influence of the choice of scaling on the accuracy with which the PC source terms can be represented than for the state variables (shown in Tables 2.2 and 2.3). With regard to question 1 posed in section 2.4.1.2 (can the S i be parameterized by −2 0 2 4 6 0 −2 4 2 −1 0 1 2 3 4 x 10−3 h1 R2=0.853 h2 Y OH Original PCA (a) −2 0 2 4 6 −2 0 2 4 −1 0 1 2 3 4 x 10−3 h1 R2=0.994 h2 Y OH Original PCA (b) Figure 2.3: Comparison of PCA and MARS reconstructions for OH mass fraction for a two-dimensional model based on principal components 1 and 2. VAST scaling was used. (a) PCA (linear) reconstruction of YOH in ( 1; 2)-space. (b) MARS (nonlinear) reconstruction of YOH in ( 1; 2)-space. 20 Table 2.4: R2 values for MARS regression of source terms on principal components with di erent scaling methods in PCA using 3 PCs on the temporal CO/H2 dataset [14]. Number of PCs 1 2 3 Scaling method S 1 S 1 S 2 Average S 1 S 2 S 3 Average VAST 0.838 0.949 0.929 0.939 0.968 0.938 0.223 0.710 STD 0.041 0.276 0.491 0.383 0.349 0.535 0.183 0.356 Level 0.073 0.331 0.509 0.420 0.437 0.600 0.178 0.405 Range 0.369 0.603 0.551 0.577 0.698 0.661 0.291 0.550 Max 0.407 0.669 0.619 0.644 0.751 0.735 0.253 0.580 Pareto 0.877 0.960 0.956 0.958 0.966 0.963 0.973 0.967 ?) we observe that the source terms have more error in their representation than the original state variables. This suggests that the basis selected by the PCA, which seeks to identify correlations among the state variables, may not be optimal for the representation of the PC source terms. Therefore, other methods that identify a basis that simultaneously optimizes parameterization of both the state variables and the PC source terms should be explored. Nevertheless, as the number of retained PCs increases, the accuracy of the S i parameterization also increases. We should note that the definition for S i remains unchanged as n increases, i.e. S 1 for n = 1 is defined in the same manner as S 1 for n = 3: However, their definitions are di erent for di erent scaling methods. 2.5 Filtering and Turbulent Closure The techniques and results presented in section 2.4 were discussed in the context of fully-resolved quantities. For filtered/averaged quantities, several additional issues arise: 1. How sensitive is the PCA mapping to filtering? In other words, is the PCA mapping itself a ected by filtering? 2. Are the source term functions valid for filtered quantities, i.e., is Fi( 1; 2; : : : n) = Fi(¯ 1; ¯ 2; : : : ¯ n)? We consider each of these issues in the following sections. 21 2.5.1 PCA sensitivity to filtering To determine the sensitivity of PCA to filtering, we examine computational data from a fully resolved CO/H2 jet flame (see section 2.3). The data are filtered and a PCA is applied to the filtered variables. This is performed using a top-hat filter for several filter widths to determine if/how the PCA structure itself is a ected by filtering. Figure 2.4 shows the temperature field extracted along a line-of-sight and shows the e ect of the filter on the temperature profile. x is the grid spacing of the original dataset, whereas refers to the filter width so that = x = 1 implies no filtering. Figure 2.4 indicates that the largest filter width employed here ( = x = 32) has a substantial e ect on the temperature field. Figure 2.5 shows the relative size of the kinetic energy fluctuations, K􀀀¯K ¯K = K0 ¯K , at filter widths of = x = 4 and 16. Note that = x = 16 results in a significant fraction of the kinetic being unresolved, and substantiates the observation from Figure 2.4 that = x = 16 is an appreciable filter width. Figure 2.6 shows the largest five contributions to the first three eigenvectors, which define the rotated basis or the principal components. Consider the first eigenvector. The results indicate that the definition of this eigenvector/PC is almost entirely una ected by filtering. The same results are observed for the second and third eigenvectors. This shows that the PCA reduction itself is insensitive to filtering. The remaining eigenvectors, which are associated with exponentially diminishing eigenvalues (see Figure 2.2), exhibit the same behavior and are not shown for brevity. These results are of significant importance, since the PCA reduction plays a key role in the proposed modeling strategy outlined in section 2.4. 0 0.002 0.004 0.006 0.008 0.01 0.012 500 600 700 800 900 1000 1100 X [m] Temperature [K] D/Dx = 1 D/Dx = 4 D/Dx = 32 Figure 2.4: E ect of filtering on temperature profile for a specific time and realization from the temporal CO/H2 dataset [14]. 22 5 6 7 8 x 10−3 −1 −0.5 0 0.5 1 x (m) K′/ ¯K D/Dx=4 D/Dx=16 Figure 2.5: Instantaneous profile of K􀀀¯K ¯K = K0 ¯K indicating the magnitude of the unresolved kinetic energy at = x = 4 and 16. O2 T H2OCO2 CO 0 0.2 0.4 0.6 0.8 Eigenvector component weight 1st Eigenvector T O2 H2 CO HCO 0 0.2 0.4 0.6 0.8 Eigenvector component weight 2nd Eigenvector HO2 H OH O HCO 0 0.2 0.4 0.6 0.8 1 Eigenvector component weight 3rd Eigenvector D/Dx=1 D/Dx=4 D/Dx=8 D/Dx=16 D/Dx=32 Figure 2.6: Changes in largest (most important) components of the first three eigenvectors with respect to changes in normalized filter width in temporal CO/H2 dataset [14]. The results in Figure 2.6 suggest that, over a substantial range of filter widths, the structure of a PCA remains unchanged. This is an important result since it implies that the definition of a (linear) manifold for the state variables is insensitive to filtering. 2.5.2 Turbulent closure We now turn our attention to the question of whether a mapping = G( ) is valid for the averaged/filtered quantities, i.e. ¯ ?= G( ¯ ): This is particularly important for the source terms that appear in the averaged/filtered PC transport equation, @¯ ˜ @t = 􀀀r ¯ ˜ ~˜u 􀀀 r ~JT + S¯ ; (2.5) where ˜ is the Favre-averaged/filtered value of and ~JT is the turbulent di usive flux. In traditional state-space parameterization approaches, one defines the parameterization variables and then the mapping between the state variables and reaction variables , 23 e.g. = G( ): Then the joint probability density function (PDF) of all , p( ), is used to obtain mean/filtered values of ; ¯ = Z ( ) p( ) d : The problem then becomes how to approximate p( ). If a function = G( ) (2.6) exists so that ¯ = G( ¯ ); (2.7) then there is no turbulent closure problem and the joint PDF of all is not required. 2.5.2.1 Ensemble averaging We first consider ensemble-averaged data from Flame D (see section 2.3). Ensemble averages are formed by number-averaging all samples from a given spatial location. Table 2.5 shows results for all of the species available for flame D. The results show a PCA reconstruction (linear) for two and three retained PCs as well as a (nonlinear) MARS reconstruction based on the same two and three PCs. The "original data" refer to the data processed directly from the flame D dataset where PCA was applied to the entire dataset. PCA and MARS regressions were performed to obtain = G( ) and the resulting R2 values reported. The "ensemble data" used the PCA and MARS regression obtained from the original data and applied it to the ensemble-averaged values for the PCs. Specifically, 1) PCA was applied to the original dataset, 2) MARS was performed to obtain G( ), 3) using the PCA obtained in step 1, the PCs were computed from the original data and then ensemble-averaged to obtain ¯ , 4) ¯ = G( ¯ ) was calculated and compared with the directly averaged values of ¯ to obtain an R2 value. There are several noteworthy points relative to Table 2.5: 1. As n increases from 2 to 3, the R2 value uniformly increases, indicating the increase of accuracy of a PCA-based model as the number of retained components increases. This has been discussed in detail elsewhere [19-22]. 2. The MARS representation of the data is more accurate than the corresponding direct PCA reconstruction, indicating that there is an underlying nonlinear relationship 24 Table 2.5: R2 values for di erent variables in flame D showing that no closure is required to reconstruct the original variables, from the principal components, . Approach Data Type T O2 CO2 NO H2O N2 H2 CH4 CO OH PCA, n = 2 Original 0.990 0.981 0.966 0.860 0.979 1.000 0.385 0.914 0.439 0.495 Ensemble 0.995 0.995 0.987 0.899 0.990 1.000 0.539 0.987 0.611 0.687 PCA, n = 3 Original 0.991 0.991 0.988 0.911 0.989 1.000 0.944 0.916 0.890 0.599 Ensemble 0.996 0.998 0.998 0.939 0.998 1.000 0.982 0.990 0.972 0.650 MARS, n = 2 Original 0.997 0.991 0.976 0.926 0.989 1.000 0.625 0.941 0.666 0.650 Ensemble 0.998 0.999 0.995 0.950 0.999 1.000 0.917 0.992 0.933 0.744 MARS, n = 3 Original 0.997 0.997 0.991 0.974 0.993 1.000 0.938 0.941 0.904 0.745 Ensemble 0.998 0.999 0.999 0.899 0.999 1.000 0.971 0.993 0.979 0.755 25 between the and that the linear PCA-based reconstruction cannot accurately capture. 3. The ensemble-averaged data shows R2 values that are nearly always higher than their corresponding original data values. This suggests that the PCA based models = G( ) do not incur any additional error when evaluated using mean values, ¯ = G( ¯ ): This is true for the linear reconstruction as well as the nonlinear (MARS) reconstruction. 2.5.2.2 Spatial filtering We next consider spatial filtering with the CO/H2 dataset discussed in section 2.3. Figure 2.4 illustrates the e ect of di erent filter lengths on an extracted line-of-sight repre-sented by the ODT data for the temporal CO/H2 dataset. For this particular dataset, a filter width of = x = 16 induces substantial filtering on the data. First, a PCA was performed on the fully resolved data and either n = 2 or n = 3 PCs were retained. This provides a linear mapping between the PCs ( ) and the state variables ( ). Using this mapping, we then compute ¯ directly from the dataset and then use the mapping to approximate ¯ , which is then compared to ¯ calculated directly from the dataset. These results are shown in Figure 2.7. Finally, a MARS regression was performed to map the original variables onto the PCs at the fully resolved scale, providing i = Gi( ): Then, ¯ was calculated directly from the data and ¯ i was approximated as ¯ i = Gi( ¯ ); and this was compared to the value of ¯ i calculated directly from the dataset. The profiles in Figure 2.7 show these results. From Figure 2.7, it is apparent that the error is well-controlled as the filter width is increased, indicating that the PCA-based models require no explicit closure. This is consistent with the results for the ensemble-averaged analysis performed in section 2.5.2.1. Figure 2.8 shows extracted spatial profiles (over a small portion of the domain corre-sponding to an active flame region) for CO2 and OH mass fractions for two di erent filter widths ( = x = 1 and 16) and n = 2. The solid lines represent the profiles extracted directly from the data, whereas the dashed lines are the reconstructed profiles using the PCA/MARS model. These results demonstrate the ability of the PCA/MARS modeling approach to reconstruct the unfiltered and filtered quantities, and also indicate the strong 26 10 20 30 0.9 0.95 1 D/Dx R2 Temperature PCA nh=2 PCA nh=3 MARS nh=2 MARS nh=3 10 20 30 0.9 0.95 1 D/Dx R2 YOH 10 20 30 0.9 0.95 1 D/Dx R2 YCO 2 Figure 2.7: R2 value changes with respect to the changes in normalized filter width (normalized with grid spacing length) for several variables in temporal CO/H2 dataset [14]. 7 7.5 8 8.5 9 9.5 x 10−3 0 0.05 0.1 0.15 x (m) Y CO2 Original: D/Dx = 1 Original: D/Dx = 16 Reconst: D/Dx = 1 Reconst: D/Dx = 16 7 7.5 8 8.5 9 9.5 x 10−3 0 0.5 1 1.5 2 x 10−3 x (m) Y OH Original: D/Dx = 1 Original: D/Dx = 16 Reconst: D/Dx = 1 Reconst: D/Dx = 16 Figure 2.8: Original (solid) and reconstructed (dashed) profiles for CO2 and OH for no filtering and a filter width of = x = 16 and n = 2. filtering that is occurring at = x = 16. It is particularly remarkable that the OH profiles are reconstructed so well by a two-parameter model, and that the filtered profiles are also reconstructed with reasonable accuracy. 2.5.2.3 Source term parameterization We now turn our attention to the parameterization of the PCA source terms in Eq. (2.3) and Eq. (2.5), and seek to answer question 2 posed in section 2.4.1.2 and question 2 in section 2.5: can a function S i = Fi( ) be found, and is S¯ i = Fi( ¯ )? To ascertain the performance of the PCA-based model in representing S¯ i = Fi( ¯ ), we first calculate S i and then obtain the regressing function S i = Fi( ) via MARS. Next, ¯ and S¯ i are calculated directly from the data, and compared against Fi( ¯ ). Figure 2.9 illustrates the results of this in state space while Figure 2.10 shows the associated R2 values. 27 (a) x = 1 (b) x = 4 (c) x = 8 (d) x = 16 Figure 2.9: S¯ 1 (¯ 1; ¯ 2) obtained directly from the data (points) as well as the prediction based on PCA/MARS (surface) for various filter widths using the temporal CO/H2 dataset. Figure 2.10 shows the corresponding R2 values. From these results, as well as those previously presented in Table 2.4, several conclu-sions may be drawn: 1. S i is parameterized with less accuracy than i. This is not surprising given that the PCA was designed to parameterize well, and it is well-known that S i is a highly nonlinear function of so that S i will also be a highly nonlinear function of : 2. The error in the approximation S¯ i = Fi( ¯ ) is bounded and well behaved with the moderate range of filter widths considered in this study. Indeed, the structure of S¯ i ( ) is largely una ected by filtering as shown in Figure 2.9 for S¯ 1 (¯ 1; ¯ 2), and more quantitatively in Figure 2.10. 28 1 4 8 16 32 0.92 0.94 0.96 D/Dx R2 Sh 1 MARS nh = 2 MARS nh = 3 1 4 8 16 32 0.92 0.94 0.96 D/Dx R2 Sh 2 Figure 2.10: R2 value changes with respect to the changes in normalized filter width, = x for the source term of the first and the second PCs, in temporal CO/H2 dataset [14]. Figure 2.11 shows extracted spatial profiles for S 1 and S 2 for a model with n = 2 and filter widths of ( = x = 1 and 16). These results were obtained using Pareto scaling, as this provided the best reconstruction of the source terms as shown in Table 2.4. The solid lines represent the profiles extracted directly from the data, whereas the dashed lines are the reconstructed profiles using the PCA/MARS model. These results correspond to the R2 values reported in Figure 2.10 at = x = 16. While the general trend for S is captured in both cases, it is clear that the detailed profiles for S are not captured fully, and this is also reflected in the relatively low R2 values shown in Table 2.4. Another interesting feature of Figure 2.11 is the structure of S 1 and S 2 are very similar, although their magnitudes are di erent. Figure 2.12 shows the weights that define the 7 7.5 8 8.5 9 9.5 x 10−3 −16 −14 −12 −10 −8 −6 −4 −2 0 x 105 x (m) sh1 Original: D/Dx = 1 Original: D/Dx = 16 Reconst: D/Dx = 1 Reconst: D/Dx = 16 7 7.5 8 8.5 9 9.5 x 10−3 −2.5 −2 −1.5 −1 −0.5 0 0.5 x 104 x (m) sh2 Original: D/Dx = 1 Original: D/Dx = 16 Reconst: D/Dx = 1 Reconst: D/Dx = 16 Figure 2.11: Original (solid) and reconstructed (dashed) profiles for S 1 and S 2 for no filtering and a filter width of = x = 16 and n = 2. 29 T H2 O2 O OH H2O H HO2 CO CO2HCO −1 −0.5 0 0.5 1 Weights 1st EigVec 2nd EigVec 3rd EigVec Figure 2.12: Weights for the first three eigenvectors (which define the first three PCs) for Pareto scaling. first three PCs using Pareto scaling. The first PC ( 1) is defined almost exclusively by temperature (this is commonly the case when Pareto scaling is used). The second PC is defined primarily by the reactants CO, H2 and O2. Therefore, S 1 S T while S 2 is primarily comprised of S CO, S H2 and S O2 . Since these reaction rates are all spatially correlated, it is not surprising that S 1 and S 2 are also highly correlated spatially. As expected, the profiles for 1 and 2 are quite di erent from one another since 1 follows the temperature profile (which peaks in the reaction zone) whereas 2 follows the reactants (which are strongly depleted in the reaction zone). It is important to note that a di erent choice of scaling (resulting in a di erent PC structure) will have a major influence on the resulting source term profiles. This needs to be investigated further and will be the subject of future research. 2.5.3 Comparison with SLFM To provide a reference point with a very common combustion modeling approach, we compare the SLFM model with the PCA model in their respective abilities to reproduce the CO/H2 dataset described in section 2.3. At the outset, we note that the SLFM model (for the purposes of this paper) is parameterized by three parameters: the averaged/filtered mixture fraction (¯Z), its variance ( 02 Z ), and the scalar dissipation rate ( ). The mapping = G(Z; 02 Z ; ) is obtained through solving the flamelet equations [2] and convoluting them with a -PDF for the mixture fraction. Formally, this implies that Z and are statistically independent and that the PDF of the scalar dissipation rate is approximated as ( ¯ ). This is a common modeling assumption, and may be justified in part by observations in previous 30 DNS studies that suggest errors in (Z; ) overshadow errors in approximating p( ) = ( ¯ ) [8]. Figure 2.13 shows parity plots and R2 values comparing the observed and reconstructed values of T and OH for the SLFM and PCA/MARS models at the fully resolved scale (i.e. no filtering). In this case, the PCA/MARS model employed two PCs, consistent with the number of parameters in the SLFM model (Z; ). While one would not expect the SLFM model to perform well in this case since extinction and reignition are present, this does demonstrate the ability of the PCA/MARS model to identify and parameterize the state space e ectively. Figure 2.14 shows realizations of the original and reconstructed data for the PCA/- MARS and SLFM models in state space. Here it is much more clear that the PCA/MARS identifies a manifold in state space whereas the SLFM model does not. Again, while the (a) MARS with n =2 (R2=1.000) (b) MARS with n =2 (R2=0.994) (c) SLFM (R2=0.456) (d) SLFM (R2=0.055) Figure 2.13: Parity plots for temperature (left) and OH (right) reconstructions using PCA/MARS (top row) and SLFM (bottom row). 31 0 0.5 1 10−2 100 102 104 500 1000 1500 2000 R2 = 0.456 Z cmax (1/s) Temperature Original data SLFM reconstruction 0 0.5 1 10−2 100 102 104 0 1 2 3 4 x 10−3 R2 = 0.055 cmax Z (1/s) YOH Original data SLFM reconstruction Figure 2.14: Reconstruction of temperature (left) and OH mass fraction (right) for a 2D PCA/MARS model (top) and the SLFM model (bottom) for the temporal ODT dataset [14]. The R2 value of these reconstructions are reported on each plot as well. SLFM model is not expected to perform well in this situation, the comparison is illustrative of the di erences between the models with same number of parameters. When averaging/filtering is applied, the mixture fraction variance is typically intro-duced as an additional parameter, with a presumed PDF for the mixture fraction parame-terized in terms of ¯Z and 02 Z so that the state variables are obtained via = Z (Z; ¯ ) p(Z; ¯Z ; 02 Z ) dZ: For comparison purposes, we explore the performance of several models, summarized in 32 Table 2.6. Figure 2.15 shows the performance of these models for several state variables as a function of the filter width. There are several important observations to be made: The SLFM model error is bounded if the -PDF model is also used, but the error increases (indicated by the decrease of R2 with increasing ) if no closure is made. The PCA based models demonstrate no sensitivity to filtering, requiring no explicit closure. The two-parameter PCA model is more accurate than the three parameter SLFM/PDF model. While SLFM is not expected to accurately capture the thermochemical state of this system that involves extinction and reignition, these results illustrate the degree of accuracy Table 2.6: Summary of models compared in Figure 2.15. Model Parameters Comments PCA, n = 2 ( 1; 2) Linear reconstruction using two PCs PCA, n = 3 ( 1; 2; 3) Linear reconstruction using three PCs PCA/MARS, n = 2 ( 1; 2) Nonlinear reconstruction using two PCs PCA/MARS, n = 3 ( 1; 2; 3) Nonlinear reconstruction using three PCs SLFM ¯Z ; ¯ Reconstruction using SLFM without closure SLFM/PDF ¯Z ; ¯ ; 02 Z Reconstruction using SLFM with a presumed PDF on Z 10 20 30 0 0.5 1 D/Dx R2 Temperature 10 20 30 0 0.5 1 D/Dx R2 Y OH 10 20 30 0 0.5 1 D/Dx R2 Y CO 2 PCA n h =2 PCA n h =3 MARS n h =2 MARS n h =3 SLFM SLFM b−PDF Figure 2.15: R2 value changes with respect to the changes in normalized filter width ( = x) for several state variables comparing PCA and MARS results with SLFM and SLFM- - PDF results for the temporal CO/H2 dataset. See Table 2.6 for more information. 33 that can be obtained using PCA to identify parameterizing variables for use in defining models, and illustrates that proper selection of parameterizing variables can lead to sig-nificant improvements in model accuracy even for a relatively small number parameters, n . 2.6 Considerations for Model Generation Using PCA as the basis for combustion modeling is still a new concept, and there are several outstanding issues regarding its application. In traditional combustion modeling approaches, a canonical reactor configuration is adopted and solved parametrically to obtain a mapping between the state variables and the parameterizing variables . Such models are, therefore, inherently limited by the assumptions inherent in the canonical reactor. Although PCA can be applied to canonical systems such as the flamelet equations [2], we favor using models such as ODT that allow for a wide range of coupling in length and time scales and provide statistical sampling of the state over a wider range of applicable conditions than traditional canonical reactors. PCA is particularly well-suited for application to such systems because it allows "adaptive" selection of the optimal parameterizing variables. However, it remains to be seen how sensitive the PCA is to the chosen canonical reactor. This may not be critically important in the case of using ODT as the model since it can provide reasonable results for combustion systems [14]. Additionally, the sampling density for the data used to obtain the PCA may be an im-portant consideration. For example, the PCA may be influenced by the over-representation of fuel and oxidizer and relatively small number of observations from the flame regions. Identifying biasing from over/under-sampling is not a trivial task and future work will seek to address this issue. 2.7 Conclusions and Future Work This paper discusses a novel state-space parameterization method. It belongs to the same family as other parameterization methods such as equilibrium, steady laminar flamelet (SLFM) [2], flamelet-prolongation of ILDM [11, 12], etc., but extracts the parameterization directly from data rather than presuming a functional form for the parameterization. We use PCA to identify the model parameters and then MARS (a multidimensional adap- 34 tive regression technique) to obtain the functional form between the progress variables (principal components), and the state variables, . For the jet flame considered in this paper, we observe that the structure (definition) of progress variables thus identified is independent of spatial filtering, and that the functional dependency between and is likewise independent of filter width. The same observations hold for Flame D [16] in the context of Reynolds-averaging rather than filtering. To the extent that this is a universal feature of PCA-based models, these results imply that no explicit closure model is required for the thermochemistry. Further investigation into this observation, particularly using data at higher Reynolds numbers, is certainly warranted to corroborate these observations. The "principal components" that form the independent variables to which the state variables are mapped, are not conserved scalars and their source terms, S , must be pa-rameterized as functions of . We have explored using MARS to achieve this parameter-ization, and found reasonably accurate mappings. However, further work is required here to achieve mappings that are su ciently accurate for predictive modeling. A significant finding presented here is that the functional form is independent of filter width so that given S i = Fi( ); S¯ i = Fi( ¯ ): We have also considered the e ects of scaling (preprocessing the data) on the accuracy of the resulting PCA-based models and have shown that there can be significant influence, particularly on the accuracy with which PC source terms, S i , can be represented. For a point of reference, we have compared the PCA-based models with SLFM and demonstrated that the PCA model is able to, with the same number of parameters, achieve significantly higher accuracy than SLFM. Future work will focus on improving the accuracy with which the PCA transformation can parameterize source terms and consider a posteriori analysis of the modeling approach. Also, it remains to be seen how universal a given PCA definition is. There are several factors that influence this, including sampling density in state space, the dataset from which the PCA is obtained, etc. These important issues will be considered in future work. 2.8 Acknowledgements The authors gratefully acknowledge support from the Department of Energy under award number DE-NT0005015 and the National Science Foundation PetaApps project 0904631. 35 2.9 References [1] H. Tennekes and J. Lumley, A First Course in Turbulence. MIT Press, 1972. [2] N. Peters, "Laminar di usion flamelet models in non-premixed turbulent combus-tion," Progress in Energy and Combustion Science, vol. 10, pp. 319-339, 1984. [3] --, "Laminar flamelet concepts in turbulent combustion," Proceedings of the Com-bustion Institute, vol. 24, pp. 1231-1250, 1986. [4] H. Pitsch and N. Peters, "A consistent flamelet formulation for non-premixed com-bustion considering di erential di usion e ects," Combustion and Flame, vol. 114, pp. 26-40, 1998. [5] J. van Oijen and L. de Goey, "Modelling of premixed laminar flames using flamelet-generated manifolds," Combustion Science and Technology, vol. 161, pp. 113-137, 2000. [6] J. van Oijen, "Flamelet-generated manifolds: Development and application to pre-mixed flames," Ph.D. dissertation, Eindhoven University of Technology, 2002. [7] J. van Oijen and L. de Goey, "Modelling of premixed counterflow flames using the flamelet-generated manifold method," Combustion Theory and Modelling, vol. 6, pp. 463-478, 2002. [8] J. Sutherland, "Evaluation of mixing and reaction models for large-eddy simulation of nonpremixed combustion using direct numerical simulation," Ph.D. dissertation, University of Utah, May 2004. [9] H. Bongers, J. van Oijen, L. Sommers, and L. de Goey, "The flamelet generated manifold method applied to steady planar partially premixed counterflow flames," Combustion Science and Technology, vol. 177, no. 12, pp. 2373-2393, 2005. [10] J. Sutherland, P. Smith, and J. Chen, "A quantitative method for a priori evaluation of combustion reaction models," Combustion Theory and Modelling, vol. 11, no. 2, pp. 287-303, 2007. [11] O. Gicquel, N. Darabiha, and D. Th´evenin, "Laminar premixed hydrogen/air coun-terflow flame simulations using flame prolongations of ILDM with di erential di u-sion," Proceedings of the Combustion Institute, vol. 28, pp. 1901-1908, 2000. [12] B. Fiorina, R. Baron, O. Gicquel, D. Thevenin, S. Carpentier, and N. Darabiha, "Mod-elling non-adiabatic partially premixed flames using flame-prolongation of ILDM," Combustion Theory and Modelling, vol. 7, pp. 449-470, 2003. [13] B. Fiorina, O. Gicquel, S. Carpentier, and N. Darabiha, "Validation of the FPI chemistry reduction method for diluted nonadiabatic premixed flames," Combustion Science and Technology, vol. 176, no. 5, pp. 785-797, 2004. 36 [14] N. Punati, J. Sutherland, A. Kerstein, E. Hawkes, and J. Chen, "An evaluation of the one-dimensional turbulence model: Comparison with direct numerical simulations of CO/H2 jets with extinction and reignition," Proceedings of the Combustion Institute, vol. 33, pp. 1515-1522, 2011. [15] E. Hawkes, R. Sankaran, J. Sutherland, and J. Chen, "Scalar mixing in direct numeri-cal simulations of temporally-evolving plane jet flames with detailed co/h2 kinetics," Proceedings of the Combustion Institute, vol. 31, no. 1, pp. 1633-1640, January 2007. [16] "International workshop on measurement and computation of turbulent nonpremixed flames." [Online]. Available: http://www.sandia.gov/TNF/DataArch/FlameD.html [17] I. Jolli e, Principal Component Analysis, 2nd ed., ser. Springer Series in Statistics. New York: Springer, 2002. [18] J. Shlens, "A tutorial on principal component analysis," Systems Neurobiology Lab-oratory, Salk Insitute for Biological Studies La Jolla, CA 92037 and Institute for Nonlinear Science, University of California, San Diego La Jolla, CA 92093-0402, Tech. Rep., December 2005. [19] A. Parente, "Experimental and numerical investigation of advanced systems for hydrogen-based fuel combustion," Ph.D. dissertation, Universit`a di Pisa, 2008. [20] A. Parente, J. Sutherland, L. Tognotti, and P. Smith, "Identification of low-dimensional manifolds in turbulent flames," Proceedings of the Combustion Institute, vol. 32, pp. 1579-1586, 2009. [21] J. Sutherland and A. Parente, "Combustion modeling using principal component analysis," Proceedings of the Combustion Institute, vol. 32, pp. 1563-1570, 2009. [22] A. Parente, J. Sutherland, B. Dally, L. Tognotti, and P. Smith, "Investigation of the MILD combustion regime via principal component analysis," Proceedings of the Combustion Institute, vol. 33, no. 2, pp. 3333-3341, 2011. [Online]. Available: http://dx.doi.org/10.1016/j.proci.2010.05.108 [23] H. Keun, T. Ebbels, H. Antti, M. Bollard, O. Beckonert, E. Holmes, J. Lindon, and J. Nicholson, "Improved analysis of multivariate data by variable stability scaling: application to NMR-based metabolic profiling," Analytica Chimica Acta, vol. 490, no. 1-2, pp. 265-276, August 2003. [24] R. van den Berg, H. Hoefsloot, J. Westerhuis, A. Smilde, and M. van der Werf, "Cen-tering, scaling, and transformations: improving the biological information content of metabolomics data," BMC Genomics, vol. 7, p. 142, 2006. [25] I. Noda, "Scaling techniques to enhance two-dimensional correlation spectra," Jour-nal of Molecular Structure, vol. 883-884, pp. 216-227, July 2008. [26] U. Maas and S. Pope, "Implementation of simplified chemical kinetics based on in-trinsic low-dimensional manifolds," Proceedings of the Combustion Institute, vol. 24, 37 pp. 103-112, 1992. [27] J. Friedman, "Multivariate adaptive regression splines," Annals of Statistics, vol. 19, no. 1, pp. 1-67, 1991. [28] --, "Fast MARS," Stanford University Department of Statistics, Tech. Rep. 110, 1993. [29] J. Friedman and C. Roosen, "An introduction to multivariate adaptive regression splines," Statistical Methods in Medical Research, vol. 4, pp. 197-217, 1995. CHAPTER 3 IDENTIFICATION OF INVARIANT MANIFOLDS IN TURBULENT COMBUSTION SYSTEMS2 3.1 Abstract Principal Component Analysis (PCA) has been demonstrated as an e ective method-ology to identify attractive manifolds in state-space for turbulent combustion systems. This study explores the invariance of the manifolds identified by PCA. Data from several sources involving di erent fuels and a variety of Reynolds numbers are considered. Results indicate that, for given reactant compositions and temperatures, PCA-derived manifolds are invariant with respect to Reynolds number. Additionally, based on one dataset, results indicate relative insensitivity to filtering as well. 3.2 Introduction Turbulent combustion is characterized by a spectrum of length and time scales for both the chemistry and fluid dynamics. At large Dahmkoler number (Da), when these scales are segregated and mixing timescales are much slower than chemical timescales, chemical equilibrium prevails. In situations, the mixing must be modeled (via, e.g., computational fluid dynamics) while the chemical reactions can be assumed to be in equilibrium and a simple state relationship can be used to relate a passive scalar such as the mixture fraction to the state of the system. On the other hand, when mixing timescales are extremely fast relative to chemical timescales (Da ! 0), then the system is "well-mixed" and the problem can be modeled by a well-stirred reactor, focusing on kinetics while ignoring mixing. Manifolds have been investigated from a dynamical systems perspective, where a rate analysis is performed to identify timescale segregation that implies an attractive manifold 2This chapter's work benefits from a jointwork with Alessandro Parente. 39 (where fast timescales relax the system onto a slowly varying manifold) [1]. Alternative methods postulate the manifold a priori as being described by some controlling variables such as the mixture fraction and dissipation rate [2]. One of the major challenges in turbulent combustion lies in modeling the situation where the mixing and chemical timescales overlap. As the degree of overlap increases (Da ! 1), more timescales become coupled and finite-rate chemistry e ects become increasingly important. Most combustion modeling approaches begin at the large Da limit and then attempt to incorporate some coupling between mixing and reaction. The steady laminar flamelet model is such an example that introduces the scalar dissipation rate as a mixing timescale that perturbs the state of the system away from chemical equilibrium [2]. Many variations on this have been proposed to model regimes where there is increased overlap in timescales. These models frequently add a chemical progress variable to account for this additional coupling [3-7]. The definition of progress variables is a major challenge. One would like each addi-tional progress variable (parameterizing variable) to be "orthogonal" to the previous ones so that it captures information not already represented. Additionally, progress variables should be chosen so that they represent as much of the variation in the system (departure from equilibrium) as possible. Principal Component Analysis (PCA) has been shown as a viable technique to identify progress variables for use in identifying manifolds in combustion [8-10]. However, a significant question remains as to whether the structure of the PCA (i.e. the definition of progress variables) is invariant with respect to system parameters such as Da, Re, and filter width. This chapter explores manifold variation for several datasets to determine the universality of the manifolds identified by PCA in combustion data. 3.3 Principal Component Analysis Here we present a brief overview of PCA. More details can be found elsewhere [8-11]. Consider m observations of n variables arranged in an n m matrix X whose columns rep-resent individual observations and rows correspond to di erent variables. PCA determines a basis for the data X such that the data are well-represented by a truncated basis [12, 13]. 40 The covariance matrix is defined by1 R = 1 n􀀀1X>X and the eigenvector decomposition of R may be obtained as = Q􀀀1RQ, where Q are the orthonormal eigenvectors of R, with Q􀀀1 = Q>. The eigenvectors (columns of Q) form a new basis, and the principal components (PCs) of the data in X are defined as = QX. The full set of PCs exactly reproduces all observations in the original data, by definition. The real utility in PCA comes by exploiting the fact that PCA maximizes the variance of the data in each PC direction. The rotated coordinate system has the property that the first dimension (corresponding to the largest eigenvalue) is selected to best represent the variance in the data. Subsequent directions each represent the next-largest variance in the data. Therefore, a truncated basis, i.e. a subset of the columns in Q, can approximate the original data remarkably well. We define a transformation matrix A as a rank-deficient subset of the Q matrix with n rows and n columns. The columns of A correspond to the columns of Q with the n largest eigenvalues. We may then approximate X as X A>: (3.1) This approximation by a reduced set of variables is precisely what the modeling approaches discussed earlier are trying to accomplish. While we will focus on the linear reconstruction given by equation (3.1) here, e orts to employ PCA to identify parameters for the lower-dimensional representation and use nonlinear reconstruction techniques have also been applied with PCA to provide greater accuracy [14]. In the context of combustion applications, the n variables comprising the rows of X are the ns+1 variables [T; p; Y1; Y2; : : : ; Yns􀀀1]. Performing a PCA on this set of variables yields a new (ns + 1)-dimensional basis, , which is a rotation of the original basis. Retaining n < (ns + 1) columns of Q with the largest eigenvalues defines a basis (A) for a n - dimensional parameterization of the thermochemical state of the system. According to the definition of the PCA model, it needs an empirical dataset in order to be trained on. Therefore, before utilizing PCA in actual simulations we need to make sure that it is possible to use one dataset to train the PCA mapping on, and then use it for the simulation of other cases with di erent situations without injecting a significant amount 1Here we have assumed that the data are centered (their mean is zero) and scaled by constant factors i. These are common procedures that can strongly influence the results of the PCA. For results here, we chose standard scaling (see [11, 13]). 41 of error to PCA reconstruction. So we should study the invariant of the PCA structure with respect to di erent physical parameters of the problem. Our objectives here are to determine: 1. Is the structure of the PCA (defined by Q) sensitive to system parameters such as Reynolds number (Re) and Da? 2. How sensitive is the PCA reconstruction to changes in Q? These questions are very important for PCA to be used as the basis of a practical combustion modeling approach. If the manifolds identified by PCA are invariant with respect to system parameters such as Re, then high-fidelity simulations at low Re can be used to extract manifold information for use in simulations at higher Re. Invariance of the manifold implies invariance of the optimal progress variables for the system. 3.4 Results To investigate whether PCA can identify invariant manifolds, we consider data from several sources: CH4/H2 jet flames in hot co-flow [15], CH4/air Sandia flames [16], One-Dimensional Turbulence (ODT) simulations of CO/H2-air jet flames at varying Re. This corresponds to DNS [17] and the ODT model has been shown to reproduce the DNS statistics with reasonable accuracy [18]. Computational data (the ODT dataset) has the advantage of being "complete" in the sense that the full thermochemical state is available and is not subject to measurement error. On the other hand, experimental data have the advantage of being observations of real systems, absent approximations inherent in thermochemical data, etc. that simulation relies on. By considering both, we hope that the conclusions drawn in this chapter will be more reliable. In all of the results shown below, we will use the R2 value (2.1), as a measure of the accuracy with which PCA reconstructs the original data. In this chapter, we are primarily 42 concerned with the variability of these R2 values when PCAs used from systems with di erent reactant compositions are used. 3.4.1 Dependence on chemical composition Here the e ect of changes in reactant composition on the resultant manifold structure is considered. PCA is performed on a jet in hot co-flow (JHC) burner, designed to emulate flameless conditions [15]. It consists of a central fuel jet (80% CH4 and 20% H2) within an annular co-flow of hot exhaust products from a secondary burner mounted upstream of the jet exit plane. Two di erent O2 levels in the co-flow are considered: 3 and 9% (by vol.), while the temperature and exit velocity are kept constant. The variation of the oxygen content in the co-flow from 9 to 3% is used to control the transition from the traditional non-premixed di usion flame regime to flameless conditions [10]. It is, therefore, interesting to investigate the e ect of the combustion regime on the manifold identification to determine whether an invariant manifold can be still identified for the systems. Figure 3.1 shows the PCA structure (the Q matrix) of the first four PCs (first four eigenvectors/columns of Q) for the 3% and 9% datasets, and indicates significant variability in the PC structure. CO component of the first PC, T, and CO2 components of the second PC, O2, and NO components of the third PC and CO component of the forth PC are few examples of these di erences. Table 3.1 summarizes the results for the JHC datasets, when retaining 4 PCs. Table Figure 3.1: Structure of the first four PCs for the JHC data with 3% and 9% O2 in the co-flow, showing the variability of the PC structure.1 1This figure shows the results of a study conducted by Alessandro Parente. 43 Table 3.1: R2 values for the JHC 3% and 9% O2 datasets. R2 3% and R2 9% indicate reconstruc-tion using the PCA obtained from the 3% and 9% O2 cases, respectively.1 R2 R2 9% R2 3% 3% 9% 3% 9% T 0.940 0.956 0.947 0.910 YO2 0.950 0.988 0.949 0.948 YN2 0.988 0.996 0.981 0.981 YH2 0.997 0.997 0.997 0.997 YH2O 0.970 0.959 0.968 0.956 YCH4 0.998 0.998 0.998 0.998 YCO 0.996 0.981 0.780 0.994 YCO2 0.958 0.968 0.966 0.860 YOH 0.999 0.995 0.968 0.998 YNO 0.812 0.947 0.715 0.481 3.1 shows that, in this case, the manifold invariance does not hold. If the A matrix from the 3% case is used, the main species and temperature are correctly recovered; however, the accuracy in CO and NO reconstruction dramatically drops, confirming the sharp character of the transition from the flameless to the traditional flame regime, which mainly a ect the distribution of pollutant and intermediate species such as NO and CO [10]. Similarly, the application of the manifold obtained from the 9% case to the 3% one leads to very unsatisfactory results for NO and CO2 mass species distributions. 3.4.2 Invariance to Reynolds number The invariance of the chemical manifold is demonstrated for a set of four simple jet flames (denoted C - F) with equal fuel composition (25% CH4, 75% air) but di erent Reynolds numbers: Re=13400 (C), 22400 (D), 33600 (E), 44800 (F) [16]. Here we studied the sensitivity of the reconstruction given in equation (3.1) with respect to the variation in the PCA structure. To investigate this, we use a manifold (Q) obtained on one flame to reconstruct all of the other flames via (3.1). The R2 values are obtained keeping three PCs of the covariance matrix (out of 10 original variables). Table 3.2 shows: the R2 values obtained from the PCA obtained directly on a given dataset, 1This table shows the results of a study conducted by Alessandro Parente. 44 Table 3.2: Reconstruction accuracy, R2, for the flames C, D, E and F datasets by the PCA reduction. Note that R2C and R2 F refer to the accuracy by which variables are reconstructed (via 2.1) using the PCA obtained for flame C and F, respectively R2 R2 F R2C C D E F C D E D E F T 0.985 0.985 0.976 0.971 0.982 0.981 0.974 0.984 0.977 0.974 YO2 0.987 0.986 0.980 0.979 0.983 0.983 0.979 0.986 0.980 0.977 YN2 0.982 0.982 0.980 0.980 0.983 0.981 0.980 0.981 0.979 0.979 YH2 0.975 0.969 0.964 0.970 0.973 0.966 0.966 0.966 0.964 0.965 YH2O 0.989 0.989 0.986 0.984 0.987 0.986 0.984 0.988 0.985 0.984 YCH4 0.987 0.987 0.984 0.984 0.986 0.985 0.984 0.986 0.985 0.985 YCO 0.972 0.968 0.962 0.969 0.970 0.983 0.962 0.964 0.962 0.970 YCO2 0.987 0.986 0.976 0.974 0.983 0.980 0.975 0.985 0.977 0.975 YOH 0.999 0.999 0.995 0.978 0.990 0.990 0.984 0.999 0.999 0.998 YNO 0.945 0.932 0.887 0.892 0.942 0.931 0.895 0.933 0.877 0.850 R2C - the R2 values obtained when using the PCA from flame C to reconstruct data in flames D-F, R2 F - the R2 values obtained when using the PCA from flame F to reconstruct data in flames C-E. Table 3.2 shows the e ect of variability in the PCA structure. The important comparison in Table 3.2 is between R2 values of a given variable on a given flame. For example, the temperature for E is reconstructed with R2 = 0:976. If the PCA from C is used for the reconstruction, we obtain R2 = 0:974 and if the PCA from F is used, we obtain R2 = 0:977. From reconstruction result of some of the species like NO and CO, we can see that training the PCA model on a dataset with higher Re has slightly better reconstruction result. This can be because of the more statistical information in the higher Re case and information about extinction and reignition that occur more in higher Re numbers. The variability in the PCA structure will be shown in the next chapter where the invariance to the Re and filter width is investigated simultaneously. Analysis of data from this experimental dataset at various Re indicate that PCA can be used to identify manifolds that are invariant with respect to Re. This conclusion is further supported by the results from computational studies that will be presented below. 45 3.4.3 Invariance to Reynolds number and filter width Finally, we consider the sensitivity of PCA structure to the filter width ( ), a relevant question in the context of LES. Data from three di erent CO/H2 flames at increasing Re (here denoted cases L, M, and H for Re=2510, 4478, and 9079, respectively) are considered. These data were obtained via ODT simulations with detailed thermochemistry [18] and correspond directly to three di erent DNS datasets [17]. By using computational data with detailed chemistry, we have access to the full set of species involved in the chemical mechanism (11 in this case). Figure 3.2 shows the variation in the structure of the first three principal components when varying Re and . The filter widths vary from = 1 (no filtering) to = 32 times the DNS grid spacing. The first two PCs show remarkable invariance to both Re and over the ranges considered here. Note that for this situation, = 32 corresponds to quite aggressive filtering due to the relatively low Re [14]. For the third PC, some variation is observed, particularly in the weighting on HCO. However, the variation is stronger with Re than with (for the range of Re and considered here). Table 3.3 shows the R2 values to illustrate the accuracy of linear reconstruction (via equation (3.1) with three retained eigenvectors) and its sensitivity to the minor variation in the PCA structure that occurs across a range of Re and . R2d irect indicates the base-line R2 value obtained from a PCA performed directly on the data for that Re and . R2 H;1 indicates R2 errors using a PCA obtained on case H (the high Re case) at a filter width Figure 3.2: Largest five contributions to the first three PCs in the PCA for the three Re cases at various filter widths. 46 of = 1 and R2 L;32 indicates R2 errors using a PCA obtained from case L at = 32. Reconstruction is performed using equation (3.1). Results in Table 3.3 indicate that the minor variability in the PCA structure observed across these Re numbers and 's do not impact the reconstruction accuracy to a large degree. For example, consider the recon-struction accuracy for OH in case H with = 8. The three entries in Table 3.3 show R2d irect = 0:947 (the reconstruction accuracy using the PCA obtained directly on case H with = 1), R2 H;1 = 0:947 (reconstruction accuracy using PCA from case H with = 1) and R2 L;32 = 0:913 (reconstruction accuracy using PCA from case L with = 32). Analyzing the data presented in Figure 3.2 and Table 3.3, we can draw a few conclu-sions: The PCA structure (see Figure 3.2) is relatively invariant across the range of Reynolds numbers (Re from 2510 to 9079) and filter widths ( = 1 to 32 times the DNS grid spacing) considered in these data. The third eigenvector shows more variation than the first two, however. The reconstruction error (see Table 3.3) is relatively invariant across Re. For exam-ple, the temperature R2 for case L (low Re) at = 1 is R2d irect = 0:986. If a PCA obtained from case H with = 1 is used, we find R2 H;1 = 0:986. The reconstruction error is relatively invariant across . Specifically, comparing entries in the R2 H;1 or R2 L;32 block to the R2d irect block in Table 3.3, we see that the R2 values are largely invariant when a PCA obtained at one filter width is used to reconstruct data at another filter width. 3.5 Conclusion PCA has been shown to be e ective in identifying low-dimensional manifolds that exist in turbulent combustion [8, 10]. An attractive property of PCA is that the accuracy of a representation can be easily controlled by increasing the number of retained eigenvectors (PCs) in the analysis. This chapter has explored the universality of manifolds identified by PCA on turbulent combustion data. Standard scaling has been used in the data preprocess-ing [11], but we observe invariance using other scaling techniques as well. We considered two experimental datasets and a computational dataset, each of which varied the Reynolds 47 Table 3.3: R2 values across a range of filter widths and Reynolds numbers for the CO/H2 datasets. R2d irect indicates the R2 values for reconstructions from PCA obtained directly on the data at the given eRe and . R2 H;1 indicates reconstruction using a PCA obtained on case H (high Re) with = 1. R2 L;32 indicates a reconstruction using a PCA obtained on case L (low Re) with = 32. R2d irect R2 H;1 R2 L;32 = x L M H L M H L M H T 1 0.986 0.989 0.985 0.986 0.989 0.985 0.984 0.985 0.980 4 0.987 0.990 0.987 0.987 0.989 0.987 0.985 0.987 0.983 8 0.988 0.991 0.988 0.988 0.990 0.988 0.986 0.988 0.984 16 0.989 0.992 0.989 0.989 0.991 0.989 0.988 0.989 0.985 32 0.991 0.993 0.991 0.992 0.992 0.991 0.991 0.991 0.988 YH2 1 0.950 0.964 0.950 0.913 0.952 0.950 0.951 0.957 0.934 4 0.952 0.966 0.953 0.916 0.954 0.953 0.953 0.958 0.936 8 0.954 0.968 0.955 0.918 0.957 0.956 0.954 0.959 0.938 16 0.957 0.972 0.960 0.923 0.961 0.961 0.957 0.962 0.941 32 0.961 0.976 0.966 0.926 0.965 0.967 0.961 0.966 0.946 YO2 1 0.997 0.994 0.993 0.996 0.995 0.993 0.996 0.993 0.991 4 0.997 0.995 0.994 0.997 0.995 0.994 0.996 0.994 0.992 8 0.997 0.996 0.995 0.997 0.996 0.995 0.997 0.995 0.993 16 0.998 0.996 0.996 0.998 0.997 0.996 0.997 0.996 0.994 32 0.998 0.997 0.996 0.999 0.998 0.997 0.998 0.997 0.995 YO 1 0.861 0.932 0.914 0.846 0.934 0.914 0.861 0.887 0.861 4 0.864 0.934 0.917 0.847 0.936 0.917 0.863 0.889 0.866 8 0.867 0.937 0.922 0.849 0.938 0.921 0.867 0.892 0.871 16 0.874 0.943 0.928 0.853 0.943 0.927 0.875 0.898 0.879 32 0.893 0.954 0.941 0.862 0.952 0.939 0.893 0.910 0.895 YOH 1 0.927 0.949 0.939 0.927 0.953 0.939 0.928 0.933 0.901 4 0.929 0.952 0.944 0.929 0.955 0.943 0.930 0.936 0.908 8 0.932 0.955 0.947 0.931 0.957 0.947 0.933 0.940 0.913 16 0.938 0.960 0.953 0.935 0.961 0.952 0.939 0.946 0.921 32 0.951 0.968 0.961 0.945 0.967 0.960 0.951 0.955 0.932 YH2O 1 0.946 0.945 0.939 0.940 0.942 0.939 0.939 0.927 0.898 4 0.947 0.946 0.942 0.941 0.943 0.942 0.942 0.930 0.905 8 0.949 0.948 0.944 0.942 0.944 0.944 0.945 0.933 0.909 16 0.953 0.950 0.946 0.945 0.946 0.947 0.950 0.938 0.916 32 0.959 0.955 0.951 0.951 0.949 0.952 0.959 0.946 0.925 YH 1 0.879 0.915 0.889 0.861 0.913 0.889 0.870 0.824 0.757 4 0.883 0.918 0.894 0.863 0.916 0.894 0.876 0.836 0.777 8 0.887 0.922 0.899 0.866 0.920 0.899 0.883 0.845 0.789 16 0.896 0.929 0.906 0.874 0.927 0.906 0.894 0.859 0.804 32 0.914 0.939 0.917 0.888 0.936 0.917 0.914 0.880 0.828 YHO2 1 0.997 0.981 0.970 0.948 0.965 0.970 0.996 0.997 0.998 4 0.997 0.983 0.974 0.952 0.968 0.974 0.997 0.998 0.998 8 0.997 0.984 0.977 0.954 0.970 0.977 0.997 0.998 0.998 48 Table 3.3 Continued R2d irect R2 H;1 R2 L;32 = x L M H L M H L M H 16 0.997 0.986 0.980 0.957 0.973 0.980 0.997 0.998 0.998 32 0.997 0.988 0.983 0.960 0.976 0.983 0.997 0.998 0.998 YCO 1 0.992 0.987 0.985 0.987 0.987 0.985 0.991 0.984 0.980 4 0.993 0.988 0.986 0.987 0.988 0.986 0.992 0.985 0.982 8 0.993 0.989 0.987 0.988 0.989 0.988 0.993 0.986 0.984 16 0.994 0.991 0.990 0.990 0.991 0.990 0.994 0.988 0.986 32 0.995 0.992 0.992 0.991 0.993 0.993 0.995 0.990 0.988 YCO2 1 0.985 0.989 0.988 0.981 0.989 0.988 0.983 0.980 0.979 4 0.985 0.990 0.990 0.982 0.990 0.990 0.984 0.981 0.980 8 0.986 0.991 0.990 0.982 0.991 0.990 0.985 0.981 0.981 16 0.986 0.991 0.991 0.983 0.991 0.991 0.986 0.982 0.982 32 0.988 0.992 0.992 0.984 0.992 0.993 0.988 0.983 0.984 YHCO 1 0.971 0.941 0.907 0.958 0.940 0.907 0.971 0.923 0.835 4 0.974 0.946 0.917 0.962 0.944 0.915 0.975 0.927 0.840 8 0.978 0.953 0.927 0.966 0.951 0.924 0.979 0.933 0.845 16 0.984 0.965 0.942 0.972 0.963 0.938 0.984 0.944 0.855 32 0.988 0.975 0.958 0.976 0.972 0.953 0.988 0.952 0.864 49 number of the flow (and, by implication, the Dahmkoller number). The computational dataset a orded the opportunity to explore e ects of filter width on the identified manifold. The primary conclusions obtained suggest: The PCA structure is relatively invariant across the range of Reynolds numbers (Re from 2510 to 9079 and from 13400 to 44800). PCA can identify manifolds that are relatively insensitive to filter width. Manifolds obtained via PCA are not invariant with respect to the pure stream states (temperature, composition). The structure of the manifold was examined by looking at the contributions to the pri-mary eigenvectors determined by PCA, which correspond to the "directions" in state space that exhibit the largest variance. These findings have important implications for PCA-based modeling strategies [9, 14], since PCA must be performed on an initial empirical dataset to construct the model. The findings here suggest that the dataset may not need to be very close in state space to the intended application to remain valid. 3.6 References [1] U. Maas and S. Pope, "Laminar flame calculations using simplified chemical kinetics based on intrinsic low-dimensional manifoldssic low-dimensional manifolds," Pro-ceedings of the Combustion Institute, vol. 25, pp. 1349-1356, 1994. [2] N. Peters, "Laminar di usion flamelet models in non-premixed turbulent combus-tion," Progress in Energy and Combustion Science, vol. 10, pp. 319-339, 1984. [3] J. van Oijen and L. de Goey, "Modelling of premixed laminar flames using flamelet-generated manifolds," Combustion Science and Technology, vol. 161, pp. 113-137, 2000. [4] B. Cuenot, F. Egolfopoulos, and T. Poinsot, "An unsteady flamelet model for non-premixed combustion," Combustion Theory and Modelling, vol. 4, pp. 77-97, 2000. [5] C. Pierce, "Progress-variable approach for large-eddy simulation of turbulent com-bustion," Ph.D. dissertation, Stanford University, 2001. [6] B. Fiorina, R. Baron, O. Gicquel, D. Thevenin, S. Carpentier, and N. Darabiha, "Mod-elling non-adiabatic partially premixed flames using flame-prolongation of ILDM," Combustion Theory and Modelling, vol. 7, pp. 449-470, 2003. 50 [7] J. Sutherland, P. Smith, and J. Chen, "A quantitative method for a priori evaluation of combustion reaction models," Combustion Theory and Modelling, vol. 11, no. 2, pp. 287-303, 2007. [8] A. Parente, J. Sutherland, L. Tognotti, and P. Smith, "Identification of low-dimensional manifolds in turbulent flames," Proceedings of the Combustion Institute, vol. 32, pp. 1579-1586, 2009. [9] J. Sutherland and A. Parente, "Combustion modeling using principal component analysis," Proceedings of the Combustion Institute, vol. 32, pp. 1563-1570, 2009. [10] A. Parente, J. Sutherland, B. Dally, L. Tognotti, and P. Smith, "Investigation of the MILD combustion regime via principal component analysis," Proceedings of the Combustion Institute, vol. 33, no. 2, pp. 3333-3341, 2011. [Online]. Available: http://dx.doi.org/10.1016/j.proci.2010.05.108 [11] A. Parente, "Experimental and numerical investigation of advanced systems for hydrogen-based fuel combustion," Ph.D. dissertation, Universit`a di Pisa, 2008. [12] J. Jackson, A User's Guide to Principal Component Analysis, ser. Wiley Series in Probability and Statistics. Wiley-Interscience, 1991. [13] I. Jolli e, Principal Component Analysis, 2nd ed., ser. Springer Series in Statistics. New York: Springer, 2002. [14] A. Biglari and J. Sutherland, "A filter-independent model identification technique for turbulent combustion modeling," Combustion and Flame, vol. 159, no. 5, pp. 1960- 1970, 2012. [15] B. Dally, A. Karpetis, and R. Barlow, "Structure of turbulent nonpremixed jet flames in a diluted hot coflow," Proceedings of the Combustion Institute, vol. 29, pp. 1147- 1154, 2002. [16] R. Barlow and J. Frank, "Piloted ch4/air flames c, d, e, and f - release 2.1," Sandia National Laboratories, Livermore, CA 94551-0969, Tech. Rep., Jun 2007. [17] E. Hawkes, R. Sankaran, J. Chen, S. Kaiser, and J. Frank, "An analysis of lower-dimensional approximations to the scalar dissipation rate using direct numerical simulations of plane jet flames," Proceedings of the Combustion Institute, vol. 32, no. 1, pp. 1455-1463, 2009. [18] N. Punati, J. Sutherland, A. Kerstein, E. Hawkes, and J. Chen, "An evaluation of the one-dimensional turbulence model: Comparison with direct numerical simulations of CO/H2 jets with extinction and reignition," Proceedings of the Combustion Institute, vol. 33, pp. 1515-1522, 2011. CHAPTER 4 AN A-POSTERIORI EVALUATION OF PRINCIPAL COMPONENT ANALYSIS-BASED MODELS FOR TURBULENT COMBUSTION SIMULATIONS2 4.1 Abstract Recently, Principal Component Analysis (PCA) has been proposed as a means to iden-tify and parameterize manifolds existing in turbulent reacting flow. PCA provides a way to systematically add dimensionality to state-space parameterizations without resorting to ad-hoc progress variables. This work presents an a posteriori study of PCA as a combustion model applied to a nonpremixed CO/H2 temporally evolving jet flame with extinction and reignition. As a basis for comparison, results from detailed chemistry calculations are compared with the PCA-transport results to verify the model and evaluate its performance. The e ect of increasing the number of retained principal components (PCs) is shown by comparing the results of three di erent cases retaining 1, 2, and 3 PCs. Invariance of the model's error to a system parameter (Reynolds number) is evaluated over a range of Re numbers from 2500 to 9000. 4.2 Introduction There are two general approaches to reduce the state-space dimensionality in combus-tion problems. The first approach employ reduction techniques such as the Quasi Steady State Approximation (QSSA) [1], or Rate Controlled Constrained-Equilibrium (RCCE) [2, 3]. The second approach attempts to identify controlling parameters and a function that maps the state of the system to these parameters. 2The material presented in this chapter has been submitted to Combustion and Flame. 52 Traditional state-space parameterization methods postulate some of the variables in the system as controlling variables and then parameterize the state variables based on them; therefore, the only equations to be solved are the transport equations for a new smaller set of controlling variables. The Steady Laminar Flamelet Method (SLFM) [4] is an example of these methods using mixture fraction and dissipation rate as the controlling variables. One disadvantage of most parameterization methods is that they are not easily extensible; adding additional parameters is not straightforward. Techniques such as Flamelet-Generated Man-ifold (FGM) [5-7] and Flamelet-Prolongation of ILDM model (FPI) [8-10] propose ad-hoc progress variables, but su er from the same problem of identifying the ‘best' set of parameters. As additional parameters are added to a model, these should be orthogonal to one another to optimally represent additional manifold dimensions. Additionally, these parameters should represent as much of the variation in the system as possible in order to maximize the e ectiveness of each parameterizing variable. The aforementioned criteria have guided recent work on dimension reduction tech-niques using Principal Component Analysis (PCA). PCA has been shown as an e ective methodology to identify the lower manifolds in turbulent combustion simulations [11-13]. One of the distinguishing features of a PCA-based combustion model is that one can select the number of parameters to achieve a desired accuracy of the model, and that each additional parameter explains the most remaining variance in the original data as possible [11, 13, 14]. There have been several a priori studies on PCA-based models for turbulent combustion simulations. One key outcome of these studies is the recognition that, while PCA provides a basis (a linear rotation of the original state) and a means to truncate this basis (providing the optimal low-dimensional representation), the state variables are still nonlinear functions of the new basis, and nonlinear regression is required [14-18]. Additionally, nonlinear parameterization methods based on PCA have been explored [12, 13, 19-21]. Most of these, however, result in a reduced representation that e ectively precludes incorporation in a CFD framework since the transport equations for the resulting parameters are very complex. A noteworthy exception is the Manifold-Generated Local PCA approach proposed by Coussement et al. [19]. Relatively few a posteriori studies of PCA-based models have been undertaken. Several 53 have been done based on variations of the original formulation proposed in [11] including work by Coussement et al. [19] Isaac et al. [22]. Recent work by Mirgolbabaei and Echekki [23] showed an a posteriori demonstration based directly on the originally proposed for-mulation [11]. In the present work, we use the approaches outlined in [11, 14, 24] to demonstrate a PCA-based turbulent combustion model which transports the PCs directly. The model is trained on simulations with detailed chemistry involving 11 species, and we