2.2.1. Model description

Single-zone models capture only the mean cylinder temperature and composition and are thereby insufficient to describe hot cylinder regions that are crucial for $\mathrm{NO}_x$ formation. In contrast to multi-zone models, the two-zone approach simplifies the combustion chamber into two distinct thermodynamic domains: the unburned and the burned (combusted) zones. These regions are separated by an infinitely thin flame front, which is typically not explicitly modelled due to its negligible volume and limited direct impact on the global kinetics of $\mathrm{NO}$ formation.

The primary differentiation among two-zone methodologies lies in their initialisation of these domains and the mechanisms used to describe their interaction. Li et al. [Reference Li, Zhu and Men26] present a two-zone reaction-based combustion model that accounts for both heat and mass exchange between the zones. The model described by Kozuch [Reference Kozuch27] includes a detailed mixing model, solving a global homogeneous isotropic $k-\epsilon$ model to resolve cylinder turbulence. Similar approaches to calculate in-cylinder turbulence to determine mixing are given by Andersson et al. [Reference Andersson, Johansson, Hultqvist and Noehre28] and Brückner et al. [Reference Brückner, Pandurangi, Kyrtatos, Bolla, Wright and Boulouchos29]. Heider et al. [Reference Heider, Woschni and Zeilinger30] instead describe a more empirical two-zone model that only requires few engine-specific input parameters. Oppenauer and Del Re [Reference Oppenauer and Del Re31] combine the Heider approach with the parallel mixing model of Kozuch [Reference Kozuch27]. The authors also highlight that pure physical models are difficult to adapt, especially to modern diesel engines, as the required calibration parameters vary during operation.

In the following, the two-zone model implementation is described, as it proved sufficiently accurate while requiring a minimum amount of unknown setup parameters. The fundamental logic originates from Heider et al. [Reference Heider, Woschni and Zeilinger30], with the implementation most closely resembling Kaiser [Reference Kaiser17].

Figure 1.

Schematic of the two-zone mass balance.

Figure 1 visualises the two-zone mass balance, which is described in mathematical terms through Equations (5–9). The cylinder mass is partitioned into a burned ( $1$ ) and unburned zone ( $2$ ),

(5) \begin{align} m_{total} = m_1 + m_2, \end{align}

where $m_1$ and $m_2$ are the two zone masses. The fuel-air ratio of the burned zone $FAR_1$ is defined by the air-excess ratio $\lambda _0$ and the stoichiometric $FAR_{stoich}$ :

(6) \begin{align} FAR_1 = \dfrac {FAR_{stoich}}{\lambda _0} \end{align}

The effect of vitiated inlet air and incomplete scavenging is represented through background fuel and water contents in the incoming charge. These background quantities define the initial $FAR_0$ and $WAR_0$ , which are constant across both zones. Directly injected fuel $m_{f,inj}= m_f - m_{f,0}$ and water $m_{w,inj} = m_w - m_{w,0}$ are added only to the combusted zone, resulting in the burned zone composition $FAR_1$ and $WAR_1$ . Consequently, a change in $FAR_0$ alone does not alter $FAR_1$ ; instead, the amount of freshly burned fuel and the corresponding heat release are determined by the difference $FAR_1-FAR_0$ . This difference defines the air mass $m_{a,1}$ participating in combustion.

(7) \begin{align} m_{a,1} = \dfrac { m_{f,inj} }{ FAR_1 - FAR_0 } \end{align}

(8) \begin{align} WAR_1 = \dfrac {m_{w,inj}}{m_{a,1}} + WAR_0 \end{align}

(9) \begin{align} m_1 = m_{a,1} \cdot (1 + FAR_1 + WAR_1) \end{align}

Heider et al. [Reference Heider, Woschni and Zeilinger30] give a simple model for the temperature difference between the burned and unburned zones.

(10) \begin{align} T_1 - T_{2} = B(\varphi ) \cdot A \end{align}

The shape function $B(\varphi )$ describes the evolution of the temperature difference along the crank angle $\varphi$ and gives an estimate for the energy interaction between the two zones based on pressure evolution,

(11) \begin{align} B(\varphi ) = 1 - \dfrac { \int _{\varphi _{SOC}}^{\varphi } ( p(\varphi )-p_0(\varphi ) ) \, m_1(\varphi ) \, d\varphi }{ \int _{\varphi _{SOC}}^{\varphi _{EOC}} ( p(\varphi )-p_0(\varphi ) ) \, m_1(\varphi ) \, d\varphi }, \end{align}

with $p(\varphi )$ being the cylinder pressure and $p_0(\varphi )$ the pressure without combustion. SOC and EOC refer to start and end of combustion.

The factor $A$ is described by Heider et al. [Reference Heider, Woschni and Zeilinger30] as an engine-specific calibration temperature difference. Oppenauer and Del Re [Reference Oppenauer and Del Re31] extend $A$ to be variable with swirl and EGR. Kaiser [Reference Kaiser17] describes the same approach as chosen here, where $A$ is defined as a temperature difference based on the adiabatic flame temperature.

(12) \begin{align} A = c_A \, (T_{flame} - T_{SOC}) \end{align}

The adiabatic flame temperature $T_{flame}$ is established via an enthalpy balance between $SOC$ conditions and the flame conditions given by $FAR_1$ , $WAR_1$ and the water injection temperature $T_{w,inj}$ . Compared with the original Heider formulation, this definition reduces the calibration burden on $c_A$ , since the temperature difference $T_{flame}-T_{SOC}$ already accounts for changes in mixture composition, vitiation, water addition and SOC thermodynamic conditions.

To calculate the zone temperatures, an enthalpy balance ensures energy conservation,

(13) \begin{align} m_{total} \, h_{total} = m_{1} \, h_{1} + m_2 \, h_2, \end{align}

with the specific enthalpies $h$ resulting from the temperature, pressure and composition in the zones.

2.2.2. Calibration and validation

The described methodology results in two calibration parameters. The air-excess ratio $\lambda _0$ captures the amount of air added to the burned zone, with $\lambda _0=1$ describing stoichiometric combustion. An increase of $\lambda _0$ will initially result in an increase in $\mathrm{NO}$ as the increased availability of oxygen outweighs the decrease in flame temperature. As $\lambda _0$ moves further away from stoichiometric conditions, the predictive capability of a single representative value decreases since in a real engine, $\mathrm{NO}$ formation is governed by a distribution of local equivalence ratios, temperatures and radical concentrations rather than by a spatially uniform-burned zone composition. Rich and lean regions will differently impact $\mathrm{NO}$ formation due to the availability of radicals. For direct injection, a single $\lambda _0$ near stoichiometry has been shown to be capable of providing a good $\mathrm{NO}$ estimate. For premixed conditions, the predictions will be less accurate.

The parameter $c_A$ scales the initial zone temperature difference, with $c_A=1$ resulting in $T_1(\varphi _{SOC})=T_{flame}$ . Brückner et al. [Reference Brückner, Pandurangi, Kyrtatos, Bolla, Wright and Boulouchos29] found that actual peak temperatures in conventional diesel combustion typically remain approximately $150 \, K$ below the adiabatic flame temperature. However, the authors also find that peak temperatures approach or even surpass this limit in premixed-dominated regimes. This convergence is primarily driven by reduced mixing rates and the compression heating of post-flame gases, which counteracts the thermal dissipation typically seen in diffusion flames. Andreasen and Mayer [Reference Andreasen and Mayer32] reduce the calculated flame temperature by approx $100\,K$ via heat radiation. During validation $c_A=0.9$ was found to result in decent $\mathrm{NO}_x$ trends while the resulting reduction in $A$ matched the above-cited temperature ranges.

While in this study $c_A$ was kept constant, one could make a case to modify it for a reason other than simply the difference between adiabatic flame temperature and actual hot gas temperature. Although the Heider two-zone model accounts for inter-zone energy transfer via the shape function $B$ (Equation 11), it does not explicitly resolve the finite-rate turbulent mixing or local mixture stratification. Under high-pressure, high-temperature conditions close to stoichiometry, $\mathrm{NO}$ formation may therefore approach equilibrium rapidly. In such cases, the simplified treatment of mixing and stratification may cause the calculated $\mathrm{NO}$ concentration to become limited by equilibrium during the peak-temperature period. Consequently, $\lambda _0$ is typically adjusted empirically to match experimental $\mathrm{NO}$ results. While this may suffice for part-load studies, it can result in an underestimation of $\mathrm{NO}_x$ reduction when evaluating advanced mitigation strategies, such as water injection. In such cases, the model may fail to show the expected $\mathrm{NO}$ reduction from flame temperature cooling because the predicted concentrations initially remain equilibrium-limited. As $c_A$ directly affects peak temperature, it could be used as a proxy in calibration to limit the extent to which $\mathrm{NO}$ formation becomes equilibrium-controlled. In the present study, however, it is kept constant, as satisfactory agreement with literature data was obtained and, given the uncertainty inherent to predesign studies, a physically well-founded justification for its adjustment would be difficult.

The chosen validation engines and calibration parameters are presented in Table 1. To validate the described piston $\mathrm{NO}_x$ model, crankshaft-based engines have been set up due to the lack of sufficiently detailed data on free piston engines. The performance setups were tuned to closely align with available data, such as temperature or pressure evolutions and integral metrics such as power or specific fuel consumption $SFC$ . Subsequently, $\lambda _0$ was tuned for best alignment with literature data, while $c_A$ was kept constant. In addition to the diesel part-load studies, previously shown by Kaiser [Reference Kaiser17], an $\mathrm{H2}$ engine and vitiated air studies were added.

Table 1.

Overview of piston engine setups used for validation

Table 1. Long description

A table comparing piston engine setups for validation. The table has four rows and four columns. The columns are labeled Engine, Stroke, Bore, and lambda_0. The row labels are Diesel MAN B&W 4T50ME, Diesel Large Generator, and H2 Lombardini 15LD500. Row 1: Diesel MAN B&W 4T50ME, 2, 0.5, 1.093. Row 2: Diesel Large Generator, 2, 0.58, 1.060. Row 3: H2 Lombardini 15LD500, 4, 0.087, 1.005.

¹[Reference Andreasen and Nyggard6, Reference Raptotasios, Sakellaridis, Papagiannakis and Hountalas15, Reference Scappin, Stefansson, Haglind, Andreasen and Larsen33].

²[Reference Lamaris, Hountalas, Zannis and Glaros14].

³[Reference Sfriso, Berni, Fontanesi, d’Adamo, Breda, Teodosio, Frigo and Antonelli34, Reference Simi35].

Figure 2 shows the load and EGR (at $75\,\%$ load) study for the MAN B&W 4T50ME-X. The model captures the power-specific $\mathrm{NO}_x$ emissions well across various loads. However, an underestimation of specific emissions decrease with increasing load can be observed. This decrease is likely a result of the improved mixing at higher engine speeds, which is not captured by the model. For the EGR study, the inlet $FAR$ was calculated for each $EGR$ value on a dry volume basis, as described by Andreasen and Nyggard [Reference Andreasen and Nyggard6]. The model captures the trend generally well, with the exception of high $EGR$ values. The experimental setup utilised an exhaust gas scrubber that cools the recirculated air and removes significant amounts of water, resulting in dry EGR. Because the $FAR$ definition of the model inherently includes combustion-product water, it cannot currently simulate a fully dry recirculation path. This is the likely reason for the higher emission reduction predicted with $EGR$ , where the additional increase in specific heat capacity of the recirculated water has a substantial effect on the flame temperature and thereby $\mathrm{NO}$ formation. For the aeronautical application only wet EGR will be considered.

Figure 2.

$\mathrm{NO}_x$ predictions for the MAN B&W 4T50ME-X (a) load and (b) EGR study.

Figure 3a shows the $\mathrm{NO}_x$ emission predictions for the MAN B&W 4T50ME-X operating with water-in-fuel emulsion at $75\,\%$ load. In contrast to the experiment, water is not premixed with the fuel in the model; instead, it is introduced into the cylinder at $298.15$ K, with the mass flow determined by the fuel flow and the specified water-fuel ratio ( $WFR$ ). The model predicts a linear increase in specific fuel consumption ( $SFC$ ) due to the thermodynamic penalty of heating the water to cylinder temperatures and the associated latent heat of vaporisation (Figure 3b). The experimental results [Reference Andreasen and Nyggard6] show a less pronounced $SFC$ increase at moderate $WFR$ and even report $SFC$ improvements at lower loads. The authors suggest this effect could result from an effective fuel nozzle optimisation due to the added mass flow. At low $WFR$ , the model predicts a gradual $\mathrm{NO}_x$ reduction, highlighting the sensitivity of the model to the transition between equilibrium-controlled and kinetically limited $\mathrm{NO}$ formation. At higher $WFR$ , $\mathrm{NO}$ formation becomes clearly kinetically limited, and the model’s predicted reduction slope closely matches the experimental trends.

Figure 3.

$\mathrm{NO}_x$ and $SFC$ predictions for the MAN B&W 4T50ME-X water-in-fuel emulsion study.

Figure 4a depicts the load study of the large diesel generator described by Lamaris et al. [Reference Lamaris, Hountalas, Zannis and Glaros14]. The model shows good alignment at high load but underestimates the $\mathrm{NO}_x$ emissions at lower load. For the small $\mathrm{H_2}$ engine at $2500$ rpm shown in Figure 4b, the model captures the trend of increasing $\mathrm{NO}_x$ emissions with equivalence ratio but fails to capture the high emissions for $\phi =0.6$ . Sfriso et al. [Reference Sfriso, Berni, Fontanesi, d’Adamo, Breda, Teodosio, Frigo and Antonelli34] show that the early injections in this engine result in a wide distribution of local equivalence ratio. The distributions change with global $\phi$ and engine speed. For $\phi =0.6$ at high engine speed, the local equivalence ratios reach up to stoichiometry. This shows one limitation of the model, where the single equivalence ratio choice for the combusted zone fails to capture the wide distribution of the actual conditions.

Figure 4.

$\mathrm{NO}_x$ prediction for (a) the large generator load study [Reference Lamaris, Hountalas, Zannis and Glaros14] and (b) $\mathrm{H_2}$ 4-stroke Lombardini 15LD500 [Reference Simi35].

Overall, the piston $\mathrm{NO}_x$ model shows to be capable of predicting general $\mathrm{NO}_x$ trends in the studied piston engines well, including water injection and EGR. For the large two-stroke diesel engines, good alignment with a higher $\lambda _0$ of $1.06$ and $1.093$ was found. These result in a higher $\mathrm{NO}_x$ than at the stoichiometric setting due to air entrainment. However, the small four-stroke $\mathrm{H}_{2}$ engine calibration resulted in a near-stoichiometric combusted zone with $\lambda _0=1.005$ . The previous application of this two-zone model [Reference Kaiser17] showed two smaller four-stroke engines with $\lambda _0$ at $1.01$ and $1.03$ . While these ranges might imply a dependency on the engine cycle (two-stroke vs. four-stroke), here it is argued that engine scale and speed are the more physically relevant factors. In large-bore engines, larger physical mixing scales and slower speeds favour a broader, leaner bulk $\mathrm{NO}_x$ production zone. Conversely, in small-bore engines, regardless of cycle architecture, rapid mixing tends to restrict high-temperature $\mathrm{NO}_x$ production to narrow, near-stoichiometric interfaces.

Consequently, for the smaller-bore two-stroke FDP engine evaluated in this study, $\lambda _0=1$ is adopted. This choice prioritises the physical similarities in mixing dynamics with small four-stroke engines over the nominal similarity in cycle architecture with large diesel engines. While this hypothesis provides a consistent physical basis for this study, further validation, preferably with a two-stroke, small-bore, free-piston engine, is recommended. Within the defined scope, the model is considered sufficient to capture the general trends in the investigated CCE.

2.3.1. Model description

Similar to the piston engine, the fundamental difficulty of gas turbine $\mathrm{NO}_x$ modelling lies in the extreme temperature sensitivity of the thermal NO mechanism. The rate of formation scales exponentially with temperature; consequently, the majority of NO is produced in local regions where the equivalence ratio is near stoichiometric ( $\phi \approx 1$ ), even if the global mixture is lean. Traditional zero-dimensional (0D) models that assume a perfectly stirred reactor (PSR) based on average inlet conditions often fail to capture this reality. They tend to substantially underpredict NO formation at lean operating points.

To resolve the combustor mixture without solving the 3D Navier–Stokes equations, stochastic descriptions of the flow field have been introduced [Reference Fletcher and Heywood36, Reference Heywood and Mikus37]. This methodology assumes that while the spatial location of a fuel pocket is unknown (0D/1D approximation), the probability of finding a fluid element at a specific equivalence ratio can be mathematically defined. The key aspects for these models are the stochastic description of the local equivalence ratio and its time evolution.

The stochastic distribution of the local equivalence ratio ( $\phi$ ) can be described by a probability density function (PDF). Heywood and Mikus [Reference Heywood and Mikus37] assume a Gaussian distribution, while other distributions such as $\beta$ -PDFs could offer more flexibility if supporting data is available. Heywood and Mikus also consider the stochastic distribution of the local residence time in the primary zone. The fluid homogenises with time, which is of key interest to $\mathrm{NO}_x$ formation. This can be captured by sectioning of the combustion chamber into zones with different roles. Heywood and Mikus [Reference Heywood and Mikus37] consider a primary zone without any mixing and a secondary zone where cooling air dilutes the richer primary zone gases and quickly freezes the NO formation process. Therefore they only consider primary zone formation, while noting that this assumption fails if a substantial part of the primary zone flow is fuel rich. Fletcher and Heywood [Reference Fletcher and Heywood36] consider an intermediate zone after the primary zone where homogenisation of the standard deviation takes place proportional to $\sqrt {t}$ . In the following dilution zone no further NO formation is assumed. Dederichs et al. [Reference Dederichs, Zarzalis, Habisreuther, Beck, Prade and Krebs38] assume an exponential decrease of the standard deviation with respect to a characteristic time, which was established by means of a reacting large eddy simulation (LES).

The strength of such a stochastic approach is its ability to decouple the macro-scale fluid dynamics from the micro-scale mixing and describing these with a limited number of physical parameters. The assumed PDF can easily be modified to allow for more control over the flow mixedness. This allows quick calibration to different combustor designs, while keeping computational costs low. However, no description for an appropriate cooling air addition model, which is critical to e.g. rich-burn, quick-mix, lean-burn (RQL) combustors, was found in the literature. In the following an implementation for such a stochastic model is described, including the novel adaption to RQL-type combustors. The foundational model could also be adapted to other advanced combustion concepts, such as lean direct injection or staged combustion, which would be an interesting area for future research.

This model simplifies the spatial considerations by assuming only one zone with both homogenisation of the mixture and mixing with the cooling air. To account for the imperfect mixing of fuel and air, the methodology employs a Gaussian (normal) PDF to describe the distribution of local equivalence ratio $\phi$ within the combustor, following the description by Heywood and Mikus [Reference Heywood and Mikus37],

(14) \begin{align} f(\phi ) = \dfrac {1}{\sigma \sqrt {2 \pi }} \cdot exp \left [ - \dfrac {(\phi -\bar {\phi })^2}{2 \sigma ^2} \right ], \end{align}

where $\sigma$ is the standard deviation. The mixing quality is quantified by the unmixedness parameter $s$ , defined as the ratio of the standard deviation $\sigma$ to the mean equivalence ratio $\bar {\phi }$ :

(15) \begin{align} s = \dfrac {\sigma }{\bar {\phi }} \end{align}

The evolution of mixture inhomogeneity is modelled as an exponential decay, described by a characteristic mixing timescale $\tau _{mix}$ and the initial unmixedness $s_0$ [Reference Dederichs, Zarzalis, Habisreuther, Beck, Prade and Krebs38]:

(16) \begin{align} s(t) = s_0 \cdot exp \left ( - \dfrac {t}{\tau _{mix}} \right ) \end{align}

The mean equivalence ratio $\bar {\phi }$ is determined by the injected fuel mass and the cumulative air addition profile $\alpha (t)$ . Only a part of the inlet air is directly routed into the primary zone. This split is described by the assumed inlet air split $\beta$ :

(17) \begin{align} \beta = \dfrac {w_{t=0}}{w_{inlet}} \end{align}

Here, $w_{t=0}$ refers to the mass flow in the combustor at start of simulation and $w_{inlet}$ is the inlet mass flow from the upstream component.

To capture the quick quenching required in RQL combustors, the cooling air dilution is modelled with a squared exponential based on a characteristic dilution timescale $\tau _{dil}$ . This phenomenological formulation assumes all cooling air to be added within this quick mixing. In real combustors a substantial fraction of the cooling air is added continuously for dilution towards the combustor exit, but the assumption is made that this has little influence on $\mathrm{NO}$ formation, as it is frozen by the quenching. The cumulative air addition profile $\alpha$ at time $t$ is defined as

(18) \begin{align} \alpha (t) = \left\{\begin{array}{l@{\quad}l} \beta & {} {\rm for\ } 0 \leq t \leq \tau _{primary} \\[5pt] \beta + (1 - \beta ) \cdot \left ( 1 - exp\left [- \left (\dfrac {t-\tau _{primary}}{\tau _{dil}}\right )^2 \right ] \right ) {} & {\rm for\ } t \gt \tau _{primary} \end{array}\right \}, \end{align}

with $\tau _{primary}$ describing the length of the rich, primary zone ahead of quenching.

The dilution $\tau _{dil}$ , the mixing timescale $\tau _{mix}$ and the primary zone time $\tau _{primary}$ are dependent on engine operation. As done by Dederichs et al. [Reference Dederichs, Zarzalis, Habisreuther, Beck, Prade and Krebs38], the timescales are estimated with calibration parameters $c_{dil}$ , $c_{mix}$ or $c_{primary}$ , the bulk flow velocity $u_{bulk}$ and a characteristic length (here the nozzle diameter is used).

(19) \begin{align} \tau _i = c_i \, \dfrac {d_{nozzle}}{u_{bulk}} = c_i \, \dfrac {\rho \, A_{combustor} \, d_{nozzle}}{\dot {m}} \end{align}

$A_{combustor}$ refers to the annular cone flow area after the injector including the cooling ducts. The nozzle diameter $d_{nozzle}$ serves as a characteristic diameter for the premixing passage.

2.3.2. Discretisation

The flow field is discretised both temporally and stochastically with respect to the local fuel–air equivalence ratio $\phi$ . The stochastic discretisation segments the PDF given by Equation (14) into a set of parcels, with each parcel representing the same fixed fraction of the fuel mass flow (see Figure 5a). The respective air mass then results from the local $\phi$ obtained from the inverse cumulative distribution function of $f(\phi )$ . Total mass conservation is ensured through appropriately scaled parcels.

Figure 5.

Joule combustor $\mathrm{NO}_x$ model (a) discretisation visualisation and (b) discretisation sensitivity.

The discretisation sensitivity is displayed in Figure 5b. Noteworthy is the high sensitivity to the timestep highlighting the time-sensitive nature of thermal $\mathrm{NO}$ formation. The number of parcels is less impactful on $\mathrm{NO}_x$ with little difference between $30$ and $40$ parcels. Accordingly $30$ parcels and timestep intervals of $0.1$ ms are selected as they lie within $4\%$ of the highest-resolved setting while minimising computational cost. Due to the sensitive calibration, this deviation is secondary. Calculation times are then typically well within a second.

The temperature in each parcel and timestep is calculated from a flame temperature calculation via enthalpy balance to inlet conditions given the local $\phi$ . Effects of vitiated air, as in the case of the precombusted air in the CCE, are captured by the inlet conditions, which propagate into the combustor. The thermal $\mathrm{NO}$ formation at each timestep and parcel is calculated as previously described in subsection 2.1.

2.3.3. Calibration and validation

The proposed model integrates outputs from a gas turbine performance simulation as primary boundary conditions. Specifically, it requires the inlet mass flow, pressure, temperature, FAR and water–air ratio (WAR), alongside the total fuel mass added. Furthermore, the chemical and mixing evolution within the combustor is defined by a number of parameters: the inlet air split $\beta$ , representing the initial equivalence ratio; the initial unmixedness $s_0$ ; and the characteristic timescales for internal homogenisation $\tau _{mix}$ , dilution $\tau _{dil}$ and the primary zone $\tau _{primary}$ . While these timescales are derived from specific calibration constants, as described by Equation (19), their physical significance and representative values from the literature are discussed in the following section.

The inlet air split $\beta$ represents the fraction of total combustor air directed into the primary reaction zone. It establishes the initial stoichiometry, dictating the base level of thermal $\mathrm{NO}$ production before any secondary dilution or internal mixing occurs and whether later the stoichiometric region needs to be crossed. The upper bound of $\beta$ is typically constrained by the lean flammability limit. In staged combustion systems, a lightly loaded primary zone may operate at an equivalence ratio of approximately $0.8$ during low-power conditions to ensure high combustion efficiency and minimise CO and unburned hydrocarbon (UHC) emissions. Conversely, at full-load conditions, the primary zone $\phi$ is often lower, at around $0.6$ , to mitigate $\mathrm{NO}_x$ formation as the secondary zone becomes active ([Reference Lefebvre and Ballal39], pp. 24–25). RQL combustors show an even higher initial equivalence ratio (1.2–1.6 [Reference Lefebvre and Ballal39]) and rely on quick mixing. Smith [Reference Smith40] reports a primary zone equivalence ratio of about $2$ for the V2500-A5 combustor at 100% load, with only 19% and 7% of total combustor air going to the primary zone and the nozzle, respectively. It was found that such rich primary zones would typically result in $\beta \lt 0.4$ .

The unmixedness $s_0$ accounts for the spatial fluctuations in the fuel-air mixture. Due to the exponential dependence of thermal $\mathrm{NO}_x$ formation on temperature, local hot spots near stoichiometric regions produce significantly higher emissions than a perfectly homogeneous mixture at an identical mean temperature. Consequently, $s_0$ is essential for capturing the influence of these high-temperature eddies. The impact of initial unmixedness is fundamentally coupled with the mean equivalence ratio and, by extension, the air split $\beta$ . While a lower limit of $s_0=0$ represents a perfectly premixed state, which is difficult to achieve in aero-engine combustors, higher values characterise non-premixed systems. For instance, Heywood and Mikus [Reference Heywood and Mikus37] used $s_0=0.5$ for predicting GT-309 NO emissions ( $\phi (t=0)=0.81$ at 100% load). Smith [Reference Smith40] reported $s_0$ in the range of $0.2$ to $0.4$ for rich Pratt & Whitney injectors. Smith identified two primary drivers for mixing quality: the injection mechanism and the total energy available for mixing, typically represented by the pressure drop across the combustor. It is expected that unmixedness increases as engine load decreases, as both the turbulence intensity from the fuel injectors and the available pressure drop are reduced.

The mixing timescale $\tau _{mix}$ represents the time required for fuel and air to reach a homogeneous state. Heywood and Mikus [Reference Heywood and Mikus37] assume no change in umixedness with time, as they consider a stochastic distribution of the primary zone residence time. Dederichs et al. [Reference Dederichs, Zarzalis, Habisreuther, Beck, Prade and Krebs38] estimated a characteristic mixing time of $5.03$ ms, while noting lower numbers from the literature. Beyond direct estimation, $\tau _{mix}$ (in conjunction with $s_0$ ) can be calibrated by evaluating species or temperature distributions at the combustor exit plane. These distributions are frequently documented in literature due to their critical influence on high-pressure turbine cooling requirements and nozzle guide vane durability. Alternatively, as a first-order approximation, the timescale may be estimated using the geometric combustor length and the bulk average flow velocity.

The dilution timescale $\tau _{dil}$ defines the rate at which secondary air is introduced. In RQL architectures, rapid dilution is critical for minimising residence time within the high $\mathrm{NO}_x$ window, where fuel-rich pockets inevitably pass through stoichiometric conditions before reaching the final lean state. Pegemanyfar [Reference Pegemanyfar41] performed Reynolds-averaged Navier–Stokes (RANS) simulations of an RQL combustor, observing that quenching occurred between $30$ and $40$ % of the total combustor length. However, the study also identified hot spots resulting from near stoichiometric combustion significantly further downstream. This implies a relatively slow mixing process which, within the context of the current model, would be characterised by a larger $\tau _{dil}$ , leading to elevated $\mathrm{NO}_x$ emissions.

The primary zone residence time $\tau _{primary}$ describes the duration until cooling air is added. Pegemanyfar [Reference Pegemanyfar41] provides residence time distributions at the primary zone exit plane with a mean value of approximately $3$ to $4$ ms. For the GT-309, Heywood and Mikus [Reference Heywood and Mikus37] describe a decreasing mean primary zone residence time from $11.2$ to $6.8$ ms for $50$ and $100$ % load, respectively.

Model calibration and validation with an in-service combustor requires some description of key parameters. For the combustor of the V2500 turbofan engine sufficient literature statements were found. The model setup for the V2533-A5 is summarised in Table 2. All other inputs result from cycle calculation. With most model parameters originating from literature statements, two parameters were adjusted to calibrate to emission measurements. The dilution timescale parameter $c_{dil}$ was used as the main calibration parameter resulting in the overall scale of $\mathrm{NO}_x$ predictions. The primary zone length parameter $c_{primary}$ was calibrated for the $30\%$ thrust point, as it is near stoichiometric in the primary zone.

Table 2.

Model setup for the V2533-A5 RQL combustor

Table 2. Long description

A table with six rows and five columns detailing the model setup for the V2533-A5 RQL combustor. The columns are labeled Variable, Unit, Value, and Justification. The rows list specific variables such as beta, s0, cmix, cdil, cprimary, dnozzle, and Acombustor with their respective units, values, and justifications. Each row provides a detailed explanation for the chosen value of the variable, referencing literature and measurements.

Figure 6.

Sea-level static LTO cycle $\mathrm{NO}_x$ prediction of the (a) V2533-A5 compared to ICAO databank and (b) near future setup.

Figure 6a shows the sea-level static (SLS) $\mathrm{NO}_x$ emissions at $100$ , $85$ , $30$ and $7 \, \%$ of rated thrust (landing and take-off operations (LTO)) of the V2533-A5 engine (strongest variant of the V2500 engine family). Measurements, reported in the International Civil Aviation Organization (ICAO) emissions databank [42], are plotted alongside. The model is able to accurately reproduce the engine emissions across the different thrust levels, with some deviation at the lowest thrust setting.

Figure 7.

$\mathrm{NO}_x$ model evolutions including stochastic ranges for the V2533-A5 at full and $30 \, \%$ thrust.

For the combustor analysed here, the dilution timescale resulting from $c_{dil}$ is most sensitive on $\mathrm{NO}_x$ emissions. To accurately replicate the V2533-A5 $\mathrm{NO}_x$ measurements, a relatively short primary zone and very quick quenching was required. The evolutions of the equivalence ratio, temperature and the resulting $\mathrm{NO}_x$ is given in Figure 7 for full and $30\,\%$ thrust. At full thrust the primary zone shows the described $\bar {\phi }=2$ with a wide distribution due to initial unmixedness. The temperatures and $\mathrm{NO}$ formation peaks during the stoichiometric crossing in the quenching phase. The $30\,\%$ thrust point shows a near stoichiometric primary zone, resulting in quenching quickly stopping $\mathrm{NO}$ formation.

Comparison with the findings of Pegemanyfar [Reference Pegemanyfar41] reveals some important differences. Specifically, at full thrust, the temperature peak due to stoichiometric crossing occurs around $1$ ms after the initial combustion in the primary zone, a transition much more rapid than reported by Pegemanyfar [Reference Pegemanyfar41]. This discrepancy may stem from inherent differences in combustor geometry or signify specific model limitations. While the local equivalence ratio is modelled stochastically, the dilution is not. The rich jet does not uniformly experience quenching. However, since this non-uniformity is expected to exert a linear influence on $\mathrm{NO}_x$ formation, the current model presents a simplified yet accurate approach, as the calibration of $c_{dil}$ compensates the simplified physics.

Figure 8 shows the model parameter sensitivities. The highest sensitivity is seen with the parameters describing the initial state in the combustor: $\beta$ and $s_0$ . With these factors fixed by the combustor technology, the main drivers are the dilution and primary zone length timescales $c_{dil}$ and $c_{primary}$ depending on the thrust case. As the $30$ % thrust case has a near-stoichiometric primary zone, it is most affected by a variation of the primary zone setup. Under these rapid-quenching conditions, the homogenisation rate ( $\tau _{mix}$ ) becomes a secondary driver for $\mathrm{NO}_x$ production. With the relatively slow rate ( $c_{mix}=5$ ) there is still considerable variation in the local temperature towards the combustor exit plane.

Figure 8.

Joule combustor $\mathrm{NO}_x$ model parameter sensitivity for the V2533-A5 at full and $30$ % thrust.

2.3.4. Setup for near-future combustor

Due to limited data availability, the presented model was initially calibrated to the V2533-A5 engine. To account for the future entry-into-service of the investigated engine architecture, the model was subsequently adjusted to represent a future low- $\mathrm{NO}_x$ combustor. This was achieved by first calibrating an empirical $\mathrm{NO}_x$ correlation to the V2533-A5 engine and then applying a literature-based technology factor representative of a near-future combustor. The combustion-zone timescale parameters were subsequently adjusted until the physics-based model reproduced the $\mathrm{NO}_x$ formation implied by this technology factor.

A widely used empirical correlation for the estimation of aero-engine $\mathrm{NO}_x$ is described in the gas turbine simulation environment GasTurb [Reference Kurzke43], that is referenced here as the $T_3$ - $p_3$ method:

(20) \begin{align} s_{NO_x} = \left ( \dfrac {p_3}{2965 \, kPa} \right )^{0.4} \cdot exp \left ( \dfrac {T3-836\,K}{194\,K} + \dfrac {6.29-100 \cdot WAR_{amb}}{53.2} \right ) \end{align}

This correlation requires the combustor inlet temperature $T_3$ , pressure $p_3$ and optionally the ambient humidity $WAR_{amb}$ . The emission index can be estimated by multiplying the severity factor $s_{NO_x}$ with a technology factor $c_s$ :

(21) \begin{equation} EI\,NO_x = c_s \cdot s_{NO_x} \end{equation}

Kurzke [Reference Kurzke43] provides $c_s = 32$ for conventional combustors and $c_s = 23$ for dual annular combustors. For the V2533-A5 engine, calibration against the ICAO databank SLS-LTO NOx emissions yielded $c_s=18$ , which provided good agreement with the reported data (Figure 6a). The increased emissions around $30$ % thrust due to a stoichiometric primary zone are not captured by this simple empirical correlation. For near-future projected lean direct injection combustors, $c_s = 12.4$ was calculated [Reference von der Bank44]. The authors state that $70$ % of the inlet air flows through the injectors, which would enable a model setup using $\beta =0.7$ . However, for the CCE the use of a lean direct injection combustor is questionable, due to the high level of inlet vitiation and the low heat addition in the combustor during typical cruise operation. To ensure consistent flammability a rich primary zone is likely required. To account for reduced $\mathrm{NO}_x$ emissions of a near-future combustor the model was calibrated to $c_s = 12.4$ (Figure 6b) by only reducing the dilution and primary zone timescale parameters ( $c_{dil}=0.17$ , $c_{primary}=0.19$ ). It is noted, that this level of quick quenching may be unrealistically fast, and it is recommended to study the model adaptation to future combustors in more detail.

2.3.5. Effect of vitiated inlet air

Similar to the piston engine, the effect of humidified air on $\mathrm{NO}_x$ emissions is twofold: it involves the linear reduction of adiabatic flame temperature with increasing water content and a direct chemical impact on reaction kinetics through the suppression of key radical species, such as atomic oxygen. Schimek et al. [Reference Schimek, Göke and Paschereit45] experimentally demonstrated that the chemical influence is stronger than the purely thermal effect, by comparing it to a nitrogen dilution stream with identical adiabatic flame temperature. Conversely, Hiestermann et al. [Reference Hiestermann, Konle and Guillebon46] attributed most of the $\mathrm{NO}_x$ reduction to the reduction in maximum temperature, in a numerical study of methane combustion. In ultra-wet environments (up to 50% steam), Göke and Paschereit [Reference Göke and Paschereit47] demonstrated that $\mathrm{NO}_x$ reduction is driven by the chemical suppression of the thermal, $\mathrm{N_2O}$ and $\mathrm{NNH}$ pathways via radical scavenging. However, they noted that steam injection can simultaneously increase prompt $\mathrm{NO}_x$ in natural gas flames due to elevated $\mathrm{CH}$ radical concentrations.

The present study considers only premixed water vapor (e.g. originating from the piston engine), rather than direct water injection into the gas turbine, which would likely offer a higher $\mathrm{NO}_x$ reduction potential. The predicted $\mathrm{NO}_x$ reduction as a function of the water–air ratio $WAR=\dot {m}_{water}/\dot {m}_{air}$ is presented in Figure 9a, showing strong correlation with literature data. [Reference Schimek, Göke and Paschereit45] analysed the effect of steam dilution at different equivalence ratios and nozzles, leading to multiple data points for the same $WAR$ . The prediction will slightly vary depending on baseline $\mathrm{NO}_x$ , due to the exponential nature of thermal $\mathrm{NO}$ formation. The relatively high initial $\mathrm{NO}_x$ of the V2533-A5 at full thrust therefore shows a very steep decrease.

Figure 9.

Relative $\mathrm{NO}_x$ reduction with (a) WAR and (b) precombusted air.

The use of precombusted inlet air represents a potent strategy for mitigating $\mathrm{NO}_x$ emissions through various architectural implementations. In EGR systems, a fraction of the exhaust gas is recirculated to the inlet; similarly, in staged combustion, sequential combustors, such as inter-turbine burners, utilise vitiated streams. Comparable conditions are found in the CCE, where the exhaust gases of the piston engine are funnelled into the downstream combustor, leading to a significant amount of inlet air vitiation depending on system operation. The underlying physical mechanisms mirror those of water injection: a reduction in adiabatic flame temperature and a lowering of the inlet oxygen concentration.

The relative reduction in $\mathrm{NO}_x$ emissions resulting from precombusted inlet air is depicted in Figure 9b. To simplify, similar to Muley and Lear [Reference Muley and Lear8], it has been assumed that the inlet temperature remains constant (e.g. via a recirculation cooler), with adjustments made only to the inlet mass flow, the inlet FAR and the exit FAR. Compared to the study performed by Muley and Lear [Reference Muley and Lear8], the current model predicts a more significant reduction in $\mathrm{NO}_x$ . This discrepancy is likely attributable to the differing modelling approaches. While Muley and Lear [Reference Muley and Lear8] use a single-zone PSR at stoichiometric conditions, the stochastic model used in this study captures the local temperature peaks and their sensitivity to oxygen depletion, possibly thereby predicting a higher relative reduction in thermal $\mathrm{NO}$ formation. The Joule combustor in the CCE will see high inlet equivalence ratios, where the model aligns reasonably well with the literature reportings.

The stochastic model shows solid capability in capturing the $\mathrm{NO}_x$ formation in gas turbine combustors. While relying on five model parameters, most of these can be well defined from a general description of the combustor. The model is also capable of capturing the effects of vitiated air on $\mathrm{NO}_x$ emissions. For the purpose of emission estimation in the CCE, this model is deemed effective.