A Simplified Chemical Reactor Network Approach for Aeroengine Combustion Chamber Modeling and Preliminary Design

Abstract

In a time when low emission solutions and technologies are of utmost importance regarding the sustainability of the aviation sector, this publication introduces a reduced-order physics-based model for combustion chambers of aeroengines, which is capable of reliably producing accurate pollutant emission and combustion efficiency estimations. The burner is subdivided into three volumes, with each represented by a single perfectly stirred reactor, thereby resulting in a simplified three-element serial chemical reactor network configuration, reducing complexity, and promoting the generality and ease of use of the model, without requiring the proprietary engine information needed by other such models. A tuning method is proposed to circumvent the limitations of its simplified configuration and the lack of detailed geometric data for combustors in literature. In contrast to most similar frameworks, this also provides the model with the ability to simultaneously predict the combustion efficiency and all pollutant emissions of interest ($N{O}_x$, $CO$ and unburnt hydrocarbons) more effectively by means of implementing a detailed chemical kinetics model. Validation against three correlation methods and actual aeroengine configurations demonstrates accurate performance and emission trend predictions. Integrated within two distinct combustion chamber low-emission preliminary design processes, the proposed model evaluates each new design, thereby displaying the ability to be employed in terms of optimizing a combustor’s overall performance given its sensitivity to geometric changes. Overall, the proposed model proves its worth as a reliable and valuable tool for use towards a greener future in aviation.

1. Introduction

The aviation sector as a whole contributes today between 2–5% of the total greenhouse gasses (GHGs) produced by human activity [1,2,3,4]. The plan laid out to implement carbon neutrality by 2050 [5,6,7] includes the design of more fuel-efficient propulsion systems that produce lower GHG emissions [8] while burning either conventional fuel or new greener fuels, also known as sustainable aviation fuels (SAFs), such as drop-in biokerosene or non-drop-in carbon-free fuels such as hydrogen [5,9]. Therefore, the development of robust and reliable models of the combustion process, incorporating a better understanding of pollutant formation mechanisms, has a dual advantage. The implementation of SAFs in aircraft engines requires a robust method for predicting their effects on the engine’s GHG emissions and performance. The incorporation of such models in multidisciplinary preliminary design frameworks may substantially decrease both risk and costs while expediting the implementation of low-emissions technologies.

The most accurate yet computationally expensive solution for the in-depth simulation of the combustion process, in aeroengines, is through the use of computational fluid dynamics (CFD) packages. The complexity of combustion chamber geometry (rarely depicted in detail), the tedious meshing and configuration process, and the costly execution time are significant drawbacks, thus making 2D or 3D CFD analysis more suitable for later design stages, thus taking place after the preliminary stage [10,11]. Yet, modern advances in high-performance computing and CFD efficiency show that detailed chemical analysis is feasible, as depicted in the work of Oliveira et al. [12], in which the CFM56-3 engine’s burner $C{O}_2$, $N{O}_x$, $CO$ and UHC emissions were predicted for different working conditions and different fuels, including drop-in biofuels. In the aforementioned work, and in accordance with similar experimental research [9,13,14,15] it has been confirmed that the usage of different fuels causes changes in $CO$ and UHC emissions, which in turn alter the combustion efficiency of the burner.

Another approach for the estimation of gas turbine burners’ performance parameters are correlation methods, which can be empirical or semiempirical. Due to their simplicity and their wide range of comprehensive implementation, correlation emission models are widely adopted in gas turbine cycle analysis and preliminary design frameworks [16,17,18]. Yet, given that these methods are data-driven from specific up-to-the-present-day engines, they cannot forecast the effect of design change, while extrapolations outside of their respective data range are of questionable accuracy and are better to be avoided. Such methods exist for both emission and efficiency predictions, yet they are mostly independent of one another, even though (contrary to $N{O}_x$) $CO$ and UHC emissions are tightly linked with the burner’s efficiency.

Recently, physics-based reduced-order models are gaining a lot of attention, especially for their implementation in the conceptual design stage of the burner and the engine [8,19,20]. These models follow a control volume approach with the different zones of the burner divided into finite volumes, in which relatively similar flow phenomena and chemical species concentrations occur. In each of these volumes, the governing physical and chemical equations are solved iteratively in an integral form [10,21]. They require higher computational resources compared to correlation-based models, and are less dependent upon the combustor configuration while yielding greater consistency in their outcomes. Compared to respective CFD-based methods, they mostly exhibit lower accuracy, but considerably less computational cost. The most frequently adopted method of this type is the chemical reactor network (CRN) approach, where each control volume is modeled as one or more chemical reactors in which reduced order transport and finite chemistry equations with fuel-specific kinetic mechanisms simulate the process for a predetermined fuel-to-air ratio (FAR). The different reactor types in CRN models are the plug flow reactor (PFR), the perfectly stirred reactor (PSR), and the partially stirred reactor (PaSR) [10,22]. The PFR assumes a 1D axial flow with no radial diffusion and is considered suitable for simulating regions of the burner where turbulent mixing is low. The PSR and PaSR consider a 0D flow in which the species, prior or subsequent to the main reaction, are uniformly or partially mixed, respectively. Yet, none of the aforementioned reactors can simulate by themselves the axial change of chemical composition of the working gas along the length of the burner, which is caused by the combustion, as well as the dilution and cooling air streams added downstream of the main flame front. This problem is circumvented through the implementation PSRs, meaning networks of PSRs or PaSRs, each with a predefined mean FAR and mixing quality, so as to better replicate the mass transport and temperature distribution within each zone, thereby creating a discretized pseudo-1D flow representation.

State-of-the-art CRN models consist of PSRs, especially when it comes to models with potential to be incorporated within design loops performed for combustion chamber preliminary sizing [8,23,24]. In particular, concerning the simulation of the most physically complex region of the combustor, the primary zone (PZ) within which the flame front occurs, the authors favor the use of PaSR due to the need to model the significant fluctuations occurring in the mixing quality within this turbulent zone. In the literature, these PaSRs are often constructed as a series of parallel [8,10,24] or recirculating [23] PSRs, with each reactor network trying to replicate the corresponding flow field’s stream lines. Each reactor, forming the aforementioned PSRs for the PZ, is simulating a fraction of the overall combustion, with a FAR following a discretized Gaussian distribution, thus recreating the effect of the imperfect mixing of fuel–air mixture [24]. These variations in the FAR within the PSRs can capture fluctuation in emission production, especially $N{O}_x$, given that its formation is primarily temperature dependent, by reproducing post-combustion temperature perturbations occurring between different subregions of the PZ. The partial mixing in the PZ is dependent upon many factors such as burner geometry, operating conditions, fuel injection, etc. and, therefore, design-specific semiempirical correlations are used [24]. Yet, the employment of partial mixing models for the PZ in CRNs causes increases in computational cost and the diminishment of generality while not always yielding an increased accuracy of the predictions [8]. The use of PFRs for modeling the whole burner is less frequent but still implemented [25] in cases with weak turbulent mixing. For the secondary (SZ) and the dilution zones (DZs), up to two serial PSRs or one PFR are used [8,23,24], thereby considering that an almost perfectly mixed gas flows through them.

When it comes to the reaction mechanisms applied in CRNs to simulate the chemical subreactions of the combustion taking place within the burner, two options can be distinguished. On the one hand, simplified, also known as surrogate reaction mechanisms with less potential species and possible reactions are implemented to decrease the range of estimated pollutant emissions to also reduce its cost [10]. On the other hand, detailed mechanisms, taking into account all possibly occuring species and subreactions occurring within the simulated global reaction, have also been employed in recent work [8,23,24,25].

Furthermore, when evaluating the monitored output of CRN models found in corresponding research, it is deduced that it mainly consists of emission estimations of the gas turbine burners, considering mostly $N{O}_x$ [8,23,24] and sometimes also $CO$ [10], always accomplishing the expected trends. Very few research implements CRN combustor models for the estimation of combustor efficiency with respect to working conditions and fuel [24], while even rarer is the evaluation of UHC emissions with such models.

In the current work, a simplified CRN combustor model for aeroengines is developed, validated, and applied within a respective design process. The proposed CRN consists of a generalist, yet design-sensitive physics-based reduced-order model capable of simultaneously predicting the trends of an aeroengine combustion efficiency and its produced emissions. To verify the prediction capabilities of the model under consideration, the results produced from it are compared with corresponding ones from correlations methods, aero-thermodynamic models, and emission test data from the International Civil Aviation Organization (ICAO) databank for three engines utilized in civil aviation: the CFM56-7B2, the CFM LEAP-1A2, and the RR TRENT 772. Finally, in order to prove the competence of the model in preliminary design, a combustion chamber design application is also employed, thereby incorporating the model in question for the purposes of emission and performance evaluation.

2. Methodology

In this section, the assumptions made, the numerical methods used, the proposed model concept, the virtual structure, and the configuration are discussed.

2.1. Overview

The model proposed here was coded in the source C++ environment of the CANTERA open-source software package in order to produce a tool, that would be as least computationally demanding as possible, within the capabilities of the current package and easy to incorporate, as an external function form, within object-oriented gas turbine thermodynamic analysis and design frameworks, such as PROOSIS [26]. For the purposes of this study, a detailed kerosene–air reaction mechanism was employed so as to produce predictions for all emissions of interest in aviation: $N{O}_x$, $CO$ and UHC. Yet, the plug-in configuration of the mechanism in CANTERA makes the model suitable for possible simulation of alternative fuel combustion should such mechanisms for SAFs become available. In contrast to common practice found in literature, the partial mixing in the PZ was neglected, and a simplified PSR structure was devised to require the least possible geometric and flow data by following a proposed three-zone division approach (primary, secondary, and dilution), thereby making it as generalized as possible for conventional isobaric combustors while still yielding reliable emission and efficiency tendencies with respect to power setting variations. Lastly, the output state of the combustion products was utilized to make a fuel-dependent estimation of the burner’s combustion efficiency, thus opening the possibility to address the need to investigate the burner’s efficiency variations when using SAFs.

2.2. Modeling Approach

The CRN method for modeling gas turbine combustion chambers assumes that the overall volume of the burner, in which the combustion of the fuel and the dilution of its products with cooler air streams take place, can be divided into a number of subvolumes, with each of them having homogeneous flow characteristics and containing perfectly mixed chemical species [24]. When considering a conventional gas turbine burner, it is common practice for most authors [10,11,18,23,24] to divide it into three regions of constant volume, with reference to the spatial distribution of the air streams inlets, also known as dilution holes, across the liner. These three regions, being referred to as the PZ, the SZ and the DZ, are each assumed to have different yet constant mean ratios of total (initially unburnt) fuel flow to air flow introduced in each zone according to their respective positions with respect to each dilution hole row, as visualized in Figure 1. These mean ratios in each zone are referred to from now on as $FA{R}_{PZ},FA{R}_{SZ},$ and $FA{R}_{DZ}$, respectively, for reasons of simplicity. In practice, while the overall volume of the burner remains unchanged, the volume of each zone may slightly change due to differentiations in gas density, axial velocity, and flame temperature under partial loads. Yet, the inclusion of semiempirical correlations to model the reactor volume change with respect to the load or the consideration of the zone length ratio as an extra design variable, which can independently fluctuate for each examined working condition, would add to the complexity of the model, which is not desired.

In the present work, the physical processes in each zone were simulated through the use of a single PSR, with each interconnected with its respective inlet and outlet to form a serial PSRs configuration to model the combustor as a whole. In order to safely assume the implementation of a PSR to model a volume, it is necessary for the characteristic flow time, also known as the mixing time ${\tau }_m$, to be significantly smaller than the characteristic reaction time ${\tau }_r$, which in turn requires high turbulent mixing [23,24]. This criterion implies a very low Damkohler number, which is defined as the ratio of the mixing time to the reaction time ($Da={\tau }_m/{\tau }_r$), and is met in aeroengine combustors, given that mixing time is, generally, sufficiently low to allow for near-perfect completion of the combustion. The proposed partitioning is as simplistic and generalistic as possible so as to require the lowest amount possible of burner-specific gas flow data, thus not requiring the partition to be redesigned if the geometry is modified. The inclusion of parallel PSRs for the PZ or SZ modeling would require more engine-specific data or ad hoc assumptions to determine the cooling and secondary air mass flow or the corresponding equivalence ratio found in each PSR of the respective zone, which comes in conflict with the simplistic modeling strategy. This is the main reason for which the partitioning is referred to as semi-1D, given that it does not opt to precisely reflect the exact flow processes occurring into an actual combustor but to capture the correlation between the basic burner geometry and the combustion performance. The constant volume of each zone and the constant pressure Dirichlet boundary conditions imposed at the inlet and outlet of the burner have an impact on the gas density and temperature before and after the fuel ignition in the PZ. Additionally, the simplified air flow partitioning (which correlates to the hole geometry) defines the equivalence ratio in each zone. which has a direct impact on the post-deflagration temperature. We argue that, for a high-generality and low-order model such as the one proposed, the aforementioned zone partitioning sufficiently captures the influence of the zone volume and the basic air streams on the post-deflagration temperature, which in turn has a direct impact on the combustion efficiency and pollutant formation.

During the combustion of a fuel, the hydrocarbon chains break due to localized high temperatures (cracking), thus forming smaller organic compounds, which in turn react with the air or with each other to form products of combustion. Simultaneously, chemical species in the air react with each other (e.g., $O_2$ with $N_2$), thereby creating a wide range of oxides with varying volatilities. The aforementioned reactions, along with several other intermediate reactions, occur for every global reaction and regulate the proportions of reactants and products until chemical equilibrium is reached in the system. All the above processes govern the chemical phenomena transpiring within a combustor and can be used to define the system of ordinary differential equations (ODEs) governing the molar concentration of each species with respect to time, which is known as the chemical reaction mechanism or combustion mechanism [10]. These mechanisms are implemented in CRN models to simulate the chemical kinetics of various reactants.

Each mechanism is specifically designed and extensively tested to simulate a specific reaction, and it is even possible to require two distinct mechanisms for the same reaction due to the fact that different chemical processes occur with respect to different conditions (i.e., rich or lean mixture) [27]. Mechanisms that include the widest known ranges of species and reactions transpiring within the referred reaction are known as detailed mechanisms, with their cost of use scaling with the number of both species and reactions. Surrogate mechanisms simplify the ODEs by reducing the species and reactions taken into account, and thus their cost by focusing on the inclusion of formation reactions regulating the pollutants to be monitored (i.e., $N{O}_x$). In the literature, both detailed [25] and surrogate mechanisms [8,10,23,24] have been utilized in CRN model applications by making trade-offs between emission accuracy and speed of execution. In this paper, the detailed combustion mechanism of J. Luche [27] for jet fuel was employed, within the CANTERA framework, as a plug-in feature.

2.3. CRN Formulation and Structure

CANTERA version 2.6.0 was utilized to solve the chemical kinetics, thermodynamics, and transport equations in the PSRs. The components and the structure of the model are described below, together with the postprocess and the tuning method.

2.3.1. CRN Virtual Components

In our model, the implemented components, also known as the object classes of the CANTERA object-oriented modeling library, are Reservoirs, MassFlowControllers, Valves, and Reactors.

The Reservoir is a space of infinite volume and constant pressure in and out of which mass can be transported. This serves to keep unchanged the overall state and composition of the Reservoir’s gas in order for the component in question to keep a chemically consistent output when needed.

The MassFlowController is an object that determines the mass flow from an upstream component to a downstream component. The mass flow rate ${\dot{m}}_c$ can be time-dependent or constant.

The Valve is an object that links two components in order to regulate the pressure drop between them, thereby determining the mass flow rate from the upstream to the downstream component.

The Reactor is a PSR object; namely, it is a finite volume with homogeneously mixed chemical species in which reactions between those species can occur.

More about the governing equations and the theoretical background of each of the aforementioned components, other than the Reservoir, can be found in the relevant literature [22,28].

2.3.2. Layout

The overall layout of the proposed CRN model is shown in Figure 2. The fuel and air are introduced into the different zones through the aforementioned MassFlowControllers. Three air ratio parameters are defined as the mass fractions of the air added into each zone, which are assumed constant and independent with respect to the burner working conditions in accordance with the literature [10,23], given that the diameters of the holes remain unchanged. The formulation of the air ratios reads as follows:

$$A{R}_i=\frac{{\dot{m}}_i}{{\dot{m}}_a}=const.\&\sum _{i}A{R}_i=1,i=\left\{PZ,SZ,DZ\right\}$$

(1)

The value of the air flux and, therefore, the value of the corresponding FAR and equivalent ratio ($\varphi$) in each zone are determined by the $AR$ parameters, as well as the known overall air ${\dot{m}}_a$ and fuel flow ${\dot{m}}_f$ of the burner. Hence, the mass flow and FAR in each zone, coming from the air reservoir, are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{m}}_i=A{R}_i·{\dot{m}}_a\&& FA{R}_i={\dot{m}}_f/\left({\dot{m}}_a\sum _{n=PZ}^{i}A{R}_n\right)=FAR/\sum _{n=PZ}^{i}A{R}_n\hfill \\ \hfill & fori=\left\{PZ,SZ,DZ\right\}\hfill \end{array}\end{array}$$

(2)

where $FAR={\dot{m}}_f/{\dot{m}}_a$, and, thus, the equivalence ratio in each zone is defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\varphi }_{PZ}& =\frac{FAR}{FA{R}_{st}·A{R}_{PZ}},{\varphi }_{SZ}=\frac{FAR}{FA{R}_{st}·\left(A{R}_{PZ}+A{R}_{SZ}\right)},\hfill \\ \hfill {\varphi }_{DZ}& =\frac{FAR}{FA{R}_{st}·\left(A{R}_{PZ}+A{R}_{SZ}+A{R}_{DZ}\right)}=\frac{FAR}{FA{R}_{st}}=\varphi \hfill \end{array}\end{array}$$

(3)

where $FA{R}_{st}$ is the fuel-to-air ratio for which stoichiometric combustion of the given fuel occurs (no oxygen excess after a perfect combustion).

A Valve is placed after the DZ reactor in order to define the outlet pressure of the component, thus achieving the recreation of an isobaric combustor by assuming that $p_{out}^{\left(DZ\right)}=p_{exhaust}$. By including a pressure drop user-defined parameter $dPqP$, the pressure of the exhaust Reservoir is assumed as $p_{exhaust}=dPqP·p_{in}^{\left(PZ\right)}$.

The igniter configuration is one of the two ways available to initialize the combustion, as deduced from the combustion modeling CANTERA documentation [22,28]. If the igniter option is discarded, then it is possible to initiate the reaction by introducing in the PZ the reactants of the complete combustion of the air–fuel mixture within it so that the heat released from the ideal reaction in the PZ may keep the fuel burning in the other zones until a steady solution is achieved. Otherwise, the igniter introduces a time-dependent Gaussian pulse of highly volatile atomic hydrogen (H) in the PZ, at the same temperature and pressure as the fuel, so as to instantly react with the air and release enough heat to cause the fuel to react: firstly in the PZ and then in the rest. The igniter approach is preferred for the combustion modeling in this paper, given that it better emulates the combustion of the partially mixed air–fuel mixture in the PZ for a wider range of the burner’s working conditions. The hydrogen mass flux passing through MassFlowController connecting the igniter Reservoir with the PZ assumes the form

$${\dot{m}}_{ign}=m_f·g\left(t\right),g\left(t\right)=A{e}^{-{\left(\frac{t-t_0}{\tau }\right)}^2}$$

(4)

where its parameters are set: $m_f={\dot{m}}_f$ is the constant fuel mass flow, $A=$ 0.04–0.12 is the scalar defining the dimensionless height of the Gaussian mass flow distribution, $t_0$ refers to the ignition start time, and $\tau =\frac{P{W}_{1/2}}{2\sqrt{2}}$, with $P{W}_{1/2}$ being the width of the distribution at the half of its height, which is set between $0.2{\dot{m}}_f$ and $0.7{\dot{m}}_f$.

An example of the time-dependant profile of the temperature at the outlet of each zone when using the hydrogen igniter, for a burner simulation from an engine operating in cruise state (see Section 3.1), is plotted in Figure 3.

The inputs of the model are the following: $p_{t,in}$, $T_{t,in}$, ${\dot{m}}_a$, $FAR$, $A_{ref}$, L, $t_{all}$, $N_t$, $A{R}_{i}i=\{PZ,SZ,DZ\}$, the combustion mechanism file, and the fuel composition. The outputs of the model under consideration are $p_i$, $T_i$, ${\varphi }_i$, and ${\dot{m}}_i$ per zone of the combustor, as well as the $EI$ of $N{O}_x$, $CO$, UHCs, $n_b$, and $LHV$ of the fuel. All symbols are described in the Abbreviations Section.

2.3.3. Postprocess

Once the main simulation of the combustion is completed and a steady state solution is achieved, the working gas state or its species concentration at the outlet of the DZ are utilized to make estimations of its combustion efficiency ($n_b$) by implementing two distinct methods. The first method is determined by means of the following equation:

$$n_b=\frac{\Delta h_{t,out}-\Delta h_{t,in}+FAR(\Delta h_{t,out}-\Delta h_f)}{FAR·LHV}$$

(5)

where $\Delta h_{t,j}=h(T_j,FAR)-h(T_{ref}=15{}^{\circ }\mathrm{C},FAR),j=\left\{in,out\right\}$ is the enthalpy per mass unit of the working gas relative to a reference temperature, and $\Delta h_f$ is the enthalpy per mass unit of the fuel relative to the same reference temperature. The process for the calculation of $LHV$ is described in Appendix A.1.

The second method assumes that the difference between the available and the released heat can be computed by subtracting from the lower heating value ($LHV$) of the fuel, the heat that would be released from the products of imperfect combustion ($CO$ and UHCs), as calculated in the following equation:

$$n_b=\frac{LHV-Q_{CO}-Q_{UHC}}{LHV}$$

(6)

where $Q_{CO}$ and $Q_{UHC}$ define the heat that could be released by the combustion of the remaining corresponding species in the product gasses. The algorithm used for the calculation of these quantities is given in Appendix A.2.

The mass fractions $Y_k$ (in $kg$ of species $-k$ per $kg$ of total products) of the products of the combustion at the burner’s outlet provided by the CANTERA solver are employed to calculate the $EI$ of each pollutant of interest. The general formula for the computation of the $EI$ yields

$$E{I}_k=\frac{(1+FAR)}{FAR}·Y_k·{10}^3\left[\frac{grofspecies-k}{kgoffuel}\right],k=\left\{N{O}_x,CO,UHC\right\}$$

(7)

In particular, $Y_{N{O}_x}=Y_{NO}+Y_{N{O}_2}$ and $Y_{UHC}$ define the sum of all the hydrocarbon species accounted for in the current mechanism and remaining in the working gas as the products of incomplete combustion, respectively.

2.4. Model Tuning

In order to produce burner performance and emission prediction results throughout the whole flight envelope of an engine, data related to the design of the burner are needed, which are usually not provided by the manufacturer and are rarely found in the literature. The liner length (L), the reference area and height of the burner ($A_{ref}$ and $h_{ref}$, respectively), the length of each zone, etc. are geometrical quantities that can be easily assumed by reading the cross-section of an engine. The mass fraction of air introduced in each zone, quantities which are proportional to the number and diameter dilution holes along the length of the liner and the dome, as well as the swirler, also known as the inlet of the PZ, are far harder to assume, if not given, without some knowledge of the flow field around the dome and the liner. Yet, these quantities define the FAR and the equivalence ratio in each zone, which were described previously in Equations (2) and (3), respectively. Therefore, the need to define the mass fractions of the air added into each zone ($A{R}_{PZ},A{R}_{SZ},A{R}_{DZ}$) arises.

According to the literature related to isobaric combustor modeling [10,29], in the typical working conditions of the engine, also known as the cruise phase, the equivalence ratios in the burner’s zones should be bound within ${\varphi }_{PZ}=1-1.09$ for the PZ and within ${\varphi }_{SZ}=0.58-0.8$ for the SZ, while in the DZ equivalence ratio, the zone is equal to the overall FAR, as defined in Equation (3). The values of ${\varphi }_{PZ}\&{\varphi }_{SZ}$ can be modified within these respective ranges in order to produce consistently accurate emission estimations at different working conditions of the engine (e.g., take-off, idle, etc.). Once the appropriate values of ${\varphi }_{PZ}\&{\varphi }_{SZ}$ in the cruise conditions are selected, they define the air flow in each zone through parameters $A{R}_j$ by means of Equation (3), which in turn define the dilution hole area for each zone, which remains constant for all other examined working conditions.

3. Verification Examples

When evaluating models such as the one proposed, given that their implementations are opted mainly for preliminary design, the validity of the trends of their results with respect to the working conditions are prioritized over the actual accuracy of these results. Therefore, the burner emissions and efficiency estimations computed through the current CRN model were compared with corresponding reference values from experimental engine testing from the International Civil Aviation Organization (ICAO) databank [30] for the emission reference, fine-tuned correlations, and aerothermodynamic engine models for the efficiency reference.

3.1. Engine Test Cases

For the verification of GHG emission estimations according to the proposed model, emission data for three commercial civil turbofan engines were used. The selection of the engines was made with the prospect of covering a significant spectrum of technological advancements and aircraft sizes in order to validate the model’s generality. Two twin-spool high-bypass turbofans for narrow body aircraft were studied, which are an early version of the CFM56 and the CFM LEAP. A three-spool high-bypass turbofan engine, from the RR Trent 700 series, was also assessed to encompass an engine implemented in wide-body aircraft.

The burner inlet conditions used as input for the evaluated CRN model were produced through thermodynamic models of each engine in the PROOSIS [26] modeling environment, and they were calibrated to yield the expected thrust and fuel consumption for a range of working conditions. Such data, used as input for the CRN model, are given in Appendix B. These models were coupled with an empirical efficiency estimation method [29], which correlates the efficiency of a combustor with the loading parameter $\omega$ as described in the following equation. Therefore, the efficiency predictions computed through the aero-thermodynamic engine models were also compared with the corresponding CRN-generated estimations.

$$n_b=1-exp\left(ln(1-n_b^{\left(ref\right)})+PLC·ln\left|\frac{\omega }{{\omega }_{ref}}\right|\right),\omega =\frac{{\dot{m}}_a}{{\left(\frac{p_{t,in}}{p_{t,std}}\right)}^{1.8}·{10}^{0.00145(T_{t,in}-400)}}$$

(8)

$p_{t,in}$ and $T_{t,in}$ are measured in $bars$ and $Kelvins$, respectively, while according to the corresponding literature [29], the $PLC$ is assumed to be $1.6$. For this specific model, the reference values of the aforementioned inputs were assumed for the take-off working condition of the engine, where the $n_b$ is expected to be closest to 1 (e.g., $n_b^{\left(ref\right)}=0.999$).

$E{I}_{NOx}$ is also computed by implementing the correlation method of Odgers et al. [31] for the sake of comparing the CRN model results with a widely accepted correlation method. The origin of data for this model implies that it can be utilized for the examined CFM56-7B2 and RR Trent 772 engines and not for the CFM LEAP-1A2. This method is semiempirical, given that it predicts the $E{I}_{NOx}$ by use of the parameters $p_{t,in}$ (in Pa), $T_{fl}$ in (in Kelvin), and $t_{form}$ (in s), as defined in the following equation:

$$E{I}_{NOx}=29·exp\left(\frac{-21670}{T_{fl}}\right)·p_{t,in}^{0.66}·\left[1-exp\left(-250·t_{form}\right)\right]$$

(9)

where $T_{fl}$ is assumed to be equal to the PZ output temperature, $T_{PZ}$ is predicted from the corresponding simulation of the proposed model, $t_{form}$ is assumed equal to 1 ms and the $E{I}_{NOx}$ is computed in terms of gr of $N{O}_x$ per kg of fuel.

In Figure 4, Figure 5 and Figure 6, the model simulation results for the engines CFM56-7B27, LEAP-1A26, and RR Trent 772, respectively, are presented. All the results generated by the CRN model employed the igniter configuration. In accordance with the tuning method depicted in Section 2.4, the values of ${\varphi }_{PZ}$ and ${\varphi }_{SZ}$ in cruise were modified to best-fit the ICAO emission values for all the evaluated conditions: take-off (100% PS), climb (85% PS), approach (30% PS), and idle (7% PS). The selected $\varphi$ values and their corresponding $A{R}_j$, as well as all other geometrical inputs of the CRN model to produce the following predictions, are displayed in

Table A2. Burner geometry for the assumed ${\varphi }_{PZ}$ and ${\varphi }_{SZ}$.

Parameter	CFM56-7B27	CFM LEAP-1A26	RR TRENT 772
$A_{ref}$ [m${}^2$]	0.160	0.225	0.180
L [m]	0.178	0.157	0.191
${\varphi }_{PZ}$ *	1.057	1.090	1.024
${\varphi }_{SZ}$ *	0.605	0.600	0.580
$A{R}_{PZ}[\%]$	30.82	29.14	28.46
$A{R}_{SZ}[\%]$	23.02	23.80	21.78
$A{R}_{DZ}[\%]$	46.15	47.05	49.76

* in the cruise conditions.

of Appendix B. Accordingly, the ICAO trends were met in all the evaluated predictions. The $N{O}_x$ emissions and the efficiency displayed a notable absolute accuracy for all the power settings in all simulations. The $CO$ and UHC emissions estimations were of varying accuracy values depending on the simulation. For the cases in which the comparison with the correlation method was possible, the results of both methods were consistently comparable in terms of accuracy. The proposed model seems capable of reproducing correct trends with varying absolute accuracy over a wide range of different applications, while simultaneously computing all emissions of interest and the combustion efficiency.

The computational cost of an evaluation of one of the examined working conditions with the proposed model required a mean time of approximately 15.6 s on an HP Prodesk 400 G7 Microtower PC equipped with a CPU Intel(R) Core(TM) i5-10400 6-core CPU with a base speed of 2.90 GHz.

3.2. Efficiency Study

Further studying the combustion efficiency predictions of the proposed model comparison are held with two correlation methods. The first implemented correlation model is empirical and is the one previously described by Equation (8) in Section 3.1. The second method used is semiempirical and correlates the efficiency with a $\theta$ parameter that is dependent on ${\dot{m}}_a$ (in lb/s), $p_{in}$ (in psi), $A_{ref}$ and $h_{ref}$ (in $i{n}^2$ and $in$, respectively), and $T_{t,in}$ (in ${}^{o}R$), as well as ${\varphi }_{PZ}$, as formulated in the equation below:

$$\theta =\frac{p_{t,in}^{1.75}{A}_{ref}{h}_{ref}}{{\dot{m}}_a}exp\left(\frac{T_{t,in}}{b}\right){10}^{-5},b=\left\{\begin{array}{cc}\hfill 382\left[\sqrt{2}+ln({\varphi }_{PZ}/1.03)\right]& for{\varphi }_{PZ}<1.03\\ \\ \hfill 382\left[\sqrt{2}-ln({\varphi }_{PZ}/1.03)\right]& for{\varphi }_{PZ}\ge 1.03\end{array}\right.$$

(10)

In Equation (10), all the variables are measured in the aforementioned imperial units. To calculate the $n_b$ of a burner once the $\theta$ and b parameters are computed, through Equation (10), interpolations are made on an $n_b$ versus the $\theta$ map, which includes $iso-b$ lines corresponding to different ${\varphi }_{PZ}$ values. This map has been produced through the accumulation of actual engine data over the years and can be found from multiple literature sources, such as [18].

The proposed CRN model estimations of the efficiency were produced through Equations (5) and (6) for the steady state of the simulation during the postprocess. The two methods utilized for its calculation yielded results with minor deviations.

3.2.1. Parametric Analysis

To validate the accuracy of the sensitivity of the predicted efficiency with respect to variations in the working conditions, a parametric analysis was held. The reference value of the efficiency was assumed at the ICAO cruise condition of the CFM56-7B27 engine, and the relative divergence from from the engine was computed for a $\pm 10\%$ perturbation of the inlet pressure, temperature, and mass flow, respectively. This process took place for the proposed model and for both of the correlation methods discussed above, and the results of the analysis are plotted in Figure 7. When it comes to combustor efficiency estimation, the model in question seems to follow the same trends as the two implemented correlation models for the evaluated fluctuations of the inlet conditions.

3.2.2. Efficiency Map Generation

In order to verify the efficiency estimations of the CRN model, not only for small perturbations but for a wider range of operating conditions, a re-generation of the semiempirical $n_b$ versus $\theta$ map was attempted. To achieve this, we attempted to recreate an $iso-$b line by first accordingly setting the ${\varphi }_{PZ}$ from Equation (10), in the cruise input configuration for the CFM56, and then by executing the CRN model for an increasing set of inlet pressure values to attain the wanted span of $\theta$ values on the reference map of the semiempirical model. This process was repeated until four $iso-$b lines were successfully generated for $b=\left\{300,350,400,530\right\}$. The obtained results are presented in Figure 8, where the generated $iso-b$ lines display significant proximity to the reference lines, especially for $\theta \ge 2$. For $\theta <2$, the CRN results were slightly overestimated, yet this could be expected given that the semiempirical map is generalist in nature and it should not be expected to perfectly replicate the efficiency of a specific engine’s burner. The correct trend of the efficiency predictions over a wide range of input conditions has been, therefore, demonstrated.

4. Preliminary Design

Having shown that the combustor analysis method produces reasonable predictions, it was decided to incorporate it in a combustor preliminary design process. A simplified one-dimensional design algorithm was implemented in order to produce an initial geometry. That geometry was then fitted to the CRN burner model, which evaluated the design by producing a reliable estimate of the combustor performance and emissions. The geometry was then modified in order to reach optimal design, and to achieve this, two methods were employed:

• A correlation algorithm based upon Mellor et al. [32], opting to decrease $N{O}_x$ emissions and keep the $CO$ emissions constant, yielded a slightly modified geometry, which, in turn, the CRN model evaluated with respect to its landing and take-off (LTO) cycle emissions.
• An optimization loop was employed in which the CRN model was integrated so as to produce optimized LTO cycle emissions.

The benefit offered to the overall design process is that it allows dimensioning of the combustor, as well as the possibility to calculate its weight, instead of relying on empirical formulations [18,33].

4.1. Combustor Design Algorithms

The methodology utilized is the one shown in Figure 9. In order to calculate the dimensions of each component of the burner, semiempirical relations and data were used. Such relations and data were drawn from Fletcher et al. [29], Lefebvre et al. [18], Mattingly [33], and Mellor [34]. The first component tp be designed was the diffuser. The total pressure, temperature, and mass flow at the exit of the compressor were used as inputs, which came from aero-thermodynamic cycle analysis of the engine. The Mach number at the exit of the diffuser needed to be known, as well as information about the geometry of the compressor exit. Once the diffuser length and height were derived, the total pressure at its exit was also evaluated.

Subsequently, a flow partition algorithm calculated the distribution of air flow into the different zones of the combustor (PZ, SZ, and DZ). The parameters needed to calculate the flow distribution are the mass flow of the fuel, the equivalence ratios at the SZ and DZ, the cooling mechanism (effusion, transpiration, or film), and the cooling efficiency. The reference area can be calculated from the equation proposed by Lefebvre et al. [18], or the annulus area can be calculated from the equation proposed by Fletcher et al. [29]. The liner area can then be calculated via the calculation of the optimum ratio. That optimum ratio can be calculated either from Lefebvre et al. [18] or Mattingly [33].

Dimensions of the different components of the combustor were then defined, namely, the swirler, dome, and snout. For the design of the swirler, an equation and values of the design variables were obtained from Lefebvre et al. [18]. The design variables are the angle of vanes, the number of vanes, the thickness of the vanes, and the type of the vanes (curved or flat).

The final step of the initial sizing process is the calculation of the lengths of the different zones in the liner and the calculation of the number and diameter of the holes that are placed on the liner. For the calculation of the lengths, different methodologies can be implemented. All of them have in common that the length of the PZ is calculated via the equation provided from Mattingly [33]. Lastly, for the sizing of the holes, the equation of the maximum depth penetration of the jets was used. This equation was derived from Lefebvre et al. [18].

Once a first estimation of the combustor dimensions was achieved, the two aforementioned optimization methods were implemented. The first algorithm was based upon Mellor et al. [32] and is visualized in Figure 10a. It recalculates the length of the PZ and SZ in order to achieve lower $N{O}_x$ and $CO$ emissions while maintaining the length (L) and reference area of the combustor ($A_{ref})$ by adding as additional inputs the desired $E{I}_{NOx}$ at take-off and the $E{I}_{CO}$ at idle. This algorithm acts as a tuning of the initial preliminary design in order to achieve lower emissions. It must be mentioned that this algorithm defines differently the zones lengths in the liner, which will be seen in Section 4.2. The second method comprises a typical design optimization loop, shown in Figure 10b, using some or all of the geometrical inputs of the proposed CRN model as design variables. Different optimization cases are conducted for both single- and multiobjective optimization of the LTO emissions of the designed burner.

4.2. Application

In order to assess the fidelity of the produced combustor initial designs, data published in the open literature were used, and the obtained design variables were compared to those that have been published. Data for three engines (JT9D, TF41, and J79) were drawn from Mattingly [35].

By implementing the equation of Lefebvre et al. [18] for calculating the $A_{ref}$ within the initial design process described in Figure 9, the computed liner length L displayed a relative variation from the actual value, which is bound within $\pm 3\%$.

Then the initial design algorithm was applied to the input data corresponding to the CFM56 burner according to published data taken from [10]. The produced initial design is visualized in a simplified manner in Figure 11.

Finally, four optimization scenarios were held using the aforementioned geometry, either through Mellor’s process or by incorporating the CRN model in an optimization loop for different objective functions and objectives:

• Mellor: The algorithm for low emissions from Mellor et al. [32] was implemented. The aim was to decrease by 15% the $N{O}_x$ emissions while retaining the same $CO$ emissions. The chosen objective combination was the same as the one applied in Mellor’s publication [34].
• Opt1: The selected design variables were ${\varphi }_{PZ}$ and ${\varphi }_{SZ}$ at cruise, which are bound between 0.8–1.09 and 0.58–0.8, respectively, as well as the length ratio of the PZ ($L{R}_{PZ}$) and the SZ ($L{R}_{SZ}$), which are both bound between 5–70%. The objective was the minimization of the $N{O}_x$ emissions while keeping the $CO$ emissions unchanged, thus replicating Mellor’s [34] objective.
• Opt2: The selected design variables were exactly the same as the case of Opt1, while the objective was the minimization of all three pollutant emissions: $N{O}_x$, $CO$, and UHCs.
• Opt3: The selected design variables were the ones selected in the cases Opt1 and Opt2, plus L and $A_{ref}$, which can be modified within a $\pm 10\%$ range from their initial value (from the initial design).The objective, as in case Opt2, was the minimization of all three pollutant emissions: $N{O}_x$, $CO$, and UHCs.

For all four optimization scenarios, the $EI$ for each pollutant, yielded by the CRN model when evaluating each design, was incorporated in an ICAO LTO cycle analysis so as to determine the overall mass of the emissions produced for the LTO data corresponding to the CF56 engine.

The geometrical data for the final solution of each of the aforementioned optimization scenarios can be found in

Table A3. Initial and improved combustor geometrical data.

Parameter	Initial	Mellor	Opt1	Opt2	Opt3
$A_{ref}$ [m${}^2$]	0.155	0.155	0.155	0.155	0.148
L [m]	0.224	0.224	0.224	0.224	0.224
${\varphi }_{PZ}$ *	0.842	0.818	0.913	1.077	1.089
${\varphi }_{SZ}$ *	0.656	0.639	0.605	0.589	0.580
$A{R}_{PZ}[\%]$	38.69	39.81	35.67	30.48	29.92
$A{R}_{SZ}[\%]$	10.94	11.17	18.15	25.43	26.21
$L{R}_{PZ}[\%]$	20.29	15.19	14.18	19.29	35.12
$L{R}_{SZ}[\%]$	19.79	22.92	30.99	48.76	41.46

* In cruise.

of Appendix C, while the computed LTO emissions for each final design are presented in

Table 1. LTO cycle emission results for each final design of each optimization case.

	Initial	Mellor	Opt1	Opt2	Opt3
$E{I}_{NOx}$ TO	78.02	75.56	49.68	30.85	26.06
$E{I}_{NOx}$ Cl	59.66	50.26	49.29	21.39	18.83
$E{I}_{NOx}$ Ap	2.14	1.67	3.51	12.83	18.19
$E{I}_{NOx}$ Id	0.31	0.20	0.47	1.47	2.09
LTO $N{O}_x$ (gr)	12,999.3	11,469.2	10,102.0	6090.6	6042.0
$E{I}_{CO}$ TO	1.83	1.46	1.15	1.45	1.88
$E{I}_{CO}$ Cl	1.51	1.19	1.00	1.41	1.91
$E{I}_{CO}$ Ap	2.04	2.30	2.48	3.33	3.84
$E{I}_{CO}$ Id	25.33	34.55	25.77	18.39	17.03
LTO $CO$ (gr)	5131.0	6778.7	5129.8	3940.2	3824.6
$E{I}_{UHC}$ TO	4.0 $\times {10}^{-6}$	7.0 $\times {10}^{-6}$	5.0 $\times {10}^{-6}$	6.0 $\times {10}^{-6}$	6.0 $\times {10}^{-6}$
$E{I}_{UHC}$ Cl	6.7 $\times {10}^{-5}$	1.1 $\times {10}^{-4}$	6.7 $\times {10}^{-5}$	2.6 $\times {10}^{-5}$	2.2 $\times {10}^{-5}$
$E{I}_{UHC}$ Ap	0.040	0.052	0.036	0.015	0.013
$E{I}_{UHC}$ Id	2.077	3.360	1.516	0.233	0.134
LTO $UHC$ (gr)	384.0	620.5	279.9	44.4	25.7

. The corresponding relative difference of the LTO analysis for each final geometry is shown in the following Figure 12.

Based on Table 1 and Figure 12, it can be deduced that, for this specific application, Mellor’s low $N{O}_x$ algorithm was capable of decreasing the LTO $N{O}_x$ emissions by approximately 11.8%, yet this was achieved by simultaneously increasing the $CO$ and UHC emissions. Adversely, the Opt1 analysis solution showcased the ability to achieve a better $N{O}_x$ reduction (22.3%) by also decreasing the UHC by 27.1% and keeping the $CO$ practically unchanged. Still, even better LTO emission reductions could be achieved for all three pollutants through optimization scenarios Opt2 and Opt3.

The improved burner geometries are shown in Figure 13 and Figure 14. It is clear that, regarding the air mass flow in each zone, greater reductions were achieved by decreasing the hole diameter before the SZ and by increasing it before the DZ. When it comes to the length of each zone, the only safely assumed conclusion is that the increase of the length of the SZ yielded a consistent reduction trend with respect to the three examined pollutants.

5. Discussion and Conclusions

A simplified model implementing the CRN approach for combustion chamber modeling has been presented. A particular feature that is different from most approaches in the literature is that each zone of the burner was simulated through one PSR formulating three-element serial PSRs for the whole burner. This served to decrease the complexity, the computational cost, and the amount of engine-specific data needed to operate it mainly concerning the partial mixing within the PZ. The limitations of the proposed model’s PSRs layout were circumvented through the implementation of a tuning method, thereby assuming expected values of equivalence ratios for each zone at cruise and an igniter mechanism, which were used to emulate the spark ignition to kick-start the combustion of the fuel.

The implemented PSRs configuration of the model provided the ability to utilize effectively detailed combustion mechanisms without yielding a prohibitive computational cost for each simulation. The employment of detailed mechanisms enabled the model to simultaneously generate predictions for the combustion efficiency and a wide range of its pollutant emissions in contrast to most similar research in the literature that focus mostly on $N{O}_x$ and $CO$ estimations. Its generality and accuracy were put to the test against three different correlation methods (two for predicting the combustion efficiency and one for estimating the $E{I}_{NOx}$) and three actual aeroengine configurations. In all cases, the proposed model correctly simulated the expected performance and emission trends compared to both the correlation models and the reference ICAO emission data, without taking into account the volume change of each zone under partial loads. Even so, the consideration of the zone volume as a separate degree of freedom, in each examined condition during the emission trend matching process, may be encompassed in future work on the proposed model to achieve greater prediction accuracy without modifying the current CRN configuration.

The model was also integrated within a combustion chamber preliminary design process and managed to generate combustor designs with better LTO emission performance compared to Mellor’s correlation-based, algorithm for low $N{O}_x$ [34]. The design application validates the proposed model’s sensibility to geometric modification, thus opening the possibility of further use in overall engine design processes. The objective of our team is to pursue the extension of our current state-of-the-art aeroengine preliminary design methods [16] with the inclusion of design-sensitive physics-based models for the simulation of the combustion chamber, such as the one proposed in the current work.

Overall, the proposed model is believed to be capable of fulfilling the current demand for reliable yet generalist computational tools for a more in-depth simulation of the combustion process when compared to the widely applied correlation methods.