Design and simulation of high temperature rise triple-swirler combustor

Abstract

To meet the development needs of aeroengines for high thrust-to-weight ratios and fuel-air ratios, a high temperature rise triple-swirler main combustor was designed with a total fuel-air ratio of 0.037, utilising advanced technologies including staged combustion, multi-point injection and multi-inclined hole cooling. Fluent software was used to conduct numerical simulations under both takeoff and idle conditions, thereby obtaining the distribution characteristics of the velocity and temperature fields within the combustor, as well as the generation of pollutants. The simulation results indicate that under takeoff conditions, the high temperature rise triple-swirler combustor achieves a total pressure loss coefficient of less than 6% and a combustion efficiency exceeding 99%. Under takeoff conditions, the OTDF and RTDF values are 0.144 and 0.0738, respectively. The mole fraction of NOx emissions is 3,700ppm, while the mole fraction of soot emissions is 2.55×10⁻⁵ppm. Under idle conditions, the triple-swirler combustor maintains a total pressure loss coefficient of less than 6% and a combustion efficiency greater than 99.9%. The OTDF and RTDF values are 0.131 and 0.0624, respectively. The mole fractions of CO and UHC emissions are both 0×10⁻³²ppm at the calculation limit of Fluent software.

Nomenclature

Roman letters

$a$

absorption coefficient, m⁻¹

${a_s}$

speed of sound, m/s

$b_{nuc}^ * $

normalised radical nuclei concentration, particles ×10⁻¹⁵/kg

${D_i}$

diffusion coefficient of the component $i$ , m²/s

$f$

mixture fraction

${f_b}$

body force term, m/s²

${g_i}$

component of the gravitational vector in the $i$ th direction, m/s²

${G_b}$

generation of turbulence kinetic energy due to buoyancy, kg/(m·s³)

${G_k}$

generation of turbulence kinetic energy due to the mean velocity gradients,kg/(m·s³)

$I$

radiation intensity which depends on position ( $\vec r$ ) and direction ( $\vec s$ ), W/(m²·sr)

$k$

turbulent kinetic energy, m²/s²

${M_i}$

molar mass of the component $i$ , kg/mol

${M_t}$

turbulent Mach number

$n$

refractive index, sr^−1/2

$p$

pressure, Pa

${\Pr _t}$

turbulent Prandtl number for energy

$\dot q$

heat received by a fluid per unit mass per unit time via radiation, J/(kg·s)

$R_{nuc}^ * $

normalised net rate of nuclei generation, particles ×10⁻¹⁵/(m³·s)

${R_{soot}}$

net rate of soot generation, kg/(m³·s)

${R_u}$

universal gas constant, =8.314J/(mol·K)

$R{R_i}$

net mass generation rate of the component $i$ due to chemical reactions, kg/(s·m³)

$S$

strain rate tensor, s⁻¹

${S_k}$

user-defined source term, kg/(m·s³)

${S_m}$

source term due solely to transfer of mass into the gas phase from liquid fuel droplets or reacting particles, kg/(m³·s)

${S_u}$

user-defined source term, kg/(m³·s)

${S_\varepsilon }$

user-defined source term, kg/(m·s³)

$t$

time, s

$T$

temperature, K

$\hat u$

internal energy of the fluid per unit mass, J/kg

$v$

velocity, m/s

${Y_M}$

contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, kg/(m·s³)

${Y_i}$

mass fraction of the component $i$

Greek letters

$\beta $

coefficient of thermal expansion, K⁻¹

$\Gamma $

surface forces(pressure and viscous forces), N/m²

${\delta _{ij}}$

Kronecker delta function

$\varepsilon $

turbulent dissipation rate, m²/s³

$\lambda $

thermal conductivity, W/(m·K)

$\mu $

dynamic viscosity, Pa·s

${\mu _l}$

laminar viscosity, Pa·s

${\mu _t}$

turbulent viscosity, Pa·s

$\rho $

density, kg/m³

$\sigma $

Stefan-Boltzmann constant, =5.669×10⁻⁸W/(m²·K⁴)

${\sigma _k}$

turbulent Prandtl numbers for $k$

${\sigma _{nuc}}$

turbulent Prandtl number of nuclei transport

${\sigma _s}$

scattering coefficient, m⁻¹

${\sigma _\varepsilon }$

turbulent Prandtl numbers for $\varepsilon $

$\Phi $

phase function

${\Omega}'$

solid angle, sr

Sub-Superscripts

${\rm{ - }}$

Favre mean (density-averaged) variable

$ox$

oxidizer

$ \to $

vector

$3$

diffuser inlet section

$4$

flame tube outlet section

1.0 Introduction

As aeroengine performance continues to advance, the design of the main combustion chamber is being pushed towards achieving a higher thrust-to-weight ratio, greater temperature rise, increased fuel/air ratio (FAR) and reduced pollution. In recent decades, significant advancements have been made in both military small-bypass-ratio turbofan engine combustion chambers and commercial large-bypass-ratio turbofan engine combustion chambers. The fourth-generation military small-bypass-ratio turbofan engines, such as the F119 and EJ200, have achieved remarkable performance indicators. The FAR within the combustion chamber has reached approximately 0.030, while the turbine inlet temperature has reached 1850K. Moreover, the combustion efficiency exceeds 0.99, and the temperature distribution coefficient at the outlet is around 0.25. Moving onto the fifth-generation military small-bypass-ratio turbofan engines, such as the F135, the FAR within the combustion chamber has surpassed 0.037, and the turbine inlet temperature stands at approximately 1,980K. Commercial large-bypass-ratio turbofan engines, including the GE90, PW4084, TRENT800 and CFM56, have achieved notable outcomes by adopting technologies such as radial staging, axial staging and rich combustion-quenching-lean combustion (RQL). These advancements have resulted in enhanced reliability, improved economy and better emission characteristics.

In order to achieve superior agility, supersonic cruise, high stealth and other characteristics of fighter jets [1], aeroengines have been subjected to higher design requirements. One of the most urgent requirements is to improve the engine’s thrust-to-weight ratio. Increasing the temperature in front of the turbine ${T_4}$ is the most direct and effective method to enhance the thrust-to-weight ratio of military aeroengines [2]. Furthermore, increasing the design FAR of the main combustion chamber can significantly elevate the temperature in front of the turbine. For example, the U.S. Air Force’s 2030 development plan is projected to have a main combustion chamber FAR of 0.062 [3]. It is evident that the advancement and practical implementation of high FAR (or high temperature rise) main combustion chamber technology are critical in the research and development of advanced military aerospace engines.

An immediate consequence of the high FAR is the corresponding increase in the equivalence ratio within the primary combustion zone of the combustion chamber. To mitigate the occurrence of visible smoke, it is imperative to maintain the average equivalence ratio within the primary combustion zone of the combustion chamber at 1.4 or below. As a result, the volumetric gas flow within the primary combustion zone will also increase accordingly [4]. The augmented air flow within the primary combustion zone not only reduces the equivalence ratio during high operating conditions but also during low operating conditions. This increase in primary combustion zone airflow velocity deteriorates the lean blowoff performance of the combustion chamber [5]. Furthermore, the crucial performance requirements related to lean blowoff, such as ignition during startup and high-altitude re-ignition, both of which are indispensable requirements for reliable operation, cannot be guaranteed [6, 7].

The increased combustion air in the primary combustion zone can be obtained from the air for flame tube cooling and the dilution air [8]. However, reducing the quantity of dilution air would deteriorate the temperature distribution quality at the exit of the combustion chamber. Therefore, the amount of dilution air cannot be reduced. In conclusion, the only way to increase the combustion air quantity is by reducing the cooling air quantity in the flame tube. However, due to the high temperature rise corresponding to the high heat load of the flame tube, the cooling of the flame tube will face significant challenges [9, 10].

The present study introduces the concept of partitioned combustion to address the main issues of high temperature rise (HTR) combustor. To cool the flame tube, an advanced multi-inclined holes cooling method [10] is utilised. The diffuser design adopts the annular curved-wall expanding-angle flow-facing design method [11], which not only minimises the total pressure loss but also facilitates the manufacturing process. The primary and dilution holes are structured in an ‘air intake spoon-bucket-shaped’ configuration [12], which is conducive to increasing the combustion air. Moreover, the design of the flame tube converging section utilises the ‘convex arc flow-facing design’ approach [13], which is applicable for increasing the combustion air. These design choices successfully fulfill the requirements of high FAR and high thermal performance in the combustion chamber.

2.0 Methodology

2.1 Numerical model

2.1.1 Reynolds averaged Navier-Stokes

This paper conducts three-dimensional numerical simulations of the designed HTR combustor using ANSYS Fluent software. Essentially, any fluid flow phenomenon of Newtonian fluids follows the laws determined by the Navier-Stokes (N-S) equations, which are as follows:

Mass conservation equation (continuity equation):

(1) \begin{align}\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \vec v} \right) = 0\end{align}

Momentum conservation equation:

(2) \begin{align}\frac{{D\vec v}}{{Dt}} = {\vec f_b} - \frac{1}{\rho }\nabla p + \frac{\mu }{\rho }{\nabla ^2}\vec v + \frac{1}{3}\frac{\mu }{\rho }\nabla \left( {\nabla \cdot \vec v} \right)\end{align}

where $\frac{{D\vec v}}{{Dt}}$ represents the change in fluid momentum over time, known as the inertial force term; ${\vec f_b}$ represents the body force term; $ - \frac{1}{\rho }\nabla p$ represents the pressure gradient force term; and $\frac{\mu }{\rho }{\nabla ^2}\vec v + \frac{1}{3}\frac{\mu }{\rho }\nabla \left( {\nabla \cdot \vec v} \right)$ represents the viscous force term.

Energy conservation equation:

(3) \begin{align}\rho \frac{{d\left( {\hat u + \frac{{{v^2}}}{2}} \right)}}{{dt}} = \rho {\vec f_b} \cdot \vec v + \nabla \left( {\vec v \cdot \vec \Gamma } \right) + \nabla \left( {\lambda \nabla T} \right) + \rho \dot q\end{align}

where $\hat u$ represents the specific internal energy of the fluid per unit mass, $\rho {\vec f_b} \cdot \vec v$ denotes the work done by body forces on fluid parcels, $\nabla ( {\vec v \cdot \vec \Gamma } )$ signifies the work done by surface forces (pressure and viscous forces) on fluid parcels, $\nabla \left( {\lambda \nabla T} \right)$ indicates the heat received by fluid parcels from the surroundings via thermal conduction, and $\rho \dot q$ denotes the heat received by fluid parcels from the surroundings via radiation.

From Equations (1) to (3), it is evident that the N-S equation set includes four unknowns: $\rho $ , $\vec v$ , $p$ and $T$ . Therefore, an additional relationship is required to complete the solution.

For incompressible flow, where density $\rho $ is a known quantity, there are effectively only three unknowns, corresponding to three equations, which can be solved. For compressible flow, typically dealing with ideal gases, it satisfies the state equation of ideal gases:

(4) \begin{align}p = \rho T\sum\limits_{i = 1}^N {\frac{{{Y_i}{R_u}}}{{{M_i}}}} .\end{align}

The system of equations becomes solvable with the addition of the state Equation (4), resulting in four equations and four unknowns. The flow of Newtonian fluids represents a particular solution of this system of differential equations, and different flows only differ in their initial and boundary conditions. For multi-component flow systems involved in combustion flow, mass conservation is divided into total mass conservation and conservation of each component. The conservation equation for component $i$ is:

(5) \begin{align}\frac{{\partial \left( {\rho {Y_i}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \vec v{Y_i} - \rho {D_i}\nabla {Y_i}} \right) = R{R_i}\end{align}

In theory, given appropriate initial and boundary conditions, one can directly obtain any physical quantity at any position and time in the flow field through numerical solution of the above equations. This includes quantities such as velocity, pressure, temperature and mass fractions of components.

This paper employs the Reynolds averaged Navier-Stokes (RANS) method to time-average the N-S equation system, decomposing the solution variables in the instantaneous (exact) N-S equations into the mean (ensemble-averaged or time-averaged) and fluctuating components. After time-averaging, the N-S equation becomes:

(6) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho {u_i}} \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho {u_i}{u_j}} \right) = - \frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {u\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}} - \frac{2}{3}{\delta _{ij}}\frac{{\partial {u_l}}}{{\partial {x_l}}}} \right)} \right] + \frac{\partial }{{\partial {x_j}}}\left( { - \rho \overline {u'_{\!\!i}u'_{\!\!j}} } \right)\end{align}

called RANS equations. The time-averaged N-S equations introduce the Reynolds stress term $ - \rho \overline {u'_{\!\!i}u'_{\!\!j}} $ to account for turbulence effects compared to the original N-S equations. Under conditions with a sufficiently large Reynolds number and strong turbulence effects in the combustion chamber, the flow becomes highly turbulent. Therefore, neglecting the Reynolds stress term is not permissible in such scenarios. In summary, a crucial task in RANS simulations is to model the Reynolds stress term introduced through averaging the N-S equations. Typically, the Boussinesq’s hypothesis relates Reynolds stresses to mean velocity gradients, transforming the resolution of Reynolds stresses into determining turbulent viscosity ${\mu _t}$ . Furthermore, the turbulent viscosity ${\mu _t}$ can be determined by solving for the turbulence kinetic energy $k$ and the turbulence dissipation rate $\varepsilon $ . Simultaneously, transport equations for $k$ and $\varepsilon $ are solved within the system to close the RANS equations. Specifically, the Boussinesq’s hypothesis is expressed as follows [14]:

(7) \begin{align} - \rho \overline {u'_{\!\!i}u'_{\!\!j}} = {\mu _t}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) - \frac{2}{3}\left( {\rho k + {\mu _t}\frac{{\partial {u_k}}}{{\partial {x_k}}}} \right){\delta _{ij}}\end{align}

This approach offers the advantage of relatively low computational cost in calculating the turbulent viscosity ${\mu _t}$ .

2.1.2 Turbulence model

This study utilised the Realizable k-ε (RKE) two-equation turbulence model in Fluent simulations, which is derived from Boussinesq’s hypothesis and is suitable for strong swirling flows. The RKE model was developed to overcome the limitations of the standard k-ε model in capturing swirling effects. Research findings show that the RKE model outperforms other k-ε models in handling flow separation and complex secondary flows [15]. Moreover, because turbulence models suitable for fully developed turbulence are incompatible with the viscous sublayer dominated by molecular viscous forces within the boundary layer, this study employs the standard wall function approach to link the wall boundary to the fully developed turbulent region.

The modeled transport equations for $k$ and $\varepsilon $ in the RKE model are

(8) \begin{align}\frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + {G_k} + {G_b} - \rho \varepsilon - {Y_M} + {S_k}\end{align}

and

(9) \begin{align}\frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \varepsilon {u_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + \rho {C_1}S\varepsilon - \rho {C_2}\frac{{{\varepsilon ^2}}}{{k + \sqrt {v\varepsilon } }} + {C_{1\varepsilon }}\frac{\varepsilon }{k}{C_{3\varepsilon }}{G_b} + {S_\varepsilon },\end{align}

(10) \begin{align}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\textrm{where}\; {C_1} = \max \left[ {0.43,\frac{\eta }{{\eta + 5}}} \right] \eta = S\frac{k}{\varepsilon } S = \sqrt {2{S_j}{S_j}}\;\; \end{align}

Furthermore, the model constants ${C_{1\varepsilon }}$ , ${C_2}$ , ${\sigma _k}$ and ${\sigma _\varepsilon }$ have been established to ensure that the model performs well for certain canonical flows. The model constants are

${C_{1\varepsilon }}$ =1.44, ${C_2}$ =1.9, ${\sigma _k}$ =1.0 and ${\sigma _\varepsilon }$ =1.2.

In Equations (8) and (9), ${G_k}$ represents the generation of turbulence kinetic energy due to the mean velocity gradients, defined as

(11) \begin{align}{G_k} = - \rho \overline {u'_{\!\!i}u'_{\!\!j}} \frac{{\partial {u_j}}}{{\partial {x_i}}}\end{align}

To evaluate ${G_k}$ in a manner consistent with the Boussinesq hypothesis, ${G_k} = {\mu _t}{S^2}$ , where $S$ is the modulus of the mean rate-of-strain tensor, defined as $S = \sqrt {2{S_j}{S_j}} $ . And ${G_b}$ is the generation of turbulence kinetic energy due to buoyancy, given by

(12) \begin{align}{G_b} = \beta {g_i}\frac{{{\mu _t}}}{{{{\Pr }_t}}}\frac{{\partial T}}{{\partial {x_i}}}\end{align}

For the RKE model, the default value of ${\Pr _t}$ is 0.85. The coefficient of thermal expansion, $\beta $ , is defined as

(13) \begin{align}\beta = - \frac{1}{\rho }{\left( {\frac{{\partial \rho }}{{\partial T}}} \right)_p}\end{align}

${Y_M}$ represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, modeled according to a proposal by Sarkar [16]:

(14) \begin{align}{Y_{\rm{M}}} = 2\rho \varepsilon M_t^2\end{align}

(15) \begin{align}\textrm{where}\,{M_t}\,\textrm{is the turbulent Mach number, defined as}\, {M_t} = \sqrt {\frac{k}{{{a_s}^2}}} \end{align}

2.1.3 Discrete phase model and combustion model

In the designed HTR combustion chamber, liquid fuel (using C₁₂H₂₃ as an alternative fuel to aviation kerosene) is initially atomised by the fuel nozzle and then directly injected into the primary combustion zone. After further fragmentation and evaporation, it becomes gaseous fuel and participates in combustion chemical reactions. The numerical simulation of this process can be accomplished through two models: the discrete phase model (controlling processes such as movement, fragmentation, aggregation and evaporation of the liquid phase represented by droplets) and the combustion model (determining processes such as changes in component mass fractions and heat release caused by chemical reactions).

The discrete phase model (DPM) is utilised in simulating the interaction between the liquid and gas phases using the Euler-Lagrange method. Assuming that the liquid phase ejected from the fuel nozzle consists of discrete liquid droplets, the motion of the droplets is tracked under the conditions of the current gas phase continuous field (Lagrangian method). This allows for further fragmentation, coalescence and evaporation of the droplets, and the impact of this process on the gas phase (source terms in the governing equations) is subsequently taken into account in order to solve the governing equations for the gas phase continuous field (Eulerian method). This model has two important assumptions: firstly, the density of the discrete phase particles is much greater than that of the continuous phase; secondly, the interaction between particles can be neglected. By defining initial parameters such as particle size and velocity distribution at the nozzle exit, the DPM can significantly reduce the computational burden associated with two-phase flow simulations, while effectively capturing the size, distribution and evaporation characteristics of the droplets within the core region downstream of the nozzle. This study utilises the DPM for simulating fuel spray, employing the concept of parcels and the algorithm of O’Rourke [17] for collision and aggregation processes. The fragmentation process incorporates the K-H/R-T (Kelvin-Helmholtz/Rayleigh-Taylor) model [18, 19].

Due to the typical diffusion combustion where the fuel and air enter the reaction zone in different phases, this study adopts a non-premixed combustion model based on the assumption of chemical equilibrium. Upon entry of fuel and air separately into the combustion chamber, it is assumed that their reaction is rapid enough to achieve chemical equilibrium within the solution’s time scale. The mixed gas follows the ideal gas state equation. Hence, the thermal-chemical state at any position during combustion (comprising component mass fraction ${Y_i}$ , temperature $T$ and density $\rho $ ) can be correlated with the local mixture fraction $f$ . In the non-premixed combustion model based on the assumption of chemical equilibrium, direct solution of the chemical reaction process is bypassed. Instead, the thermal-chemical process of fuel combustion is simplified to a mixing problem involving fuel, oxidiser, and intermediate product components. Thus, the energy equation (Equation 3) lacks source terms for chemical reactions, and chemical reactions are controlled by local fuel-rich limits and strain rates. The entire reaction system consists of 11 components, namely C12H23, CH4, CO, CO2, H2O(g), H2O(l), H2, O2, N2, C(s) and OH. Unlike finite-rate methods, the non-premixed combustion model simplifies the solution of the thermal-chemical state by requiring transport equations solely for mixture fractions, rather than for all components. The mixture fraction distribution is obtained by solving the transport equation of the mixture fraction, while the instantaneous mass fractions of the components, temperature, and other parameters are uniquely related to the mixture fraction. Additionally, the non-premixed combustion model also requires solving the transport equation for the mixture fraction variance and linking the instantaneous values of thermochemical parameters with the mean values through the probability density function (PDF) method [20] to introduce the effect of turbulence on chemical reactions. In this paper, the $\beta {\rm{ - }}$ function closest to the PDF observed in experiments has been selected, and the specific shape of the PDF generated by this function is entirely determined by the mean mixture fraction and mixture fraction variance controlled by the aforementioned transport equations.

Because of atomic conservation in chemical reactions, the mixture fraction is expressed in the form of atomic mass fractions. For the combustion reaction involving carbon, hydrogen and oxygen components in aviation kerosene, the mixture fraction is formulated as follows:

(16) \begin{align}f = \frac{{Y{}_C + {Y_H} + {Y_O} - {Y_{O,ox}}}}{{{Y_{C,fuel}} + {Y_{H,fuel}} - {Y_{O,ox}}}},\end{align}

where the subscripts ‘ $ox$ _’ and ‘ $fuel$ ’ denote the inlet positions of the oxidiser stream (i.e., air stream) and the fuel stream, respectively, in terms of mass fractions. In highly turbulent environments, due to the fact that turbulent diffusion is typically much stronger than molecular diffusion, the diffusion coefficients of each component are considered equal. Additionally, for the solution of multiple species and temperature in turbulent flames, the diffusion of components and enthalpy in turbulence can also be assumed to be equal [21], that is, the Lewis number ( $Le$ ) = 1.

The Favre mean (density-averaged) mixture fraction equation is:

(17) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho \bar f} \right) + \nabla \cdot \left( {\rho \vec v\bar f} \right) = \nabla \cdot \left( {\frac{{{\mu _l} + \mu {}_t}}{{{\sigma _t}}}\nabla \bar f} \right) + {S_m} + {S_u}.\end{align}

In addition to solving for the Favre mean mixture fraction, ANSYS Fluent solves a conservation equation for the mixture fraction variance, $\overline {{f^{'2}}} $ [22]:

(18) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho \overline {{f^{'2}}} } \right) + \nabla \cdot \left( {\rho \vec v\overline {{f^{'2}}} } \right) = \nabla \cdot \left( {\frac{{{\mu _l} + {\mu _t}}}{{{\sigma _t}}}\nabla \overline {{f^{'2}}} } \right) + {C_g}{\mu _t} \cdot {\left( {\nabla \overline f } \right)^2} - {C_d}\rho \frac{\varepsilon }{k}\overline {{f^{'2}}} + {S_u}\end{align}

where ${f}' = f - \bar f$ . The default values for the constants ${\sigma _t}$ , ${C_g}$ and ${C_d}$ are 0.85, 2.86 and 2.0, respectively. The mixture fraction variance is used in the closure model describing turbulence-chemistry interactions.

2.1.4 Radiation model and emission model

The discrete ordinates (DO) model is chosen for radiation calculation in this paper, and the radiation model requires solving the radiation transfer equation (RTE). The RTE for an absorbing, emitting and scattering medium at position $\vec r$ in the direction $\vec s$ is:

(19) \begin{align}\frac{{dI\left( {\vec r,\vec s} \right)}}{{ds}} + \left( {a + {\sigma _s}} \right)I\left( {\vec r,\vec s} \right) = a{n^2}\frac{{\sigma {T^4}}}{\pi } + \frac{{{\sigma _s}}}{{4\pi }}\int_0^{4\pi } I \left( {\vec r,{{\vec s}'}} \right)\Phi \left( {\vec s \cdot {{\vec s}'}} \right)d{\Omega '}\end{align}

where $\vec r$ = position vector, $\vec s$ = direction vector, ${\vec s'}$ = scattering direction vector, $s$ = path length, $T$ = local temperature. The DO model considers the RTE in the direction $\vec s$ as a field equation. Thus, Equation (19) is written as

(20) \begin{align}\nabla \cdot \left( {I\left( {\vec r,\vec s} \right)\vec s} \right) + \left( {a + {\sigma _s}} \right)I\left( {\vec r,\vec s} \right) = a{n^2}\frac{{\sigma {T^4}}}{\pi } + \frac{{{\sigma _s}}}{{4\pi }}\int_0^{4\pi } I \left( {\vec r,{{\vec s}'}} \right)\Phi \left( {\vec s \cdot {{\vec s}'}} \right)d{\Omega '}\end{align}

Furthermore, the RTE also accounts for the impact of discrete phase particles on the radiation process.

It is generally recognised that there are three main mechanisms for NOx formation in combustion: thermal, prompt and fuel. Fuel-NOx refers to the oxidation of nitrogen compounds present in the fuel during the combustion process to produce NOx. In this study, aviation kerosene is used as the fuel, with the molecular formula of C₁₂H₂₃, which does not contain nitrogen compounds. Therefore, NOx emissions are modeled using thermal NOx and prompt NOx mechanisms. For thermal and prompt NOx mechanisms, only the NO species transport equation is needed:

(21) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho {Y_{NO}}} \right) + \nabla \cdot \left( {\rho \vec v{Y_{NO}}} \right) = \nabla \cdot \left( {\rho D\nabla {Y_{NO}}} \right) + {S_{NO}}\end{align}

where ${Y_{NO}}$ is mass fractions of NO in the gas phase, and $D$ is the effective diffusion coefficient. The source term ${S_{NO}}$ is to be determined next for different NOx mechanisms.

The formation of thermal NOx is determined by a set of highly temperature-dependent chemical reactions known as the extended Zeldovich mechanism. The principal reactions governing the formation of thermal NOx from molecular nitrogen are as follows:

(22) \begin{align}O + {N_2} \Leftrightarrow N + NO\end{align}

(23) \begin{align}N + {O_2} \Leftrightarrow O + NO\end{align}

(24) \begin{align}N + OH \Leftrightarrow H + NO\end{align}

A third reaction has been shown to contribute to the formation of thermal NOx, particularly at near-stoichiometric conditions and in fuel-rich mixtures.

Under the quasi-steady-state assumption of nitrogen atoms, the net rate of NO generation depends not only on N₂ and O₂, but also on the mass fractions of intermediate products O and OH [23]. This rate is determined by the local temperature and mass fractions of relevant components in the flow field, and is used to compute the source term in the NO transport equation. Solving the NO transport equation subsequently yields the distribution of NO throughout the flow field.

Prompt NOx primarily forms in fuel-rich zones with high hydrocarbon content and low oxygen levels. Specifically, this mechanism entails hydrocarbon ions (CH, CH₂, CH₃ and C₂) generated during combustion colliding with N₂ in the combustion air to produce HCN and CN. These species then react with abundant O and OH in the flame to form NCO, which is further oxidised to NO. Moreover, when HCN concentrations are high in the flame, significant amounts of ammonia compounds (NH_i) are present, reacting rapidly with oxygen atoms to generate NO.

The soot emissions are evaluated using the two-step model, in which the stoichiometry for fuel combustion is set to 3.4012 due to the selection of C₁₂H₂₃ as the fuel. In the one-step Khan and Greeves model [24] ANSYS Fluent solves a single transport equation for the soot mass fraction:

(25) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho {Y_{soot}}} \right) + \nabla \cdot \left( {\rho \vec v{Y_{soot}}} \right) = \nabla \cdot \left( {\frac{{{\mu _t}}}{{{\sigma _{soot}}}}\nabla {Y_{soot}}} \right) + {R_{soot}}\end{align}

where ${Y_{soot}}$ = soot mass fraction, ${\sigma _{soot}}$ = turbulent Prandtl number for soot transport, ${R_{soot}}$ = net rate of soot generation. Furthermore, ${R_{soot}}$ is the balance of soot formation ${R_{soot,\,\,\,\,form}}$ and soot combustion

(26) \begin{align}{R_{soot,\,\,comb}}:{R_{soot}} = {R_{soot,\,\,form}} - {R_{soot,comb}}\end{align}

The rate of soot formation ${R_{soot,\,\,form}}$ is given by a simple empirical rate expression:

(27) \begin{align}{R_{soot,form}} = {C_s}{p_{fuel}}{\varphi ^r}{e^{ - \frac{E}{{RT}}}}\end{align}

where $C{}_s$ = soot formation constant ( $kg/\left( {N \cdot m \cdot s} \right)$ ), ${p_{fuel}}$ = fuel partial pressure ( $Pa$ ), $\varphi $ = equivalence ratio, $r$ = equivalence ratio exponent, $\frac{E}{R}$ = activation temperature ( $K$ ).

The rate of soot combustion ${R_{soot,\,\,comb}}$ is the minimum of two rate expressions [25]:

(28) \begin{align}{R_{soot,\,\,comb}} = \min \left[ {{R_1},{R_2}} \right]\end{align}

The two rates are computed as

(29) \begin{align}{R_1} = A\rho {Y_{soot}}\frac{\varepsilon }{k}\end{align}

and

(30) \begin{align}{R_2} = A\rho \left( {\frac{{{Y_{ox}}}}{{{v_{soot}}}}} \right)\left( {\frac{{{Y_{soot}}{v_{soot}}}}{{{Y_{soot}}{v_{soot}} + {Y_{fuel}}{v_{fuel}}}}} \right)\frac{\varepsilon }{k}\end{align}

where $A$ = constant in the Magnussen model, ${Y_{ox}},{Y_{fuel}}$ = mass fractions of oxidiser and fuel, and ${v_{soot}},{v_{fuel}}$ = mass stoichiometries for soot and fuel combustion. The two-step Tesner model [26] predicts the generation of radical nuclei and then computes the formation of soot on these nuclei. ANSYS Fluent therefore solves transport equations for two scalar quantities: the soot mass fraction (Equation 25) and normalised radical nuclei concentration:

(31) \begin{align}\frac{\partial }{{\partial t}}\left( {\rho b_{nuc}^ * } \right) + \nabla \cdot \left( {\rho \vec vb_{nuc}^ * } \right) = \nabla \cdot \left( {\frac{{{\mu _t}}}{{{\sigma _{nuc}}}}\nabla b_{nuc}^ * } \right) + R_{nuc}^ * \end{align}

In these transport equations, the rates of nuclei and soot generation are the net rates, involving a balance between formation and combustion.

2.2 Structural and aerodynamic design of combustor components

2.2.1 Key parameters of combustor

The HTS triple-swirler main combustor designed in this study differs from conventional combustors by increasing the proportion of combustion air and allowing more air to enter from the dome, with an appropriate increase in dome height. The research is focused on the design of an annular combustor with 20 domes across the annular ring. The design operating condition is set for takeoff. Taking the single-dome combustor as an example, most of the techniques obtained from it can be directly applied to the full annular combustor. Selecting the single-dome combustor in the initial stage of the research design can save costs and reduce the complexity of the work. Based on the requirements of the high FAR combustion organisation in the HTR combustor, a single-dome HTR triple-swirler main combustor is designed. In addition to increasing the outlet temperature of the combustor, it also features a more uniform temperature distribution at the outlet. Key technical parameters such as flow field, temperature field, outlet temperature distribution and exhaust pollution under different operating conditions are obtained. Tables 1 and 2 present the aerodynamic and structural parameters of the designed single-dome annular combustor, while Table 3 outlines the design requirements for overall performance. In this paper, ‘ppm’ denotes measurements expressed in parts per million.

Table 1. Aerodynamic parameters of HTR triple-swirler combustor

Table 2. Structural limit size of HTR triple-swirler combustor

Table 3. Overall performance requirements of HTR combustor

2.2.2 Cooling hole design

The view of the central cross-section of the HTR triple-swirler main combustion chamber designed in this study is illustrated in Fig. 1. The flame tube cooling adopts multi-inclined holes cooling method [10], with a small hole diameter of 0.7mm and a tilt angle of 22°. The hole pattern is in a long rhombus arrangement, with the flow distance between holes, $s$ , being equal to the lateral hole spacing, $p$ , as shown in Fig. 2. In comparison to a regular rhombus arrangement, the long rhombus arrangement results in stronger interference between columns under the same unit area porosity conditions. The upper wall of the flame tube is designed with 44 rows of diverging small holes, with 24 holes evenly distributed in each row along the circumference. The lower wall is designed with 47 rows of diverging small holes, with 16 holes evenly distributed in each row along the circumference. No small holes are designed within 5mm of the main combustion holes and dilution holes. The total number of diverging small holes is 1,653, and the flow coefficient is determined to be 0.78 after multiple iterations.

Figure 1. Model of single-dome HTR triple-swirler combustor.

Figure 2. Schematic diagram of cooling hole arrangement of flame tube.

The dome cooling adopts a direct hole impact cooling method, as shown in Fig. 3, with double rows of concentric circles arranged uniformly around the circumference, each row with 72 holes evenly distributed. The adjacent holes are misaligned by a circumferential angle of 2.5°, with a cooling hole diameter of 0.8626mm. The flow coefficient is determined to be 0.85 after multiple iterations.

Figure 3. Arrangement of dome cooling holes.

2.2.3 Diffuser design

When designing the diffuser of the main combustion chamber, several design methods are commonly used, such as ‘annular straight-wall expansion angle axisymmetry’, ‘annular straight-wall inner expansion angle of zero’, ‘annular straight-wall outer expansion angle of zero’ and ‘dump diffuser’. However, these methods often result in the formation of separated vortices near the wall surface close to the pre-diffuser exit, leading to significant pressure losses. For instance, with the ‘annular straight-wall expansion angle axisymmetry’ method, separated vortices are illustrated in Fig. 4(a). Additionally, considering pressure loss as a key factor, design methods such as ‘constant pressure gradient diffuser’, ‘constant velocity gradient diffuser’, ‘lemniscate diffuser’ and ‘mixed-shape diffuser’ are also commonly employed. Nevertheless, their application is limited due to excessively large exit wall angles and difficulties in manufacturing.

Figure 4. Two methods of diffuser design.

In conclusion, the combustion chamber of this study employs the annular curved-wall expanding-angle flow-facing design method [11] for diffuser design that minimises total pressure losses and facilitates ease of processing, as depicted in Fig. 4(b). The key aspect of this approach lies in the fact that both the inlet and outlet sections of the pre-diffuser are perpendicular to the airflow direction. Additionally, the average radius of the inlet cowl and the centre streamline of the outlet airflow of the pre-diffuser align with the central axis of the combustion chamber, ensuring smooth intake of the cowl. This design method exhibits simplicity and expediency, with the 2D geometric structure composed of straight line segments and circular arcs, facilitating ease of processing. Furthermore, the diffuser designed using this approach eliminates the occurrence of separated vortices, as shown in Fig. 4(a), in the pre-diffuser section, resulting in lower total pressure loss coefficient and higher static pressure recovery coefficient. The average of the total pressures at the outer annular chamber inlet, the inner annular chamber inlet and the cowl inlet is taken as the diffuser exit total pressure, denoted as ${P_{td}}$ . The calculated total pressure loss coefficient of the diffuser is 0.23%. Likewise, the average of the static pressures at the outer annular chamber inlet, the inner annular chamber inlet and the cowl inlet is taken as the diffuser exit static pressure, denoted as ${P_{sd}}$ . The calculated static pressure recovery coefficient of the diffuser is 80%.

2.2.4 Design of primary hole and dilution hole

As the dome height increases with more combustion air, a key issue arises: the transverse jet penetration depth of the primary hole with the conventional simple flat-hole structure reaches only up to the height of the main combustion stage swirler, as illustrated in Fig. 5(a). This significantly weakens the truncation and joint composition on the recirculation zone, resulting in a substantial increase in the axial length of the recirculation zone and a deterioration in the mixing and uniformity of the fuel-air mixture within the primary zone. Consequently, combustion in the primary zone becomes relatively weaker, with more combustion extending beyond the primary hole, ultimately impacting combustion efficiency and the quality of the outlet temperature distribution. Therefore, this study employs the ‘air intake spoon-bucket-shaped structure’ [12] for the primary holes. This structure extends a certain distance into the flame tube, enhancing the penetration depth. It has been verified that the transverse jet penetration depth of the primary hole with this structure can reach the centre axis of the flame tube, as illustrated in Fig. 5(b). Additionally, as the dome height increases, the adjustment capability of the transverse jets formed by the dilution holes with the conventional simple flat-hole structure for the combustion chamber outlet temperature distribution gradually diminishes.

Figure 5. Comparison of penetration depths for different hole structures at the takeoff condition.

In summary, the primary and dilution holes of the combustion chamber designed in this study both use the ‘air intake spoon-bucket-shaped structure’. For the primary hole, a cross-sectional view of this structure is shown in Fig. 6, which is an enlarged view of the red-circled area in Fig. 1. As illustrated in Fig. 6, the design features a ‘short front and long rear’, which increases the proportion of the transverse jet from the primary hole entering the recirculation zone. This configuration facilitates the formation of the recirculation zone in conjunction with the swirling air and reduces obstruction to its development within the primary zone. Additionally, the inner and outer annular chamber flow directions introduced by the primary holes are opposite to the tilt direction of the ‘spoon’ end of the structure. Compared to the ‘straight end’, this design increases the contact area between the inner and outer annular chamber flows and the right long wall surface of the primary hole. This results in increased residence time, improved wall cooling and reduced erosion. Furthermore, it aids in increasing the penetration depth and allows for improved control of the transverse jet direction by adjusting the tilt angle of the ‘spoon’ end.

Figure 6. Primary hole with the air intake scoop-bucket-shaped structure.

The primary and dilution holes are distributed on the upper and lower walls of the flame tube, as depicted in Fig. 7. From Fig. 7, it can be observed that there are a total of four primary holes and four dilution holes in the single-dome combustion chamber. The upper and lower walls have two holes each, with a diameter of 13mm. Moreover, on one side of the flame tube’s wall, the primary holes and dilution holes are arranged in a uniformly intersecting pattern circumferentially. On the same cross-section of the combustion chamber, the primary holes or dilution holes on the upper and lower walls of the flame tube exhibit a uniformly intersecting arrangement circumferentially.

Figure 7. Circumferential distribution of primary holes and dilution holes.

2.2.5 Design of flame tube convergence section

The transition principle from the inlet of the dilution section to the outlet cross-section of the flame tube is to minimise pressure loss, ensure a uniform exit flow field and achieve the desired temperature distribution. When designing the convergent section of the flame tube, the convergence laws typically used for the flame tube’s tail flow area are the Vitosinski law and the constant velocity gradient rule. In cases where the combustion chamber is inclined, it has been observed during the design process that using either the Vitosinski law or the constant velocity gradient rule leads to a sharply converging region at the transition between the straight section and the convergent section of the flame tube due to the increased dome height, as shown in the red-circled area of Fig. 8. This results in a rapid increase in the velocity of the combustion flow at this location, reducing the residence time and adversely affecting the mixing in the dilution zone. In the case of a tilted chamber, Figure 9 compares the upper wall surface inner curves of the convergent section for the flame tube designed using the Vitosinski law” and the constant velocity gradient rule. The comparison shows that the convergent section designed with the constant velocity gradient rule has a greater and steeper curvature on the upper wall of the flame tube. From a mathematical perspective, when designing the upper wall of the flame tube convergent section using the Vitosinski law or the constant velocity gradient rule, an increase in dome height can easily result in inflection points on the upper wall curve, such as point P in Fig. 9 (red-circled area). The inflection point, also known as the point of reverse curvature, refers to the point that changes the upward or downward direction of the curve in mathematics. Intuitively speaking, the inflection point is where a continuous curve transitions from concave to convex. At this point, the upper wall is more susceptible to deformation.

Figure 8. Two traditional design laws for flame tube configurations.

Figure 9. Comparison of two design laws in convergent sections with upper wall surface curves.

In conclusion, the flame tube convergent section of the combustion chamber designed in this study adopts the convex arc flow-facing design ([13]) scheme applicable to the increase of combustion air, as shown in Fig. 10. In the Fig. 10, line D represents the upper wall of the straight section, line A represents the lower wall of the straight section, line G indicates the inlet cross-section of the convergent section, line I represents the outlet cross-section of the flame tube, and line H is the centreline of the combustion chamber, which is the line connecting the midpoint of the diffuser inlet cross-section with the midpoint of the flame tube outlet cross-section. The key to this method is the smooth transition between the straight section and the convergent section of the flame tube wall, preventing abrupt changes in velocity and pressure of the combustion flow at the transition, and ensuring the absence of inflection points. Upon verification [13], it has been demonstrated that the utilisation of the convex arc flow-facing design approach for the main combustion chamber results in higher combustion efficiency, lower emissions of pollutants, a more uniform outlet flow field, improved quality of outlet temperature distribution, as well as easier machining design and better durability.

Figure 10. Convex arc flow-facing design method of flame tube convergence section.

2.2.6 Dome structure and volume distribution

The lean blowoff equivalent ratio ${\Phi _{{\rm{blowoff}}}}$ of swirl cup is generally around 0.45 [27], as defined by:

(32) \begin{align}{\Phi _{{\rm{blowoff}}}} = \frac{{{{\left( {\frac{{\rm{f}}}{{\rm{a}}}} \right)}_{{\rm{blowoff}}}} \cdot 14.7}}{{{W_{{\rm{swirl cup}}}}}},\end{align}

where ${W_{{\rm{swirl cup}}}}$ represents the flow distribution ratio of the swirl cup. For typical military aircraft combustors, the blowoff FAR of ≤0.005 is required at idle condition [6]. Hence, by calculating the Equation (32), it results in ${W_{{\rm{swirl cup}}}}$ ≤16%. It can be observed from Equation (32) that as ${\Phi _{{\rm{blowoff}}}}$ decreases, ${W_{{\rm{swirl cup}}}}$ increases. The combustion chambers for a thrust-to-weight ratio of 8 (FAR of 0.027) and a thrust-to-weight ratio of 10 (FAR of 0.034) are generally designed accordingly. The overall increase in FAR in the combustion chamber has an impact on the conventional combustion chamber dome, as shown in Figs. 11 and 12 [27]. Based on the experimental results of the atomisation using a swirl cup, it can be observed that an air/liquid ratio greater than 3∼4 is sufficient for achieving effective atomization as per the design requirements [28]. From Fig. 12, it is apparent that when the total FAR is ≥0.037, the swirl cup air/liquid ratio is below the lower limit for forming good atomisation, and the atomisation begins to deteriorate, which also has a significant impact on the fuel-air mixture. On the other hand, when the swirl cup’s equivalence ratio is >3, the smoke emission will increase dramatically. Therefore, the conventional combustion chamber swirl cup structure cannot meet the design requirements for high FAR.

Figure 11. Relation between the total FAR of the normal combustor and the equivalence ratio of the dome swirl cup.

Figure 12. Relation between the total FAR of the normal combustor and the air/liquid ratio of the dome swirl cup.

The swirler of the combustor designed in this study is shown in Figs. 13, 14 and 15. The swirler consists of two-stage fuel injection and three-stage axial swirler. The pilot combustion stage fuel nozzle uses a dual-orifice swirl atomiser, with the secondary fuel passage supplying fuel at low load conditions and the primary fuel passage or both passages supplying fuel at high load conditions, with a spray cone angle of 90°. The main combustion stage fuel nozzle uses an air-atomising direct injection nozzle, with 15 evenly distributed fuel injection ports along the circumference, each located between adjacent blades of the main combustion stage swirler, at a spray angle of 90°, effectively utilising multi-point injection and air-atomising nozzles to improve exhaust smoke [29].

Figure 13. Triple-swirler swirler profile.

Figure 14. Schematic diagram of the structure of triple-swirler swirler.

Figure 15. Arrangement of fuel injection orifices of the main stage nozzle.

In the case of low load, only the pilot combustion stage nozzle is fueling, creating a locally rich fuel combustion zone in the primary zone to solve the problem of lean blowoff and to stabilise combustion [30]. In the case of high load, both the pilot combustion stage and the main combustion stage nozzles are fueling simultaneously, with the main combustion stage swirler distributing an increased flow of air, resulting in a comparatively lower FAR in the primary zone, effectively suppressing visible smoke emissions.

The volumetric distribution of the dome air in the HTR combustion chamber with a FAR of 0.037, as designed in this study, is illustrated in Table 4.

The equivalence ratio of the dome is 0.037 ^* 14.7/45.05% = 1.206, indicating that the dome is fuel rich design. Since the equivalence ratio of the dome is less than 1.4, it can effectively addresses the need for suppressing the emission of smoke in the combustion chamber [4]. Due to the unique design of the swirler, assuming that 100% of the air flow of the first stage swirler and 80% of the air flow of the second stage swirler is used to form the traditional swirl cup recirculation zone, and the FAR at idle is 0.0106, then the swirl cup equivalent ratio at idle is 0.0106^*14.7/0.055+0.125^*0.8)=1.0. If the FAR at blowoff is estimated at 50% of the idle condition, then the FAR at blowoff is 0.005, which can meet the current military aircraft requirements for lean blowoff performance. In addition, assuming that this virtual swirl cup air flow atomises the fuel, the dome swirl cup air-liquid ratio is (0.055+0.125^*0.8)/0.037=4.19, satisfying the requirement for a good atomisation air-liquid ratio for the swirl cup.

According to the definition of the swirl number, the calculation formula for the swirl number $SN$ of a bladed axial swirler is as follows:

(33) \begin{align}SN = \frac{2}{3}\left[ {\frac{{1 - {{\left( {\frac{{{R_{\rm{i}}}}}{{{R_{\rm{o}}}}}} \right)}^3}}}{{1 - {{\left( {\frac{{{R_{\rm{i}}}}}{{{R_{\rm{o}}}}}} \right)}^2}}}} \right]{\rm{tan}}\,\theta \end{align}

The intensity of swirl is characterised by the swirl number. When the swirl number is greater than 0.6, it is referred to as strong swirl and a recirculation zone appears. The existence of the recirculation zone is critical to flame stability and fuel-air mixing. When the swirl number is greater than 1.2, it is considered as extremely strong swirl.

For both the inner and outer stages of the pilot combustion stage swirler, the internal swirler airflow is mainly responsible for atomisation, while the external swirler airflow is primarily responsible for creating a recirculation zone. Therefore, it is determined that the swirl number of the external swirler should be greater than that of the internal swirler. As for the main combustion stage swirler, since the main combustion stage airflow needs to mix with a large amount of fuel to achieve premixing and pre-evaporation, a larger swirl intensity is required. However, when the swirl number exceeds 1.2, it becomes an extremely strong swirl, which easily leads to combustion oscillation issues. Based on existing design experience, the preliminary range for the main combustion stage swirl number is determined to be between 0.8 and 1.2.

Table 4. Volume distribution at the head of the combustor

One of the fundamental principles in the design of a swirler is its non-transparency [10], meaning that no light should be visible when observed from the front to the back. This principle can be mathematically expressed by the following formula:

(34) \begin{align}\frac{{\pi D}}{{\rm{n}}} \le L \cdot {\rm{tan}}\theta \end{align}

During the preliminary calculations, the total pressure loss coefficient of the flame tube component of the designed combustion chamber was determined to be 3.5%. The initial values of the flow coefficient for each stage of swirlers were set at 0.84.

The formula for windward area is as follows:

(35) \begin{align}{A_{{\rm{sw}}}} = \pi ({R_{\rm{o}}}^2 - {R_{\rm{i}}}^2){\rm{ - nt(}}{R_{\rm{o}}}{\rm{ - }}{R_{\rm{i}}})/{\rm{cos}}\theta \end{align}

In Equations (34) and (35), $n$ denotes the number of blades; $t$ represents the blade thickness; $L$ indicates the axial length of the blade; $\theta $ represents the blade exit angle; and $D$ denotes the outer diameter of the swirler flow passage.

If there is no extension at the outlet of swirler, the exit geometric flow area can be expressed by the following formula:

(36) \begin{align}{A_{\rm{s}}} = {A_{{\rm{sw}}}} \cdot {\rm{cos}}\theta \end{align}

The formula for mass flow is as follows:

(37) \begin{align}{\rm{m}} = \rho VA = \rho V{C_{\rm{d}}}{A_{\rm{s}}} = \sqrt {2\rho \Delta P} {C_{\rm{d}}}{A_{\rm{s}}}\end{align}

The key structural parameters of each stage swirler, such as the windward area, inner and outer diameters, are calculated based on the above formula. Subsequently, a swirler structural model is established for simulation calculations. The calculation results are then compared with the air distribution ratios of each stage swirler in Table 4. If the results do not match, the initial flow coefficient is adjusted based on the calculation results, and the iteration calculation is performed again until the three-dimensional simulation results match the air distribution results. After multiple iterations, the structural parameters of each stage swirler in the triple-swirler swirler are shown in Table 5.

Table 5. Structural parameters of each swirler

Based on Table 5, it can be observed that the rotation direction for the three stages of swirl flow are respectively counterclockwise, clockwise, and clockwise when viewed from the direction of the inflow to the outflow. Specifically, the counterclockwise rotation is denoted by the minus symbol (-), whereas the clockwise rotation is denoted by the plus symbol (+). The opposite rotation of the inner and outer stages of the pilot combustion stage swirler creates strong shear force, which is beneficial for breaking the fuel film of the pilot combustion stage into fine droplets, thereby enhancing atomisation efficiency. The swirl flow of the main combustion stage is in the same direction as the outer swirl flow of the pilot combustion stage, which helps to enhance the flow field in the primary zone, thereby improving combustion stability.

In the application of combustion chambers, it is common to have an extended section at the exit of the swirlers. In this case, it is necessary to calculate the effective flow area of the component with the extended section, which depends on the following factors [10]: whether the exit area of the extended section is the minimum flow area (if the exit area of the extended section is greater than or equal to the effective flow area, then the impact can be disregarded); the length of the extended section (this factor generally has little impact, mainly due to friction); whether the exit of the extended section is convergent in shape (it should be so).

If the exit area of the extended section is smaller than the effective flow area, and the difference between the effective flow area and the exit area of the extended section is significant, then the following formula should be used for calculation:

(38) \begin{align}\frac{1}{{{A_{{\rm{total}}}}^2}} = \frac{1}{{{A_{{\rm{effective}}}}^2}} + \frac{1}{{{A_{\rm{o}}}^2}},\end{align}

where ${A_{{\rm{total}}}}$ represents the total flow area of the swirler, ${A_{{\rm{effective}}}}$ denotes the effective area of the swirler, and ${A_{\rm{o}}}$ stands for the minimum outlet area of the converging section.

Overall, the total flow area of the swirler ${A_{{\rm{total}}}}$ refers to $A$ in Equation (37). When an extended section is present at the swirler outlet and the outlet is convergent, two scenarios arise: In the first scenario, if the minimum outlet area of the converging section ${A_{\rm{o}}}$ is larger than the effective area ${A_{{\rm{effective}}}}$ (= ${C_d} \cdot {A_s}$ ) of the swirler, then $A$ in Equation (37) corresponds to the effective area of the swirler ${A_{{\rm{effective}}}}$ . In the second scenario, if the minimum outlet area of the extended section ${A_{\rm{o}}}$ is smaller than the effective area of the swirler ${A_{{\rm{effective}}}}$ (= ${C_d} \cdot {A_s}$ ), $A$ in Equation (37) corresponds to ${A_{{\rm{total}}}}$ , which is calculated according to Equation (38). In Equation (38), the effective area of the swirler ${A_{{\rm{effective}}}}$ refers to ${C_d} \cdot {A_s}$ , which is ${C_d} \cdot {A_{sw}} \cdot \cos \theta $ .

According to Fig. 13, it can be observed that there is an extension section at the outlet of the third stage swirler, and the outlet of the extension section is convergent in shape. Therefore, it is necessary to calculate whether the outlet of the extension section is the minimum flow area. Calculations using the data in Table 5 show that the outlet area of the extension section is greater than or equal to the effective flow area, thus the impact can be ignored.

The pilot combustion stage fuel nozzle is used from ignition to the transition stage, after which the main and pilot combustion stage fuel nozzles operate together until reaching the maximum FAR. The formula for calculating the fuel distribution ratio, under takeoff conditions, is as follows [6]:

(39) \begin{align}{\rm{Fuel\,ratio\,of\,the\,pilot\,combustion\,stage\,fuel\,nozzle }} = \frac{{\phi * {f_{idle}}}}{{{f_{SLTO}}}},\end{align}

(40) \begin{align}{\rm{and\,fuel\,ratio\,of\,the\,main\,combustion\,stage\,fuel\,nozzle}} = 1 - \frac{{\phi * {f_{idle}}}}{{{f_{SLTO}}}},\end{align}

where $\phi $ represents the proportion factor; ${f_{idle}}$ and ${f_{SLTO}}$ represent the fuel-to-air ratios under the idle and sea-level take-off (SLTO) conditions, respectively. As the designed HTR combustion chamber lacks turbine cooling airflow, the available air quantity participating in combustion is 100%. Hence, taking $\phi $ =1.0. By utilising Equations (39) and (40), the fuel ratio for the pilot stage and the main stage is determined as 30% and 70%, respectively.

In the staged combustion process, the transition occurs when the overall FAR in the combustion chamber reaches a certain value, at which point the pilot and main combustion stages work together. For HTR combustion chambers, the blowoff FAR is 0.005 [6], which means that the blowoff equivalence ratio is 0.5. Based on the airflow and fuel allocation in the main and pilot combustion stages, the transition point can be set at 0.015, ensuring that the equivalence ratio of the combustion stage after the transition is greater than the minimum stable combustion equivalence ratio, as shown in Fig. 16.

Figure 16. Variation between the equivalence ratio of the primary and main stages with the total FAR.

3.0 Grids division and boundary conditions

3.1 Geometric model and grids division

The structure design and grid generation of the single-dome combustion chamber were conducted using ANSYS preprocessing software NX and ICEM CFD. The computational domain of the combustion chamber includes a swirler, a pre-diffuser, inner and outer bifurcated channels and a flame tube, aiming to achieve the coupled calculation of the entire flow field in the combustion chamber. Due to the adoption of dense multi-inclined holes for the cooling holes in the flame tube, using a structured grid would result in excessive time and effort consumption. Therefore, this study employs an unstructured grid generation. Considering the enormous number of grids, the complexity of the combustion chamber structure, as well as the performance of the working computer, the overall structure of the combustion chamber was divided and treated in zones during the grid generation process for the purpose of ease in searching and error correction, as shown in Figs. 17 and 18.

Figure 17. Partitioning of the HTR combustor structure.

Figure 18. Meshing of HTR combustor.

This article determines the total number of grids to be 11.7 m by comparing the radial distribution of axial velocity at the centre of the reflux zone for different numbers of grids. The conclusion is drawn when the velocity distribution no longer undergoes significant changes with increasing number of grids, suggesting irrelevance of further grid refinement. The experimental results demonstrate the effectiveness of this approach.

3.2 Boundary conditions and solution methods

In numerical calculations, air is treated as an incompressible ideal fluid. The inlet of the combustor adopts the mass-flow-inlet boundary condition, while the outlet of the flame tube is set as the outflow boundary condition. The reference pressure point is located at the centre of the diffuser inlet cross-section. Since this study uses a single dome combustor from a full-annular combustor as the computational domain, both side walls of the entire fluid domain in the single-dome combustor are set with rotational periodic boundary conditions. When solving the energy equation, thermal boundary conditions need to be defined at the wall boundaries. Since the designed combustor does not have turbine cooling air distribution, the two turbine cooling air outlets are set with adiabatic no-slip wall boundary conditions, similar to the casing and diffuser walls. Other walls are set as radiative walls, with an external emissivity of 0.8 and an internal emissivity of 0.7. Considering that the species formation is significantly effected through local temperatures, the external radiative temperature is set to 845K, the design point temperature. The mass flow rates of the intake air for the dome (including various stages of swirlers and dome cooling holes) and the flame tube wall (including primary holes, dilution holes and effusion cooling holes of the flame tube) are obtained through flow coupling calculations.

To perform numerical simulations of discrete phase trajectories and spray combustion, a hollow conical spray model is implemented at both the main and pilot combustion stage fuel nozzles, with specified structural and aerodynamic parameters. The diameter distribution for both the main and pilot stages is set to uniform. The velocity magnitude of the particle streams oriented along the spray cone is set to 80m/s, and the spray cone’s half-angle is set to 45°. The swirl fraction, which determines the relative magnitude of the swirl velocity, is set to 0.5. However, the swirl directions of the main and pilot stages are opposite. The diameter of the particles in the stream for the pilot stage is 0.04mm, whereas for the main stage, it is 0.03mm.

The solver employs the SIMPLE (semi-implicit method for pressure-linked equations) method for pressure-velocity coupling. In the spatial discretisation of the governing equations, the Green-Gauss cell-based method is used for gradient calculation, while the standard method is employed for calculating pressure interpolation. And for all other governing equations, such as those for momentum and energy, apply the second-order upwind differencing scheme.

Gradients are needed not only for constructing values of a scalar at the cell faces, but also for computing secondary diffusion terms and velocity derivatives. The gradient $\nabla \varphi $ of a given variable $\varphi $ is used to discretise the convection and diffusion terms in the flow conservation equations. When the Green-Gauss theorem is used to compute the gradient of the scalar $\varphi $ at the cell centre ${c_0}$ , the following discrete form is written as

(41) \begin{align}{\left( {\nabla \varphi } \right)_{c0}} = \frac{1}{\nu }\sum\limits_f {{{\overline \varphi }_f}} {\overline A _f}\end{align}

where ${\varphi _f}$ is the value of $\varphi $ at the cell face centroid. The summation is over all the faces enclosing the cell. For the Green-Gauss cell-based gradient evaluation, the face value, ${\overline \varphi _f}$ , in Equation (41), is typically taken from the arithmetic average of the values at the neighbouring cell centres, that is,

(42) \begin{align}{\overline \varphi _f} = \frac{{{\varphi _{c0}} + {\varphi _{c1}}}}{2}\end{align}

During computation, the residuals of monitored data and key flow field information vary with iteration steps. The monitored quantities include the average velocity at the centre section of the recirculation zone, the average temperature and maximum temperature at the outlet of the flame tube, as well as the average total pressure. Convergence is achieved when all residuals are below 0.001, and the monitored data exhibit variations within 1% over 100 iterations.

3.3 Simulation accuracy evaluation

Aiming at a HTR main combustion chamber previously studied in our research group [31], the same mathematical model as the combustion chamber designed in this paper was used to carry out the steady simulation calculation of the three-dimensional turbulent two-phase reaction flow with backflow in the whole process. Furthermore, we compared the simulation results of this combustion chamber with corresponding performance test data.

As shown in Fig. 19, a comparison of the simulation results and experimental results for the effect of the Mach number at the combustion chamber inlet on the total pressure recovery coefficient in the cold state reveals that both exhibit similar trends. However, the simulation values are lower than the experimental values, with a difference of approximately 1.4% at the design point.

Figure 19. Effect of inlet Mach number on the total pressure recovery coefficient.

Figure 20. Universal characteristic curve of combustion efficiency.

The measurement of combustion efficiency ${\eta _c}$ for the full annular combustor was conducted using a gas analysis system, and the impact of the universal criterion $\theta $ on ${\eta _c}$ was illustrated in Fig. 20. The definition formula of the universal criterion $\theta $ for combustion efficiency is as follows:

(43) \begin{align}\theta = \frac{{{P_{t3}}^{1.75}{{\rm{e}}^{\frac{{{T_{t3}}}}{{300}}}}}}{{{W_{{\rm{a}}3}}}}\end{align}

Based on Fig. 20, it can be observed that the combustion efficiency ${\eta _c}$ increases with the increase of the universal criterion $\theta $ . When the universal criterion $\theta $ is greater than 0.8, the combustion efficiency ${\eta _c}$ is close to 100%. At the design point of the combustor, the value of $\theta $ is 2.388, and it can be inferred from Fig. 20 that the experimental value of combustion efficiency is nearly 100%, while the simulated calculation shows a combustion efficiency value of 99.80%. This indicates a good agreement between the numerically calculated value and the experimental results of the component.

Figure 21. Effect of severity index S_OTDF on OTDF value.

Figure 21 shows the outlet temperature distribution characteristics from component testing. Here, ${S_{OTDF}}$ represents the severity index of the combustion chamber operating point [31]. The overall temperature distribution factor (OTDF) is one of the indicators of the temperature quality at the combustion chamber outlet, known as the hotspot indicator. It is defined as the ratio of the amount by which the maximum temperature at the combustion chamber outlet, ${T_{4\max }}$ , exceeds the average temperature, ${T_{4{\rm{ave}}}}$ , to the temperature rise across the combustion chamber, ( ${T_{{\rm{4ave}}}}{\rm{ - }}{T_{{\rm{3ave}}}}$ ). Additionally, the OTDF value in Fig. 21 represents the average for the single-dome sector area. The expressions for ${S_{OTDF}}$ and OTDF are as follows:

(44) \begin{align}{S_{OTDF}} = {\left( {\frac{{{T_{t3}}}}{{828.9}}} \right)^{ - 1.1}}{\left( {\frac{{\Delta {T_{4 - 3}}}}{{943.3}}} \right)^{ - 0.3}}{\left( {\frac{{{W_{3c}}\sqrt {{T_{t3}}} }}{{458.3{P_{t3}}}}} \right)^{ - 0.5}}\end{align}

and

(45) \begin{align}OTDF = \frac{{{T_{4{\rm{max}}}} - {T_{4{\rm{ave}}}}}}{{{T_{4{\rm{ave}}}} - {T_{3{\rm{ave}}}}}}\end{align}

where $\nabla {T_{4{\rm{ - }}3}}$ represents the temperature rise across the combustion chamber ( ${T_{{\rm{4ave}}}}{\rm{ - }}{T_{{\rm{3ave}}}}$ ), and ${W_{3c}}$ denotes the available air quantity for the entire annular main combustion chamber. Since the turbine cooling air does not participate in combustion, the available air quantity is defined as the difference between the air flow rate at the combustion chamber inlet and the turbine cooling air flow rate.

At the design point of the combustion chamber, ${S_{OTDF}}$ is 0.97, with the corresponding OTDF experimental value being 0.16, while the OTDF value obtained through simulation is 0.168. It is evident that the OTDF value obtained from numerical calculations is quite close to the OTDF value from component testing.

The comparison of the average radial temperature distribution factor (RTDF) curves obtained from numerical simulation and component testing is shown in Fig. 22. RTDF is defined as the difference between the circumferentially averaged radial temperature distribution value at the combustion chamber outlet, ${T_{{\rm{4avc}}}}$ , and the average gas temperature, ${T_{{\rm{4ave}}}}$ , divided by the difference in average temperatures between the inlet and outlet gases, as follows:

(46) \begin{align}RTDF = \frac{{{T_{4{\rm{avc}}}} - {T_{4{\rm{ave}}}}}}{{{T_{4{\rm{ave}}}} - {T_{3{\rm{ave}}}}}}\end{align}

The RTDF value is taken as the maximum value from the calculated results along the radial height of the RTDF curve. As shown in Fig. 22, the trends of the RTDF curves from simulation and experiment are consistent, though there are significant deviations near the flame tube wall. This discrepancy is due to the influence of side wall radiation and thermal conduction on the thermocouples during the experiment, resulting in lower experimental values.

Figure 22. Comparison of the nondimensional radial temperature distribution curves on the outlet of the combustor.

Overall, the selected mathematical models in this study slightly underestimate the total pressure recovery coefficient in the combustor compared to experimental values. However, they accurately predict the OTDF and combustion efficiency. The model also overestimates the RTDF, but the trend matches experimental values, with a significant deviation near the wall. Comparison between experimental and simulation data indicates a certain level of credibility in the simulation results, and the chosen mathematical models can be effectively used for combustor design and combustion performance prediction. References [32, 33] confirmed that using thermal and prompt NOx emission models is suitable for predicting NOx emissions in the combustion chamber, allowing for comparisons of NOx emission levels under various types of combustion chambers and operating conditions.

4.0 Calculation result analysis

4.1 Total pressure loss coefficient

The cold-state and hot-state total pressure loss coefficients of the HTR combustion chamber under various operating conditions are presented in Table 6.

Table 6. Total pressure loss coefficient for cold and hot fields

According to Table 6, the overall pressure loss coefficient varies with the engine state, indicating that it differs under different operating conditions. As shown in Table 1, the inlet Mach number, $Ma$ , remains constant as the operating conditions change in the combustion chamber, while the Reynolds number, ${\mathop{\rm Re}\nolimits} $ , varies with the changes in operating conditions. According to Ref. [2], an increase in ${\mathop{\rm Re}\nolimits} $ is the primary cause for a decrease in the overall pressure loss coefficient. Due to variations in air flow rate, pressure and temperature, which differ with different engine states, the Reynolds number ${\mathop{\rm Re}\nolimits} $ also varies.

Furthermore, it can be observed from Table 6 that the overall pressure loss coefficient varies with the FAR of combustor, indicating that the total pressure loss in the absence of combustion differs from that with combustion. When FAR=0, indicating a cold state without combustion, the total pressure loss includes the pressure loss through the diffuser and the total pressure loss through the flame tube. When FAR>0, indicating a hot state with combustion, apart from the fact that the overall pressure loss through the diffuser remains relatively constant, the total pressure loss in the hot state through the flame tube is significantly higher than that in the cold state, representing the heating total pressure loss.

According to Table 6, it can be seen that the total pressure loss coefficient of the HTR combustion chamber designed in this paper is less than 6% under the take-off conditions, and less than 7% under the idle conditions, meeting the design requirements of the total pressure loss coefficient for the HTR combustion chamber.

4.2 Central recirculation zone and velocity distribution

The velocity streamline maps at the center sections of the HTR combustion chamber under different operating conditions are shown in Figs. 23 and 24. Figures 25 and 26 illustrate the three-dimensional morphology of the recirculation zone at the take-off condition and the idle condition, respectively, with the left figures representing the front view and the right figures representing the right view observed from the outflow direction. A comparison of the contours of the recirculation zone on the central section of Figs. 25 and 26 is presented in Fig. 27.

Figure 23. Velocity flow diagram of HTR combustor at the take-off condition.

Figure 24. Velocity flow diagram of HTR combustor at the idle condition.

Figure 25. Three-dimensional model of the reflux zone at the take-off condition.

Figure 26. Three-dimensional model of the reflux zone at the idle condition.

Figure 27. Comparison of back-flow region in the center section at different working conditions.

From Figs. 23 and 24, it can be observed that the total air flow of the first stage swirler and the majority of the air flow of the secondary swirler for the pilot combustion stage swirler are utilised to form the traditional swirl cup recirculation zone. The remaining small portion of the air flow of the second stage swirler of the pilot combustion stage swirler and the entire air flow of the main combustion stage swirler flow downstream, adjacent to the recirculation zone, achieving the purpose of flow partitioning and air staging, laying the foundation for the combustion zone-division. Additionally, the designed HTR triple-swirler combustion chamber, whether at take-off or idle operating conditions, forms a stable combustion central recirculation zone.

The 3D model of the recirculation region under both the take-off and idle working conditions, as shown in Figs. 25 and 26, reveals similar shapes and sizes. Overall, the recirculation region takes on a torch-like morphology, with the top resembling a half-closed palm, where four fingers are brought together. This unique shape can be attributed to the uniform intersecting arrangement of transverse jet streams from the primary holes on the inner and outer annular chamber walls of the flame tube, resulting in inconsistent positions of the transverse jet of the primary holes on the inner and outer annular chamber walls where the recirculation zone is intercepted in the circumferential direction, ultimately presenting the three-dimensional model of the recirculation zone as shown in Figs. 25 and 26.

The diagram depicted in Fig. 27 reveals that the recirculation regions for both the take-off and idle conditions are nearly identical in terms of axial length and radial height within the central cross-section of the combustion chamber. However, the recirculation zone for the take-off condition exhibits a more symmetrical and fuller appearance compared to the idle operating condition, which is consistent with the velocity streamline maps shown in Figs. 23 and 24. It should be noted that the two vortices in the downstream of the swirlers in Figs. 23 and 24 are not symmetrical. This is caused by the inconsistent circumferential positions of the transverse jet cutting off the recirculation zone from the primary holes on the inner and outer annular chamber walls.

4.3 Temperature distribution of the central section

The temperature distribution of the central section of the HTR combustion chamber designed in this study under different operating conditions is shown in Figs. 28 and 29.

Figure 28. Temperature distribution contour map of the center section at the take-off condition.

Figure 29. Temperature distribution contour map of the center section at the idle condition.

The diagram in Fig. 28 illustrates that under the take-off condition, the downstream of the swirler forms a stratified combustion flame.

The outer layer of the flame is referred to as the main combustion stage flame, whereas the inner layer of the flame is known as the pilot combustion stage flame. This phenomenon, known as the double combustion zone, is attributed to the flow partitioning and fuel staging. Under high operating conditions, the pilot combustion stage ignites the main combustion stage, facilitating rapid and reliable ignition while also achieving a relatively uniform combustion temperature, as shown in the temperature distribution of the primary zone in Fig. 28. Figure 29 reveals that under low load conditions, the flame exhibits a ‘wing’ shape with a concentrated high-temperature area within the ‘wing’, implying combustion concentration in a small region. At this point, only the pilot fuel nozzle supplies fuel, leading to localised fuel-rich combustion. Such a temperature field structure ensures exceptional high-altitude ignition capability and resistance to lean blowoff during low load operation [34].

In order to address the issue of lean blowoff in the idle condition, the solution lies in deliberately creating an uneven distribution of fuel and air in the combustion zone, which can be referred to as ‘localised enrichment’ [34]. Downstream of the swirler, the flow is divided into two streams: the inner stream, which is used to create the traditional swirl cup recirculation zone, and thus known as the recirculation flow; the outer stream, which is not used in the recirculation zone and instead contributes to the stability of combustion, is referred to as the circumferential flow. Under the idle operating condition, the equivalence ratio of the recirculation flow is biased towards being rich, while the circumferential flow does not have any quenching effect on the recirculation combustion, thereby enhancing the equivalence ratio of specific local regions in the context of significantly lean fuel within the total primary combustion zone, thus meeting the requirement for lean blowoff during idle condition.

4.4 Outlet temperature distribution and combustion efficiency

The temperature distribution at the exit section of the combustion chamber under different operating conditions is shown in Figs. 30 and 31, while the radial temperature distribution at the exit section of the flame tube under different operating conditions is shown in Fig. 32. The various indicators of the outlet temperature distribution and combustion efficiency are presented in Table 7, including average temperature, maximum temperature, OTDF, RTDF and combustion efficiency.

Figure 30. Outlet temperature distribution contour map at the take-off condition.

Figure 31. Outlet temperature distribution contour map at the idle condition.

Figure 32. Comparison of outlet radial temperature distribution at different working conditions.

To characterise the combustion of aviation kerosene, one may employ the following approximate expressions:

(47) \begin{align}{\eta _B} = \frac{{C{O_2} + 0.531CO - 0.319C{H_4} - 0.397{H_2}}}{{C{O_2} + CO + UHC}}\end{align}

where the values of each component are given in volume percentage. UHC refers to the unburned hydrocarbons in the combustion products, excluding CH₄. Since the combustion calculation adopts an 11-component model for non-premixed PDF combustion, the UHC mentioned here represents C₁₂H₂₃.

From Figs. 30 and 31, it can be observed that the temperature distribution profile of the outlet section is similar between the take-off and idle conditions, with the high-temperature region appearing at approximately the same radial and circumferential positions. Figure 32 reveals that the shape of RTDF curve at the outlet section is similar for both take-off and idle conditions, exhibiting a mother-son peak profile characterised by one large and one small peak. However, the mother-son peak profile is more pronounced under the take-off condition. Figures 30 and 31 provide evidence of the presence of high-temperature regions in both the upper and lower halves of the outlet section, where the occurrence of high-temperature regions is more frequent in the upper half. Consequently, this results in the observed distribution trend of the mother-son peak.

According to Table 7, the HTR combustion chamber designed in this paper exhibits a temperature rise level of around 1160K under the take-off condition, with a combustion efficiency greater than 99%. Under the idle condition, the temperature rise level is around 410K, with a combustion efficiency greater than 99.99%. For both take-off and idle conditions, the OTDF ranges from 0.10 to 0.15, and the average RTDF ranges from 0.06 to 0.10, all meeting the design requirements of the HTR combustion chamber.

4.5 Outlet pollutant emissions analysis

The primary pollutants emitted from aviation gas turbines include carbon monoxide (CO), unburned hydrocarbons (UHC), soot and nitrogen oxides (NOx). During engine operation, NOx and soot are mainly generated under high workloads, while CO and UHC are produced in large quantities under low workloads. Therefore, the pollution calculation mainly divides into the emissions of NOx and soot under take-off conditions and the emissions of CO and UHC under idle conditions, as shown in Table 8.

Due to the adoption of an 11-component non-premixed PDF combustion model in the Fluent software, the formula for calculating UHC is as follows:

(48) \begin{align}UHC{\rm{ }}\left( {ppm} \right) = \frac{{{C_{12}}{H_{23}} + C{H_4}}}{{1 - {H_2}O}} \times {10^6}\end{align}

where the values of each component are in terms of mole fractions.

Figure 33 illustrates the spatial distribution of CO mass fraction under the idle condition. Figure 34 displays the spatial distribution of UHC mass fraction at the idle condition. Figure 35 showcases the spatial distribution of NOx mass fraction during the take-off operating condition. Lastly, Fig. 36 presents the spatial distribution of soot mass fraction under the take-off operating condition.

The distribution of CO, as shown in Fig. 33, is predominantly concentrated in the symmetrical upper and lower wing-shaped regions, which correspond to the high-temperature zones generated by localised rich fuel combustion, consistent with the temperature distribution shown in Fig. 29. Due to the local enrichment of fuel in the high-temperature zones, there is insufficient oxygen for combustion, leading to the interruption of the reaction of CO₂ generation from CO and the high-temperature decomposition of CO₂ into CO, resulting in a relatively higher production of CO.

Figure 34 demonstrates that the mass distribution of UHC is similar to that of CO, as the factors influencing the generation of UHC are the same as those influencing CO.

According to Table 8, it can be observed that the molar fractions of CO and UHC in the combustion products at the outlet section of the flame tube under the idle condition are 0×10⁻³² ppm. This magnitude level represents the computational limit of the Fluent software. Thus, it is evident that the designed HTR combustion chamber exhibits extremely low emissions of CO and UHC during the idle condition.

For military engines, stealth capability is of utmost importance as it significantly increases the possibility of infrared detection due to NOx and soot emissions. This article calculates the production of thermal NOx and prompt NOx using the nitrogen oxide generation model provided by Fluent. Thermal NOx is produced in large quantities when the temperature exceeds 1850K, while prompt NOx is produced in significant amounts at lower temperatures and with higher fuel concentrations [35].

As illustrated in Fig. 35, the distribution of NOx under the take-off condition conforms to the temperature field distribution shown in Fig. 28, where the primary mechanism for NOx generation is thermal NOx. The high temperature of up to 2,500K in the primary combustion zone leads to a substantial production of thermal NOx.

By comparing Fig. 35 with Fig. 28, it is evident that the NOx emissions in the combustion zone formed by the circumferential flow gas are lower. Thus, Fig. 35 mainly depicts the distribution of NOx in the combustion zone formed by the recirculation flow gas. This can be attributed to the fact that the circumferential flow gas is located on the periphery of the back-flow region, with a higher gas velocity and shorter residence time, which does not provide conducive conditions for the formation of prompt NOx. Additionally, the circumferential flow combustion zone of the main combustion stage is characterised by premixed combustion, with a lower and uniformly distributed combustion temperatures, thus resulting in a lower quantity of thermal NOx formation.

Table 7. Combustor outlet temperature index and combustion efficiency

Table 8. Outlet pollutant emissions under different operating conditions

Figure 33. Contour map of CO mass fraction distribution at the idle condition.

Figure 34. Contour map of UHC mass fraction distribution at the idle condition.

Figure 35. Contour map of NOx mass fraction distribution at the take-off condition.

Figure 36. Contour map of Soot mass fraction distribution at the take-off condition.

Table 8 reveals that the molar fraction of NOx in the combustion products at the flame tube outlet is calculated to be 3.700ppm, meeting the design requirements for the HTR combustion chamber.

From Fig. 36, it can be observed that the smoke generation is predominantly concentrated in the recirculation zone near the fuel nozzle. This area is characterised by high temperature, rich fuel and low oxygen levels, leading to significant smoke formation. In the dilution zone after the primary holes, due to the supplementary airflow from the dilution holes, the smoke particles further participate in combustion, resulting in the formation of CO or CO₂. As a result, the distribution of soot is essentially absent after the primary holes, as shown in Fig. 36.

According to Table 8, the molar fraction of soot in the combustion products at the flame tube outlet is only 2.55×10⁻⁵ppm, significantly lower than the visible smoke point.

5.0 Conclusion

This study independently presents the design of a novel HTR triple-swirler combustor, introducing innovative design techniques throughout the development of its various components. These include the annular curved-wall expanding-angle flow-facing design method for the diffuser, the air intake scoop bucket type structure for the primary and dilution holes, and the convex arc flow-facing design method for the convergent section of the flame tube. All of these techniques are original to this study. Furthermore, numerical simulations were conducted on the designed HTR combustor, yielding excellent computational results. These findings initially demonstrate that the combustor designed in this study is capable of meeting the performance requirements of a HTR combustor. The simulation results are outlined as follows:

1. (1) The total pressure loss coefficient varies with the engine state and the FAR of combustor. Specifically, the total pressure loss coefficient is lower at the take-off condition compared to the idle condition. Additionally, the total pressure loss coefficient is higher during hot state calculations than during cold state calculations. Furthermore, the designed HTR combustion chamber achieves a total pressure loss coefficient of less than 6% at the take-off condition, and less than 7% at the idle condition, meeting the design requirements for total pressure loss coefficient.

2. (2) Under both take-off and idle conditions, the recirculation zone exhibits a favourable morphology, with well-defined air grading and flow partitioning. The three-dimensional model of the back-flow region remains consistent in shape and size between the two conditions, resembling a torch overall. At the apex of the torch, it takes on the semblance of a semi-clenched palm, where four fingers are brought together. On the central cross-section of the combustion chamber, the back-flow region in the take-off condition appears more abundant and symmetrical compared to the idle condition.

3. (3) In the take-off operating condition, the downstream of the swirler exhibits stratified combustion, where the main combustion stage is ignited by the pilot combustion stage flame, facilitating rapid and reliable ignition of the main combustion stage while achieving relatively uniform combustion temperatures. During the idle operating condition, the high-temperature combustion is concentrated within a small wing-shaped region in the recirculation zone. Furthermore, the primary zone is overall lean and locally rich in fuel, while the peripheral circumferential flow does not quench the high-temperature combustion in the wing-shaped region, ensuring excellent ignition capability at high altitudes and resistance to lean blowoff during low load operations.

4. (4) The temperature distribution of the outlet section is generally similar in take-off and idle operating conditions, with the high-temperature region exhibiting similar circumferential and radial positions. The RTDF curve of the outlet section in both take-off and idle operating conditions shows a mother-child peak contour, comprising of one major and one minor peak; however, the contour is more pronounced in the take-off condition. The HTR combustion chamber proposed in this study exhibits a temperature rise level of approximately 1,160K in the take-off condition and a combustion efficiency of over 99%, whereas it exhibits a temperature rise level of approximately 410K in the idle operating condition and a combustion efficiency of over 99.9%. For both take-off and idle operating conditions, the OTDF is between 0.10 to 0.15, and the RTDF is between 0.06 to 0.10, all satisfying the design requirements of the HTR combustion chamber.

5. (5) Under the idle operating condition, the distribution of CO and UHC is similar, predominantly concentrated in the symmetrical wing-shaped region, consistent with the temperature distribution. Under the computational limit of Fluent software, the molar fractions of CO and UHC in the combustion products at the outlet section of the combustion chamber are both calculated as 0 ×10⁻³² ppm for the idle operating condition. For the take-off operating condition, the NOx generation is lower in the circumferential flow combustion zone while higher in the recirculation flow combustion zone. The smoke generation for the take-off operating condition mainly occurs in the recirculation zone near the fuel nozzle, with the molar fraction of NOx in the combustion products at the outlet section of the combustion chamber calculated at 3,700ppm and the molar fraction of soot calculated at only 2.55×10⁻⁵ppm. The emissions of CO, UHC, smoke and NOx from the HTR combustion chamber designed meet the emission standards.