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 [Reference Hinze14]:

(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 [Reference Mongia and Dodds15]. 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 [Reference Sarkar and Balakrishnan16]:

(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 [Reference O’Rourke17] for collision and aggregation processes. The fragmentation process incorporates the K-H/R-T (Kelvin-Helmholtz/Rayleigh-Taylor) model [Reference Beale and Reitz18, Reference Patterson and Reitz19].

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 [Reference Bazdidi-Tehrani and Zeinivand20] 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 [Reference Masri, Kalt and Barlow21], 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}}} $ [Reference Jones and Whitelaw22]:

(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 [Reference Hanson and Salimian23]. 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 [Reference Khan and Greeves24] 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 [Reference Magnussen and Hjertager25]:

(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 [Reference Tesner, Snegiriova and Knorre26] 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.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 [Reference Lin, Xu and Liu10], 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 [Reference Wang, Li, Lin and Zhang11] 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’ [Reference Wang, Li, Lin and Zhang12] 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 ([Reference Wang, Li, Lin, Zhang and Wang13]) 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 [Reference Wang, Li, Lin, Zhang and Wang13], 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 [Reference Lin, Lin, Zhang and Xu27], 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 [Reference Bahr6]. 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 [Reference Lin, Lin, Zhang and Xu27]. 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 [Reference Lefebvre28]. 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 [Reference Zhang, Zhang, Xu and Lin29].

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 [Reference Liu, Lin, Yuan, Liu, Hu and Xu30]. 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 [Reference Kress, Taylor and Dodds4]. 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 [Reference Lin, Xu and Liu10], 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 [Reference Lin, Xu and Liu10]: 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 [Reference Bahr6]:

(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 [Reference Bahr6], 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.