Development of an Open-Source Autonomous Computational Fluid Dynamics Meta-Modeling Environment for Small-Scale Combustor Optimization

Small-scale gas turbine combustors are essential for many applications. Burning all of the fuel in a small combustor introduces significant challenges. The fuel must atomize, vaporize, mix, and burn with the air in a shorter residence time when compared to medium- or large-scale gas turbine engines. In addition, the combustor must meet engine design criteria such as combustion efficiency, pressure loss, profile and pattern factor, and durability requirements. It is experimentally and computationally unfeasible to perform a full factorial design exploration because of the plethora of geometric design variables that can be optimized.

Conventional gas turbine combustors typically have been designed through ad hoc and heuristic guidelines based on trial-and-error experimentation loosely following one-at-a-time (OAT) local sensitivity analyses. This iterative combustor design process raises the development time and cost. Original engine manufacturers utilize computational fluid dynamics (CFD) extensively to guide some of their combustor design processes. Nevertheless, CFD has been limited to a few sample simulations that may not be representative of the entire design space. This leads to a localized optimum design that may not represent the global optimum.

Open-source software can be integrated within a framework to create design points and provide inputs to a cad software that dynamically creates combustor geometries. The geometry is passed to the meshing software to independently create meshes. The meshes are sent to the solver to perform CFD simulations. The results are postprocessed by another software application to compute output metrics.

Previously, it has been shown in Refs. [1–7] that a combustor can be improved from its baseline using numerical techniques. Our experience has suggested that either it is more appropriate to optimize the combustor in the engine or the test section geometry in which the combustor will be tested [7]. In this work, we optimize the combustor at the maximum engine power condition in a virtual test section that represents the actual test section. Because this approach does not consider engine component-to-component interactions, we refer to it as component-level design optimization. In contrast, if the combustor were optimized in the engine, which is beyond the scope of the current work, it would be considered a system-level design optimization.

Furthermore, our previous investigation on design optimization [7] suggested the following. The baseline design remained in the Pareto Frontier¹ throughout the optimization calculations. The optimization calculations showed improvement from an initial design point population to later iteration design points. The Euclidean distance from design points to the utopic point was used in the decision support process to select a “best” and “worst” design. Most importantly, the small Euclidean distances agree with the nondominated² Pareto Frontier design points. This observation is meaningful because it suggests that the current multiple objectives could be clustered into a cost function (i.e., the Euclidean distance to the utopic point). The difference between multiple-objective and single-objective optimization is that a single optimum design is obtained at the end of the procedure and no further decision support process is needed to down-select the better designs.

The main purpose of this work is to develop the open-source autonomous CFD metamodeling environment (OpenACME) for small-scale combustor design optimization. The specific objectives are: (1) to develop an automated numerical tool (based on open source codes) for combustor design metamodeling and optimization in order to reduce combustor design time and costs; (2) to utilize single-objective optimization (SOO) for the reacting flow cfd based on the Euclidean distance; and (3) to analyze the combustor reacting flow cfd utilizing nonparametric rank regression (NPRR) and global sensitivity analyses (GSA).

This section describes details related to open-source autonomous cfd metamodeling environment (OpenACME), combustor parametric geometry, computational domain, governing equations, combustion model, material properties, design of experiments (DoE), sensitivity analyses, and optimization.

Figure 1 illustrates the customized autonomous CFD metamodeling environment (OpenACME) for combustor optimization developed by coupling several object oriented programing (OOP) open-source codes consisting of freecad [8], cfmesh [9], openfoam [10,11], paraview [12,13], openmdao [14], and openmeta [15–18]. freecad is a Python-based parameterized cad software with many workbenches and provides user access to the mesher and solver GUIs. The PartDesign workbench enables parameterized editing of solid models that are sequentially added to form the desired geometry. cfmesh is a dedicated openfoam application for efficiently generating tetrahedral, polyhedral, and cutcell meshes. It is approximately 5 times faster than the standard openfoam SnappyhexMesh. openfoam is a C++ CFD toolkit. A custom openfoam application was developed that couples our flamelet/progress variable (FPV) model [19], conjugate heat transfer, and Lagrangian particle tracking. ParaView is a Python-based postprocessing code. oraclevirtualbox and docker are supporting software used to run both cfmesh and openfoam on Windows multicore workstations. openmdao is a framework for multidisciplinary analysis and optimization and comprises design of experiments, response surfaces (RS), uncertainty quantification (UQ), sensitivity analyses (SA), and optimization (O) algorithms. openmdao is directly coupled with openmeta. The latter is an overarching software for designing complex systems by coupling multiple applications. openmeta also provides the visualization capabilities for the metasimulation results.

The geometry of the outboard cavity combustor is illustrated in Fig. 2. The combustor contains several air jet groups to induce recirculation zones in the cavity, enhance fuel/air mixing, increase flame stability, and dilute high temperature zones. There are 48 cavity forward driver jets, 48 cavity aft driver jets, 8 outer liner jets, and two rows of 8 staggered inner liner jets. There are 8 fuel injector sites. A single-component surrogate species (n-C₁₂H₂₆) is used as fuel [20], which in turn is injected as a hollow liquid spray cone. Because of the periodicity of the combustor, only a 45 deg sector is modeled. The inlet conditions for all air jets are prescribed to total (stagnation) pressure of 396 kPa, total temperature of 497 K, turbulence intensity of 5%, and ratio of turbulent-to-laminar viscosity of 10. The outlet (fuel and air) mass flowrate (⁠$m˙exit$⁠) is 0.653 kg/s for the full annular combustor. The global equivalence ratio (ϕ) is 0.344. Note that mass flow rates through the various liner holes are automatically captured in the CFD solution and are not specified as boundary conditions. The combustor liner is made of Inconel-600. Constant solid density, thermal conductivity, and specific heat capacity are utilized for the liner material. The combustor liner thickness is fixed at 0.5 mm.

The fifteen geometric design variables are illustrated in Fig. 2 and described in Table 1. These design variables are chosen based on a combination of engineering judgment and previous experimental studies of trapped vortex combustors [1–6,21–25]. There are thirteen air inlet design variables including jet diameters, jet locations, chute aspect ratio and inclination angle; one parameter for the cavity axial length; and one parameter for fuel injector location. There is only one inlet mass flowrate and one pressure outlet boundary condition as in Ref. [7]. The lower and upper bounds for the design variables are spread as far as possible, constrained by manufacturing restrictions. The bounds ensure creation of realizable geometries and maximize the design space.

The computational domain is an assembly of two individual computational regions. The fluid region includes upstream and downstream plenums [7], the (engine shroud) enclosure, and the combustor. The solid region involves the combustor liner (cf. Fig. 2). The fluid and solid meshes are updated for each design point. The computational domains are meshed independently from each other. Cartesian meshing with boundary layer inflation is utilized for the fluid region, whereas tetrahedral meshing is used for the solid region. The latter is approximately two million cells with greater cell concentration near the combustor fluid region, whereas the former is approximately one million cells. Therefore, the two meshes are aggregated in the solver and are nonconformal with each other. Same number of processes are utilized for meshing the fluid and solid domains. All fluid cells and solid cells adjacent to each other are clustered into the same compute node determined automatically by the application load-balancing numerics. The maximum cell size in the fluid region is 0.5 mm. Three fluid boundary layer cells are utilized. The cell size in the solid region is constant and set to 0.5 mm.

A customized openfoam application was developed to simulate incompressible steady multiphase flow using the k–ω SST Reynolds averaged Navier-Stokes (RANS) model [26] with a “frozen” nonadiabatic diffusion flamelet/progress variable (FFPV) combustion model. This “frozen” version is discussed in the subsequent section. Turbulence-chemistry interaction is tabulated a priori in a tri-dimensional probability density function (PDF) table as a function of the low-dimensional variables (i.e., mixture fraction, mixture fraction variance, and progress variable). Beta and Dirac Delta presumed PDFs are used for the mixture fraction and progress variable, respectively. The products of these marginal PDFs are employed for the joint probability of low-dimensional manifold variables. This is used to convolute the thermochemical properties and Favre-averaged thermochemical variables such as molecular transport properties, species mass fractions, and progress variable source term. Multiregion conjugate heat transfer is utilized. The fluid governing equations are the momentum, pressure, turbulence, enthalpy, and low-dimensional manifold variable equations, whereas the solid governing equation is enthalpy. The fluid and solid enthalpy equations are coupled at the fluid-solid interfaces. Additional enthalpy equation iterations are utilized to match the temperature at the fluid-solid interface. The discrete phase model allows for the liquid to exchange mass, momentum, and energy with the gaseous phase. Second-order accurate schemes are used for the spatial discretization.

Traditional adiabatic flamelet models provide density, temperature, molecular transport properties, progress variable source term, and passive species mass fractions to the flow field governing equations. The same PDF table for adiabatic FPV can be turned into a “frozen” nonadiabatic FPV model by not passing neither temperature nor density, but still passing molecular transport properties, active species mass fractions, and progress variable source term from the PDF table to the CFD flow field. It is referred as “frozen” because the chemistry is “frozen” or independent of heat transfer. This bootstrapping technique requires the solution of the enthalpy equation like any other nonadiabatic model, but it avoids the complication of creating a pseudo-tetra-dimensional PDF table with absolute enthalpy as a fourth axis [28,29]. The change in species composition as a function of the low-manifold variables leads to temperature rise from the enthalpy equation. The density is obtained from the incompressible ideal gas equation. The species specific heat capacity and enthalpy of formation are obtained from Joint-Army-Navy-NASA-Air Force tables [20].

The “frozen” nonadiabatic FPV (FFPV model) and the adiabatic FPV formulation are compared with previous calculation [19] and measurements of Sandia D as presented in Fig. 3 for verification and validation purposes. Our flamelet combustion model developments are referred as common format routine. Note that there is good agreement between measurements, Fluent [19], and our adiabatic and “frozen” FPV models. Sandia D numerical setup is purely adiabatic, but the results in this figure shows that at this condition the FFPV model agrees with the FPV model, which confirms that our FFPV formulation is correct. Note that this FFPV model or any other “frozen” model is valid far away from limit-combustion phenomenon (i.e., ignition and blowout) such as the combustor limit-cycle operation.

Optimized Latin hypercube sampling (LHS) based on Morris–Mitchell sampling criteria [30] is utilized for the DoE. The genetic algorithm of LHS is allowed to evolve twice in order to find better coverage of the multidimensional design space. A pseudo-random seed is used to initiate each DoE. Four DoEs are generated with 75 design points, totaling 300 design point CFD calculations. Each CFD calculation is run on an eight-core virtual computing machine and it takes approximately six wall-clock hours for finalizing a simulation. There is also overhead time due to meshing and mesh partitioning.

Spearman's ρ and Kendall's τ order-rank correlation are utilized for the sensitivity analyses between design variables and output metrics. Both methods are computationally inexpensive and capture monotonicity and nonlinearity without an underlying model assumption. Kendall's τ order-rank correlation outperforms the Spearman's ρ correlation when the sample is small or contains outliers. The reason to use both NPRR models is to gain confidence in the model sensitivity between design variable and output metrics. Both Spearman's ρ and Kendalls' τ indicate strong positive correlation when the sensitivity coefficient is 1, strong negative correlation when the coefficient is −1, and no sensitivity when the coefficient is zero.

Constraint optimization by linear approximation (COBYLA) [32,33] is utilized because it is a gradient-free local optimizer for solving nonsmooth nonlinear programing problems. COBYLA is based on trust-region methods that utilizes a linear approximation of the objective function and constraints to determine successive trust region steps toward the minimum. Although COBYLA is a bound constraint optimizer, the algorithm could violate the constraint at initial stages. This would be detrimental for the subsequent calculations because it would initially lead to nonrealizable combustor designs. Therefore, the initial design variable variations are inspected to ensure that the constraints are not violated.

The ideal total pressure losses for the reacting flow design space is approximated as $TPLideal,space ≈min(TPL)reacting,space$⁠.

The station numbering is based on the society of automotive engineers (SAE) aerospace recommended practice (ARP)755 [35].

The baseline CFD results, the DoE for the combustor reacting flow CFD simulations, the NPRR, the GSA, and combustor optimizations are subsequently presented.

The flow field for the baseline design point (DP0) is shown in terms of velocity vector line integral convolution (LIC), equivalence ratio, and temperature in Figs. 4–6, respectively. The combustor cavity exhibits small circumferential recirculation zones that are only observed at the 0 deg and 11.25 deg planes. In addition, the streamlines indicate that there are several recirculation zones in the fluid region between the engine enclosure and combustor liner, including upstream of the combustor forward wall, immediately downstream the cavity aft wall, and downstream of the main inner and outer dilution jets. Similar to our previous work [7], the forward driver jet is tilted radially inward when it enters the combustor. Likewise, the chute exhibits a larger outer recirculation zone that diminishes its discharge coefficient.

The fuel is injected axially at the centerplane and at a radial location slightly inward from the forward driver jet location. The chute prevents rapid mixing of the fuel with the air as indicated by the saturated red color in Fig. 5. However, the fuel mixes with air in the circumferential direction as indicated by the “greenish” and “yellowish” colors of the 11.25 deg and 22.5 deg planes of Fig. 5.

There is a direct correlation between equivalence ratio and temperature as illustrated in Fig. 6. High temperatures are found in between the fuel injection sites (i.e., 22.5 deg plane) as well as downstream the chutes (cf., 0 deg plane). In addition, there is heat transfer from the combustor to the enclosure zone. Heat is transferred to the inner downstream regions of the enclosure zone as well as immediately downstream the cavity wall. The latter appears to slightly preheat the aft driver jet. In addition, the flame attaches to the chute outer wall due to the presence of a recirculation region. This recirculation region is also evident in Ref. [7].

The combustor liner temperature is presented in Fig. 7. The figure depicts the temperature on the inner surface of the liner. For comparison, the combustor liner temperature is computed with either conjugate heat transfer and nonadiabatic FPV combustion model or adiabatic walls and adiabatic FPV combustion model. The purpose for this comparison is to verify that the conjugate heat transfer model would reduce the liner temperature as expected by conducting heat along the liner and from the liner to the incoming air. The adiabatic formulation shows larger wall liner temperatures than the nonadiabatic formulation as expected. For the nonadiabatic conjugate heat transfer model, the temperature of the liner does exceed the 1300 K threshold. Some hot regions are observed on the inner combustor liner immediately upstream of the forward dilution holes.

The prediction of the performance of the baseline combustor (DP0) has declined from previous analysis in Ref. [7]. The exit average temperature is 1270 K corresponding to a combustion efficiency (⁠$η$⁠) of 90.4%. The maximum exit temperature is 1460 K corresponding to a pattern factor (⁠$PF$⁠) of 0.25. The critical liner area factor (⁠$Acritical$⁠) is 28.5%. The $TPL$ from station 3.1 to station 4 are 8.1%. The most concerning output metric is $Acritical$⁠, which has increased by a factor of ten from our previous study [7]. We found that $Acritical$ is sensitive to the molecular transport properties. In Ref. [7], we utilized constant molecular transport properties whereas in the current work, we use temperature-dependent properties as explained in Sec. 2.5. Heat transfer between the fluid and solid occurs at the molecular level. With temperature-dependent molecular thermal conductivity, the thermal resistance on the “hot” side of the liner is less than on the “cold” side of the liner. Thus, heat is transferred easily to the liner, but at a slower rate to the “cold” side. This effect raises the liner temperature. With constant thermal conductivity, the thermal resistance on both “hot” and “cold” side of the liner is the same. The rate of heat transfer to and from the liner is the same and, thus, the liner temperature is lower. This hypothesis was tested with an additional CFD calculation and $Acritical$ dropped by a factor of ten when constant molecular transport properties were utilized.

Figure 8 is the result of 300 CFD for the 4 × 75 LHS DoEs. This figure also presents the Pareto Frontier. The Euclidean distance varies between 0.25 and 2.0, with the former representing the best designs and the latter representing the worst designs. The Pareto frontier based on the maximization of the combustion efficiency and minimization of the pattern factor, critical liner area, and total pressure loss is indicated.

COBYLA optimization is performed on the most likely design points (from Fig. 8) that could lead to optimum designs. For brevity, only one local design optimization is illustrated in Fig. 9. It indicates the Euclidean distance ratio as a function of iteration. The convergence is nonmonotonic because contrary to analytical expressions, CFD calculation may fail, which would produce large Euclidean distance by default (e.g., iteration three). In other words, with analytical expressions the function exist in the entirety of the design space, but with CFD the results do not exists in the entirety of the design space for multiple reasons (e.g., the geometry or meshing procedure may fail, or the solver may diverge or crash). However, the results in Fig. 9 indicate that the Euclidean distance can be halved with just a few iterations by utilizing COBYLA.

The result of this COBYLA-optimum design in terms of velocity vector LIC, equivalence ratio, and temperature is presented in Figs. 10–12. The combustor cavity exhibits two uninter-rupted toroidal circumferential recirculation zones adjacent to the cavity aft wall. The streamlines indicate that there are several recirculation zones in the fluid region between the engine enclosure and combustor liner. Contrary to the baseline, there is one recirculation zone upstream of the combustor forward wall that restricts radially inward flow toward the chutes. Contrary to the baseline, there is a single large recirculation zone downstream the cavity aft wall. There is a smaller flow-velocity-recirculation zone downstream of the main inner dilution jets. There are differences in enclosure flow fields between the baseline and optimum design. This provides further evidence that a combustor must be optimized within the enclosure or considering other engine components. Similarly, the forward driver jet is tilted radially inward when it enters the combustor. Likewise, the chute exhibits a larger outer recirculation zone that diminishes its discharge coefficient.

The fuel is injected axially at the centerplane and at a radial location slightly outward from the forward driver jet location. As indicated by the saturated red color in Fig. 11, the chute prevents rapid mixing of the fuel with the air. The fuel, however, mixes with air in the circumferential direction as indicated by the “greenish” and “yellowish” colors of the 11.25 deg and 22.5 deg planes of Fig. 11.

Figure 12 shows that high temperatures are found in between at 11.25 deg plane as well as downstream the chutes (cf., 0 deg, 11.25 deg, and 22.5 deg plane). Heat transfer from the combustor to the enclosure zone still occurs as with the baseline. Heat is transferred to the inner downstream regions of the enclosure zone as well as immediately downstream the cavity wall. The aft driver jet is then slightly preheated. In addition, the flame attaches to the chute outer wall due to the presence of a recirculation region. This recirculation region is also evident in the baseline design. The combustor liner temperature has lowered substantially by comparing COBYLA-optimum design in Fig. 13 with that of the baseline (cf. left of Fig. 7).

The exit average temperature is 1350 K corresponding to a combustion efficiency (⁠$η$⁠) of 100.0%. The maximum exit temperature is 1640 K corresponding to a $PF$ of 0.34. The critical liner area factor (⁠$Acritical$⁠) is 27.0%. The $TPL$ from station 3.1 to station 4 are 13%. In summary, the COBYLA-optimum and baseline design are nondominated designs in terms of Pareto Optimality since the former outperforms the latter in terms of $η$ and $Acritical$⁠, whereas the latter outperforms the former in terms of $PF$ and $TPL$⁠. Moreover, the COBYLA-optimum is only marginally better than the baseline in terms of $Acritical$ (27% versus 28.5%). Nonetheless, the combustor liner temperature of the optimum is better than that of the baseline (cf. Figs. 7 and 13). This indicates there may be better metrics for liner durability than $Acritical$⁠. Perhaps, future studies should consider other metrics such as overall liner heat transfer.

Nonparametric rank regression and GSA are performed to understand the sensitivity of the output parameters to the input parameters. Figures 14 and 15 present the Spearman's ρ and DMIM results, respectively. Kendall's τ regression is not presented here because the results are similar to Spearman's ρ. The results indicate that according to Spearman's correlation, the magnitude of the coefficients are too small or the p-values are too large to suggest a strong correlation between input parameters and objective functions. Similarly, DMIM indicates that for the output design space there are no conspicuous design variables controlling $η$⁠, $PF$⁠, $Acritical$⁠, and $TPL$⁠. Therefore, compound design variables control these objective functions.

An open-source autonomous CFD metamodeling environment (OpenACME) is developed by integrating several OOP open-source codes. Currently, OpenACME can perform automated LHS DoE and local single-objective constrained optimization by linear interpolation (COBYLA).

OpenACME integrated a customized openfoam solver that can account for steady, incompressible, three-dimensional geometry using a multiphase k–ω SST RANS and “frozen” nonadiabatic flamelet progress variable (FFPV) combustion model. Conjugate heat transfer through the combustor liner is also considered.

The Euclidean distance (single-objective criterion) between a design point and the utopic point is based on the multi-objective criteria: combustion efficiency (⁠$η$⁠) maximization, and $PF$⁠, critical liner area factor (⁠$Acritical$⁠), and $TPL$ minimization. It was shown that $Acritical$ may not be the best criterion for liner durability and perhaps other metrics should be considered in future studies.

Nonparametric rank regression, GSA, and SOO optimization strategies are evaluated. SOO conducts LHS for reacting flow CFD for subsequent local constraint optimization by linear interpolation (COBYLA). The local optimization successfully improves the initial design condition. SOO is useful for guiding the design and development of future gas turbine combustors.

Nonparametric rank regression and GSA indicate that for the combustor reacting flow output design space there are no conspicuous design variables controlling $η$⁠, $PF$⁠, $Acritical$⁠,and$TPL$⁠. Therefore, compound design variables control these output metrics and the output design space is inherently nonsmooth and nonlinear. In summary, OpenACME proves to be a viable tool for combustor design metamodeling and optimization studies.

Special thanks to Scott Stouffer from UDRI for providing the resources and Adam Comer from ISSI for insightful discussions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the United States Air Force. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon.

• U.S. Air Force Research Laboratory (Agreement No. FA8650-15-D-2505; Funder ID: 10.13039/100006602).

$Acritical$ =

critical liner area factor

$dp$ =

design point

$d(dp,up)$ =

Euclidean distance

$n$ =

total number of pairs

$nc$ =

number of concordant pairs

$nd$ =

number of discordant pairs

$PF$ =

pattern factor

$T3.1$ =

combustor inlet temperature (K)

$T4,avg$ =

average exit temperature (K)

$T4,max$ =

maximum exit temperature (K)

$TPL$ =

total pressure losses (%)

$rgX$ =

rank-grade of variable X

$up$ =

utopic point

$δi$ =

design variable delta index

$η$ =

combustion efficiency (%)

$ρrgX,rgY$ =

Spearman's order-rank correlation coefficient of rank-grade variable X and Y

$σgX$ =

standard deviation of rank-grade variable X