Predicting Combustor Performance for Hydrogen-Propane Fuel Blends in Gas Turbines: A Coupled Thermofluid and Chemical Reactor Network Model

1. Introduction

Low-emissions combustors, powered by fossil fuels, are a well-developed component of gas turbines. The advent of renewable fuels and a drive towards carbon neutral power generation has revived interest in hydrogen as a fuel [1]. Currently, 98% of the available hydrogen is produced from fossil fuels [2]. To achieve carbon-neutral power production using hydrogen, it must be produced from a carbon free feedstock and renewable energy, commonly termed Green Hydrogen. However, the low availability and high cost of green hydrogen has led to increased interest in hydrogen-hydrocarbon fuel mixtures as a transitional solution, enabling partial decarbonization while maintaining fuel supply security [3]. Hydrogen enrichment in gas turbine combustors is a promising near-term decarbonization strategy. It has been shown to reduce CO and CO₂ emissions [4], improve combustor pattern factors [5], enhance flame stability and widen combustor operating limits [3]. These benefits stem from hydrogen’s higher flame speed and improved mixture reaction rate resulting from the increased oxygen (O) and hydrogen (H) radical pool [6, 7]. While hydrogen enrichment offers clear benefits, it also introduces significant challenges which include a higher likelihood of flashback, increased NOx production and higher flame temperatures. These effects may lead to pollution and increased thermal stresses that require careful consideration. Given the disparate fluid properties and chemical kinetics of hydrogen and mixture fuels compared to pure hydrocarbon fuels, may require that combustors be redeveloped to this end.

Typically, the preliminary aerodynamic design of gas turbine combustors begins with the work of Mongia [8], Melor & Fritskyk [9] and Lefebvre & Ballal [10]. This is an empirical approach using correlations derived from experimental data. Despite its simplicity, a design process using empirical correlations can be implemented to reduce a combustor’s pressure losses and pollutant formation and optimise its mass flow distribution in a timely manner. Moreover, these correlations can also be utilized to investigate off-design performance of the combustor. However, due to the data-driven nature of the correlations, working outside of the intended operating conditions may produce unreliable results and may be a limiting factor of this modelling method.

A brief review of literature pertaining to semi-empirical modelling of gas turbine combustors. Despierre, et al. [11] demonstrated the versatility of the network method by optimising an annular combustor with a genetic algorithm. The use of a genetic algorthm was facilitated by the network’s fast solving time and reached an optimum solution. Gouws [12] developed an empirical and network models to reverse engineer and predict the performance and wall temperatures of a can combustor from a T56 Rolls Royce engine. The empirical model estimated discharge coefficients for combustor flow components, but produced mass flow distribution errors of 11.8-35.7% compared to experimental data. Calibrating these coefficients reduced errors but made components non-generalizable. With a calibrated mass flow distribution, average gas and wall temperatures matched 3-D CFD simulations well, with exit temperatures differing by less than 13%. du Toit & Theron [13] developed a conjugate flow and heat transfer network model of the NASA E3 development combustor. Validation of the network model focused on the results for the combustor wall temperatures where an error range relative to 3-D CFD data was calculated to be between 1 % to 13 %. Mahto & Chakravarthy [14] presented a hybrid approach by coupling a network model and a CFD simulation. The mass flow distribtion through the combustor was determined by the network model and the CFD resolved 3-D turbulence and combustion within the flame tube. This combination significantly reduced the mesh size of the CFD simulation and decreased the time to reach convergence. The researchers then performed a complex multi-objective optimisation focussing on the effect of swirl on the blow-off region, and the improvement of combustion efficiency, total pressure drop and pattern factor. Across these studies, researchers utilized emprical pressure loss correlations for diffusers, holes and swirlers to account for the aerodynamic behavior of the combustor. Similarly the heat transfer rates to the liner and casing is performed using empirical correlations for fluid radiation, convection and surface radiation pathways.

Estimating pollutant emissions during combustor design relies heavily on CFD due to the high cost of experimental combustion characterization tests. Furthermore, the computational resources and time required for transient system-wide CFD modeling is also prohibitive. Consequently, gas turbine designers resort to low-fidelity models that represent combustors as networks of zero- and one-dimensional (1D) chemical reactors, known as chemical reactor networks (CRNs) when performing system-wide simulations. With a reduced requirement for computational resources, implementing CRNs enables the use of complex reaction mechanisms and the timeous exploration of a design space [15].

Several studies, Park et al. [16], Fichet et al. [17]; Colorado & McDonell [18]; Kaluri et al. [19], show that CRNs in the study of gas turbine combustors achieve comparable levels of accuracy when compared to CFD in the prediction of combustion completeness and pollutants. In Fichet et al. [17], NO emissions had an error of 3 %. A wide error range of 0-250 % is shown in Kaluri et al. [19] for manually and automatically extracted networks. Park et al. [16] found that the error of the CRN, relative to experimental data, increased at high combustor loading. For all simulated cases in Park et al. [16], the CRN was able to closely replicate the trends for CO though substantial errors were recorded and estimated NO to within 4 ppm. Grimm [20] automatically extracted reactor networks from CFD data and evaluated the effect of network fidelity on accuracy. It was found that networks with at least four reactors were required to reliably predict NO emissions and increasing the network size did not produce substantial improvements. Importantly, it was shown that accounting for reactor heat loss had a more substantial impact on the accuracy of predicted NO concentrations.

Most CRN studies focus solely on chemical reactor configuration and neglect momentum transport within the combustor. Aerodynamic effects (modelled via flow networks) and combustion (modelled via CRN) are typically treated separately. This separation represents a key limitation because combustors are coupled to turbomachines, where changes in mass flow rate and operating pressure can significantly affect combustor performance and pollutant formation. Many of the previous researchers which investigated the aerodynamic performance of gas turbine combustors using the 1-D network modelling approach did include combustion modelling ( Despierre, et al. [11], Gouws [12] and du Toit & Theron [13]). While the techniques implemented by these authors, chemical equilibrium and constrained adiabatic combustion, produced satisfactory outcomes, they are unsuitable for estimating pollutant emissions. Accurate emissions prediction requires complex techniques such as chemical reactor networks (CRNs) or computational fluid dynamics (CFD) modelling. Limited research literature exists that attempts to integrate both 1D aerodynamic networks and combustion kinetics through CRNs into a unified model capable of capturing sufficient detail for system-wide investigation of fuel mixture effects on gas turbine performance. Haytaysal & Yozgatligil [21] attempted to integrate chemical reactors into a flow network. However, the study focused on the solution methods for network models and provided limited results discussion, only including total pressure changes and outlet combustor temperature predictions that correlated well with benchmark CFD results.

Current literature reveals three key gaps: limited investigation of hydrogen-propane mixtures in gas turbine combustors; 1-D thermofluid network studies lacking generalizability due to oversimplified pressure loss treatment; and CRNs rarely integrated with 1-D flow networks to account for aerodynamic changes, limiting model applicability across combustor operating ranges required for system-wide design and off-design performance simulations.

In the present work, the capability of a 1-D thermofluid network model to predict the performance of a can combustor is assessed for fuel mixtures of propane and hydrogen. The network model is constructed in Flownex SE using flow elements to estimate heat transfer and total pressure changes. Well-established pressure loss correlations specific to combustor components are implemented via custom user scripts. For the first time in Flownex SE, Cantera-based well-stirred reactors (WSRs) and plug flow reactors (PFRs) elements are integrated to account for chemical kinetics within the combustion zones. The performance of the integrated thermofluid network model is validated using both CFD simulation results and experimental data for propane only firing, thereafter both the CFD combustor and thermofluid network models are used to estimate combustor performance when firing various degrees of hydrogen mixing with the propane. The novelty of the present work lies in four key contributions: CFD analysis of hydrogen-propane fuel mixture combustion within a gas turbine can combustor, detailed presentation of network modelling methodology and applicable loss models, integration of Cantera-based reactors within Flownex SE, comprehensive validation for propane-only operation, and extension of this validation to both pure and blended fuel mixtures. The present study focuses on developing a standalone combustor network model within Flownex SE, with the objective of future integration into a full-cycle gas turbine model. The resulting unified gas turbine system model will feature a fully generalizable combustor network capable of predicting engine performance and emissions across a wide range of operating conditions.

2. Materials and Methods

2.1. Case Study Gas Turbine

The analysis of combustor performance is conducted for its operation within a twin-shaft Brayton cycle gas turbine with a layout as shown in Figure 1. This is an experimental system currently at the Mechanical and Mechatronic Engineering Department of Stellenbosch University that was developed by Ssebabi [22]. In Figure 1, the alphanumeric labels demarcate the measurement stations and the control valves.

The gas-generator and power turbine configuration is comprised of two Borg Warner turbochargers, the model K31 and K44, as the high and low pressure machines. At a HP turbine rotational speed of 73000 revolutions per minute (rpm) it operates at a compressor total-to-total pressure ratio of 2.1 and combustor outlet total temperature is 840 °C. With the LP compressor as the load, at the above operating conditions, the machine was estimated to produce 6.5 kW at a thermal efficiency of 3 % [22].The system is designed for propane-only operation, and work is currently underway to evaluate its performance experimentally with hydrogen and propane fuel mixtures. This study forms part of that investigation. The combustor boundary conditions used to evaluate the network model are based on experimental data from Ssebabi [22], presented in Table 1.

2.2. Modelling Theory

2.2.1. Thermofluid Network Modelling

In a system modelled by the finite volume method (FVM) thermofluid network approach, components are represented as 1-D flow elements connected to 0-D nodes to form a network, an example of which is shown below in Figure 2 (left). This discretisation is analogous to the traditional 3-D CFD approach where, the domain is discretised into cells with faces and cell centroids (Figure 2 [right]). In the thermofluid network approach, the mass and energy balance equations are solved at nodes, and the momentum balance equations are solved at flow elements, creating a semi-staggered grid type numerical approach [24]. Closure of the mass, energy and momentum balances are through semi-empirical or empirical component characteristic models for friction (e.g. Darcy-Weisbach correlations), heat transfer (e.g. Dittus-Boelter internal convection correlations) and work (e.g. Euler turbomachine equation). All thermofluid properties such as pressure and temperature are stored at nodes while source terms such as heat transfer rates and pressure losses are stored at elements, along with mass flow rates.

Shown below in equations (1) to (3) are the differential forms for the steady-state 1-D conservation equations for mass, energy, and momentum, respectively.

$$\partial \left(\rho u\right)/\partial x=0$$

(1)

$${\displaystyle \frac{\partial \left(\rho u\left(h+\frac{1}{2}{u}^2+gz\right)\right)}{\partial x}}=q^{‴}-w^{‴}$$

(2)

$${\displaystyle \frac{\partial \left(\rho uu\right)}{\partial x}}=-{\displaystyle \frac{\partial p}{\partial x}}-\rho g{\displaystyle \frac{\partial z}{\partial x}}+{\displaystyle \frac{\partial p_{0,W}}{\partial x}}-{\displaystyle \frac{\partial p_{0,L}}{\partial x}}$$

(3)

In the above equations, $\rho$ [kg/s] is the fluid density, and $x$ [m] is the spatial dimension. On the left-hand-side (LHS) of the energy equation, $u$ [m/s] is the fluid velocity , $h$ [J/kg] is static enthalpy, $g$ [m²/s] and $z$ [m] are the gravitational acceleration constant and the elevation of the system respectively. On the right-hand-side (RHS), the volumetric heat transfer rate to the fluid and the volumetric work done by the fluid interacting with its surroundings is represented by $q^{‴}$ [W/m³] and $w^{‴}$ [W/m³]. Finally, on the RHS of the momentum balance, $p$ [Pa] is the fluid static pressure, $p_{0,W}$ [Pa] is the total pressure rise due to work interactions, and $p_{0,L}$ [Pa] is the total pressure loss across the fluid element due to friction or secondary losses.

Through the volume integration of the differential equations and applying the above-mentioned discretization approaches, the steady-state mass and energy balance equations for fictitious node i are equations (4) and (5). In these equations, ${\sum \dot{m}}_{i\to s}$ [kg/s] represents the mass flows rates from nodes i to the successor (downstream) nodes s, along that edge between these nodes. Similarly, $\sum {\dot{m}}_{p\to i}$ represents the mass flow rates along the edges between the predecessor nodes p and the node i. In the discretized integral-form energy balance equation, $h_{0,i}$ [J/kg] is the stagnation enthalpy of node i and $h_{0,p}$ [J/kg] is the stagnation enthalpy of the predecessor nodes p connected via edge $p\to i$. ${\dot{Q}}_{p\to i}$ [W] is the heat transfer rates calculated along edge $p\to i$ and similarly ${\dot{W}}_{p\to i}$ [W] is the work transfer rate along edges $p\to i$.

$${\sum \dot{m}}_{i\to s}-\sum {\dot{m}}_{p\to i}=0$$

(4)

$$\sum {\dot{m}}_{i\to s}{h}_{0,i}-\sum {\dot{m}}_{p\to i}{h}_{0,p}=\sum {\dot{Q}}_{p\to i}-\sum {\dot{W}}_{p\to i}$$

(5)

To correctly account for the compressible nature of the flow in the combustor, the differential momentum equation should be treated accordingly. As a first step, the derivative on the LHS of equation (3) is expanded using the product rule as shown in equation (6). The second term on the LHS of the below equation is zero for steady-state flow due to the divergence-free condition from the mass conservation equation.

$$\rho u{\displaystyle \frac{\partial u}{\partial x}}+u{\displaystyle \frac{\partial \rho u }{\partial x}}=-{\displaystyle \frac{\partial p}{\partial x}}-\rho g{\displaystyle \frac{\partial z}{\partial x}}+{\displaystyle \frac{\partial p_{0,W}}{\partial x}}-{\displaystyle \frac{\partial p_{0,L}}{\partial x}}$$

(6)

The momentum equation for compressible flow is written in terms of total pressure ($p_0)$ [Pa] and temperature ($T_0$) [K]. To do so, the theorem in equation (7) is substituted into equation (6) [25].

$$\rho u{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial p}{\partial x}}={\displaystyle \frac{p}{p_0}}{\displaystyle \frac{\partial p_0}{\partial x}}+{\displaystyle \frac{p{u}^2}{2T_0}}{\displaystyle \frac{\partial T_0}{\partial x}}$$

(7)

The result is the momentum equation for compressible flow shown in equation (8).

$${\displaystyle \frac{p}{p_0}}{\displaystyle \frac{\partial p_0}{\partial x}}+{\displaystyle \frac{p{u}^2}{2T_0}}{\displaystyle \frac{\partial T_0}{\partial x}}=-\rho g{\displaystyle \frac{\partial z}{\partial x}}+{\displaystyle \frac{\partial p_{0,W}}{\partial x}}-{\displaystyle \frac{\partial p_{0,L}}{\partial x}}$$

(8)

Finally, the discretised differential equation for compressible momentum is written as shown in the equation below, where $j$ and $i$ subscripts indicate the downstream and upstream nodes respectively. $p_0$ is the stagnation nodal pressure, $T_0$ the nodal stagnation temperature, $∆p_{0,W,i\to j}$ work input along edge $i\to j$ and ${∆p}_{0,L,i\to j}$ pressure loss along edge $i\to j$. The remaining overbar properties are the average properties calculated at edge $i\to j$, using a linear averaged approach and the properties at the upstream and downstream nodes.

$${\displaystyle \frac{{\overline{p}}_{i\to j}}{{\overline{p_0}}_{i\to j}}}\left(p_{0,j}-p_{0,i}\right)+{\displaystyle \frac{{\overline{\rho u^2}}_{i\to j}}{2{\overline{T_0}}_{i\to j}}}\left(T_{0,j}-T_{0,i}\right)+{\overline{\rho }}_{i\to j}g\left(z_j-z_i\right)=∆p_{0,W,i\to j}-{∆p}_{0,L,i\to j}$$

(9)

Equations 4, 5 and 9 are the generalised form of the discretised difference equations that are solved in a network modelling approach to determine changes in a system’s pressure, enthalpy, and mass flow rate. These are coupled with the necessary fluid properties and characteristic models for friction, heat transfer and work which will be discussed in the next sections.

2.2.2. Combustor Aerodynamic Loss Models

A can RQL combustor is shown in Figure 3. It is comprised of multiple aerodynamic features, highlighted in the upper portion of the diagram. These are the inlet pre-diffuser (1), snout (2), annular passage (3), axial swirler (4), annular passage (5), primary holes (6), secondary holes (7), flame tube (8) and the outlet diffuser (9).

In the current work, unless stated otherwise, we will use the loss models of Lefebvre [26] and Lefebvre & Ballal [10] to determine the pressure losses through the pre-diffuser, snout, swirler and dilution jets, and the recirculation mass flow rates due to swirl and dilution jets. The changes in pressure for the rest of the combustor can be modelled by existing components within Flownex SE. The total pressure loss through a prediffuser is calculated through equation (10) [27].

$${∆p}_{0,L,PD}={1.75P}_3{R}_a{\left({\displaystyle \frac{m_3\sqrt{T_3}}{P_3}}\right)}^2{\displaystyle \frac{{\tan }^{1.22}\left({\mathsf{\varphi }}_{\mathrm{P}\mathrm{D}}\right)}{A_3^2}} {\left(1-{\displaystyle \frac{A_3}{A_o}}\right)}^2$$

(10)

A similar equation is used, shown below in equation (11), for the pressure loss through the snout:

$${∆p}_{0,L,SN}={1.75P}_3{R}_a{\left({\displaystyle \frac{m_S\sqrt{T_3}}{P_{S_{in}}}}\right)}^2{\displaystyle \frac{{\tan }^{1.22}\left({\mathsf{\varphi }}_{\mathrm{s}\mathrm{n}}\right)}{A_S^2}} {\left(1-{\displaystyle \frac{A_S}{A_{o_s}}}\right)}^2$$

(11)

For a straight bladed axial flow swirler the total pressure drop is calculated using the following equation.

$${∆p}_{0,L,SW}={\dot{m}}_{sw}^2{\displaystyle \frac{K_{sw}\left({\left(\frac{\sec \left(\theta \right)}{A_{SW}}\right)}^2-\frac{1}{A_l^2}\right)}{2\rho }}$$

(12)

In the previous equation, $K_{sw}$ is a constant that differs based on the vane shape (flat or curved) and $\theta$ is the vane angle. For a flat vane, which is the type used in this study, it is 1.3 and 1.15 for curved. The total pressure drop is dependent on the liner cross sectional area, $A_l$ [m²], and the swirler annulus area,$A_{SW}$ [m²]. The swirling flow will cause flow reversal i.e. flow recirculation known as the entrainment mass flow ${\dot{m}}_e$ [kg/s]. As shown in equation 13, from Lefebvre [26], the average entrainment mass flow in the primary zone is dependent on the swirl number ($S_n$ in equation 14) which is a geometric indicator of the swirl intensity and the length of the primary zone, $l_{pz}$ [m].

$$\dot{m_e}=\dot{m_{sw}}\left(0.35+1.4S_N\right){\displaystyle \frac{l_{pz}^2}{2\left(D_{sw}^2-D_{hub}^2\right)}}$$

(13)

$$S_n={\displaystyle \frac{2}{3}}{\displaystyle \frac{1-{\left(\frac{D_{hub}}{D_{SW}}\right)}^3}{1-{\left(\frac{D_{hub}}{D_{SW}}\right)}^2}}\tan \left(\theta \right)$$

(14)

For the current combustor, the recirculation of mass flow rate into the primary and secondary zone may also be induced through opposing jet flow as shown in Figure 4.

The mass flow recirculated by opposing jets, $\dot{m_r}$ [kg/s], is a function of the jet’s angle of entry into the combustion zone, and the temperature ratio between the jet inlet and the respective combustion zone. To estimate the angle of entry, firstly, the jet’s loss coefficient is calculated as in equation 15.

$$C_d={\displaystyle \frac{1.25\left(K-1\right)}{{\left(4K^2-K{\left(2-\alpha \right)}^2\right)}^{0.5}}}$$

(15)

Here, $K$ is the ratio of the jet’s dynamic pressure to the dynamic pressure in the annular passage just upstream of the jet. Similarly, $\alpha$ is a ratio of a single jet’s mass flow rate to the mass flow of the annular passage upstream of the jet. The asymptotic discharge coefficient, $C_{d,n}$, is determined in equation (15) as $K$ approaches infinity. Finally, the rate of recirculated mass flow is calculated as shown in equation (16).

$${\dot{m}}_{r,pz}={0.5\dot{m}}_h\sin \left(\mathrm{asin}\left(\sqrt{{\displaystyle \frac{C_d}{C_{d,n}}}}\right)\right){{\displaystyle \frac{T_{jet}}{T_z}}}^{-0.5}$$

(16)

${\dot{m}}_h$ is the mass flow rate of a single jet and ${\theta }_{jet}$ is the jet’s angle of entry into the combustor’s flame tube. $T_r$ is the temperature ratio between the jet inlet ($T_{jet}$) [K] and the respective combustion zone ($T_z$) [K].

2.2.3. Chemical Reactor Network Models

Chemical reactor models couple the equations describing the chemical, thermodynamic and fluid dynamic processes within a chemically reacting system [28]. There are numerous formulations for chemical reactors, each with assumptions that characterise specific operating conditions. Reactor zones within a conventional combustor, can be considered as constant volume and pressure processes that may be modelled by a WSR and PFR [29, 10, 28]. The conditions assumed within the reactors are that oxidant and fuel are fed constantly and mix rapidly. Under these conditions, the reactions occur uniformly resulting in a uniform reactor temperature and composition. The following discussion will present the equations used in Cantera [30] for the simulation of combustion in a WSR and PFR.

Equation (17) presents the unsteady form of the mass conservation equation written for a single chemical species, $k$, in a WSR. It states that the change in mass of a species depends on the mass generation or consumption of the species by chemical reaction through the term ${\dot{m}}_{k,gen}={V\dot{\omega }}_{k}M{W}_k$, where ${\dot{\omega }}_k$ [mol/m³s] is the net reaction rate per unit volume, $V$ [m³] is the reactor volume and $M{W}_k$ [kg/mol] is the molecular mass of the species. The formulation also considers multiple inlet ($\sum _{in}{\dot{m}}_{in}$) and outlet ($\sum _{out}{\dot{m}}_{out}$) mass flow streams, where $Y_{k,in}$ is the mass fraction of species $k$ at the inlet and $Y_{k,out}$ is the mass fraction of said species leaving the reactor.

The energy conservation equation (18) is written in terms of temperature where $\dot{Q}$ [W] is the net rate of heat addition to the reactor from external sources. Term $\sum _k{h}_k{\dot{m}}_{k,gen }$ [W], describes the combustion heat release of the mixture as a sum of products of specific enthalpy ($h_k)$ and the rate of mass consumption/ generation (${\dot{m}}_{k,gen}$) of the species resulting from the chemical reactions.

$${\displaystyle \frac{d\left(m{Y}_k\right)}{dt}}=\sum _{in}{\dot{m}}_{in}{Y}_{k,in}-\sum _{out}{\dot{m}}_{out}{Y}_{k.out}+{\dot{m}}_{k,gen}$$

(17)

$$m{c}_p{\displaystyle \frac{dT}{dt}}=\dot{Q}-\sum _k{h}_k{\dot{m}}_{k,gen}+\sum _{in}{\dot{m}}_{in}\left(h_{in}-\sum _k{h}_k{Y}_{k,in}\right)$$

(18)

In a PFR, reactants and products flow axially under ideal frictionless conditions. The gas mixture is an ideal gas, and the composition is considered uniform perpendicular to the flow, reducing the problem to a one-dimensional analysis. Although the PFR is a steady-state system, it captures the effect of residence time on chemical kinetics by integrating the species’ net production rates over the flow-through time through the reactor. The governing differential equation for the evolution of a species mass fraction through a 1-D domain, in a PFR, is given in the equation below.

$$\rho u{\displaystyle \frac{d{Y}_k}{dz}}={\dot{\omega }}_k{MW}_k$$

(19)

In the above equation, the change in species mass fraction, $Y_k$, is determined along the flow path in direction $z$ [m] at a velocity $u$ [m/s]. Similarly, the energy equation (20) is formulated to estimate the change in temperature ($T$) [K] where $c_p$ and $\rho$ is the mixture specific heat and density, respectively. The term $-\sum _k{\widehat{h}}_k{\dot{\omega }}_k$, is the heat released due to combustion. ${\widehat{h}}_k$ is the species molar enthalpy and ${\dot{\omega }}_k$ the net reaction rate per unit volume as discussed earlier.

$$\rho u{c}_p{\displaystyle \frac{dT}{dz}}=-\sum _k{\widehat{h}}_k{\dot{\omega }}_k$$

(20)

A generalised and compact form of the net reaction rate is shown below in equation (21). This form accounts for $L$ elementary reactions and $N$ species that contribute to the reaction rate of a single species ($j$). This form is required when implementing large reaction mechanisms like the GRI-Mech.3.0 [31] for hydrocarbons.

$v_{ji}^{\prime }$ and $v_{ji}^{″}$ are matrices of the stoichiometric coefficients for the reactants and products, respectively. The term $\prod _{j=1}^N{\left[X_j\right]}^{v_{ij}}$ is the product of terms of the molar concentration $\left[X_j\right]$ of each species in the reaction mechanism raised to its stoichiometric coefficient $v_{ij}$.

$$\dot{{\omega }_j}=\sum _{i=0}^L{(v}_{ji}^{″}-v_{ji}^{\prime }){ (k}_{fi}\prod _{j=1}^N{\left[X_j\right]}^{v_{ji}^{\prime }}-k_{ri}\prod _{j=1}^N{\left[X_j\right]}^{v_{ji}^{″}})$$

(21)

In equation (21), the rate coefficient $k$ is in the modified Arrhenius form, equation (22), where the subscripts $f$ and $r$ is for the forward and reverse reaction, respectively.

$$k=A{T}^b\exp \left(-{\displaystyle \frac{Ea}{R_{u}T}}\right)$$

(22)

The magnitude of the rate coefficient is determined by the pre-exponential factor ($A$), static temperature of the reactants ($T$) and the activation energy ($Ea$). The temperature exponent ($b$), pre-exponential factor and activation energy are empirical constants specific to the reaction mechanism, and $R_u$ is the universal gas constant.

2.2.4. Heat Transfer Modelling

In this study, heat transfer is considered within the flame-tube, annular passage of the casing, and to the environment only.

Convective heat transfer between the fluid and solid in the combustor is modeled using the Dittus and Boelter heat transfer coefficient calculated in equation (23). Here, $k$ [W/m²K] is the fluid conductivity, $D_h$ [m] is the hydraulic diameter, $Re$ is the Reynolds number and $Pr$ is the Prandtl number. The constant $n$ is 0.4 if the fluid is being heated and 0.3 if the fluid is being cooled. Despite being formulated for inner tube heat transfer, equation (23) can be applied to an annular passage by using the annulus’ hydraulic diameter (D_out - D_in) [32].

$$h_c=0.023\left({\displaystyle \frac{k}{D_h}}\right)R{e}^{0.8}{Pr}^n$$

(23)

Fourier’s Law determines the rate of heat transfer through solid objects while equation (24) is used to model the radiation heat transfer rate. $\epsilon$ is the surface emissivity, $\sigma$ [W/(m²·K⁴)] is the Stefan-Boltzmann constant, T [K] is the temperature of radiating ($i$) and receiving ($j$) surfaces.

$${\dot{Q}}_{rad}=\epsilon \sigma A\left(T_i^4-T_j^4\right)$$

(24)

The emissivity of luminous flames (${\epsilon }_f$) from Lefebvre and Ballal [10] is utilized to estimate gas-to-solid radiation as shown in equation (25). This expression considers the effects of fluid pressure ($p$), and air ($\dot{m_a}$) and fuel ($\dot{m_f}$) mass flow rates on radiation. It accounts for effect of water vapour on luminosity ($L={\displaystyle \frac{336}{H^2}}$) by including the fuel hydrogen mass fraction ($H)$. The mean beam length in the combustion zone is calculated as $l_b={\displaystyle \frac{3.4V}{A_s}}$ where $V$ [m³] is the reactor volume and the $A_s$ [m²] is the area of the reactor. $T$ [K] is the temperature of the fluid ($f$) and the wall ($w$)

$${\epsilon }_f=1-\exp \left(-290P{\displaystyle \frac{336}{H^2}} {\left({\displaystyle \frac{\dot{m_f}}{\dot{m_a}}}{\displaystyle \frac{3.4V}{A_s}}\right)}^{0.5} T_f^{-1.5}\right)$$

(25)

Fluid absorptivity (${\alpha }_f$) is related to the emissivity through the temperature ratio between the fluid and the wall, as shown in Equation (26).

$${\displaystyle \frac{{\alpha }_f}{{\epsilon }_f}}={\left({\displaystyle \frac{T_f}{T_w}}\right)}^{1.5}$$

(26)

2.2.5. CFD Modelling

To simulate steady-state turbulent flow in a combustor, the Reynolds Averaged Navier-Stokes (RANS) equations are utilised. In numerical studies of similar combustors, RANS simulations using the Realizable k-ε turbulence closure is often employed to predict hydrogen and mixed fuel combustion accurately when compared to experimental results [33, 34, 35]. The Realizable k-ε formulation is more computationally economical than the k-ω Shear Stress Transport (SST) model and comparably capable of modelling flows with adverse pressure gradients, separation and reverse flow caused by swirl as seen in a combustor [36, 37]. Therefore, in the present work, turbulence closure is accomplished using the Realizable k-ε turbulence model with non-equilibrium wall treatment.

The RANS continuity, momentum and energy balance equations solved by Ansys Fluent version 24R2 are shown in equations (27), (28) and (30) respectively.

$${\displaystyle \frac{\partial }{\partial x_i}}\left(\rho u_i\right)=0$$

(27)

$${\displaystyle \frac{\partial }{\partial x_j}}\left(\rho u_i{u}_j\right)=-{\displaystyle \frac{\partial p}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left[\mu \left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}-{\displaystyle \frac{2}{3}}{\delta }_{ij}{\displaystyle \frac{\partial u_i}{\partial x_i}}\right)\right]+{\displaystyle \frac{\partial }{\partial x_j}}\left(-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$$

(28)

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\left(\rho k+{\mu }_t{\displaystyle \frac{\partial u_k}{\partial x_k}} \right){\delta }_{ij}$$

(29)

In the above equations, $\rho$ is the mean density of the fluid which is, for an ideal gas mixture, based on the mass fraction weighted sum of the species molecular weight. $u_i$ and $x_i$ are the velocity and spatial dimension in direction $i$, respectively. In the momentum equation (28), $p$ is the static pressure, $k$ is the turbulence kinetic energy and $-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ is the Reynolds stress term which, through the Boussinesq hypothesis, is modelled using equation (29). The turbulent viscosity, ${\mu }_t$, introduced by equation (29) is calculated as ${\mu }_t=\rho C_{\mu }{\displaystyle \frac{k^2}{\epsilon }}$ where $\epsilon$ is the turbulence dissipation rate. In the Realizable formulation, $C_{\mu }$ is a function of the mean strain and rotation rates, the angular velocity of the system rotation, and $k$ and $\epsilon$ [38].

When using the partially premixed combustion model in Ansys Fluent (discussed later), the total enthalpy ($H$) formulation of the energy equation is implemented as shown in equation (30) [38]. Total enthalpy includes the species enthalpy of formation and therefore accounts for the heat released during combustion. As such, the heat released due to combustion is not included in the source term, $S_h$, of equation (30).

$${\displaystyle \frac{\partial }{\partial x_i}}\left[\rho u_{i}H\right]={\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{{ k}_{eff}}{c_p}}{\displaystyle \frac{\partial H}{\partial x_j}}\right)+S_h$$

(30)

Turbulent viscosity is used to determine the effective thermal conductivity ($k_{eff}$) of the fluid seen in equation (30), $k_{eff}=k+{\displaystyle \frac{c_p{\mu }_t}{{Pr}_t}}$. Here, $k$ is the fluid thermal conductivity and the turbulent Prandtl number, ${Pr}_t$, is a constant value of 0.85.

The heat transfer interaction between fluid and solid material is considered in this study. A Conjugate Heat Transfer approach is implemented where the fluid and solid media form a single and continuous mesh. Equation (31) is the energy balance equation for the solid media only. It accounts for convective heat transfer as a result of moving walls (first term on RHS), conduction (first term on LHS) and a volumetric heat source term. $\overrightarrow{v}$ is the velocity vector of the wall, $h$ is the sensible enthalpy, $k$ is the conductivity of the solid and $S_h$ is the volumetric heat source term.

$$\nabla \ast \left(\overrightarrow{v}\rho h\right)=\nabla \ast \left(k\nabla T\right)+S_h$$

(31)

Convective heat transfer between the solid and fluid domain is modelled using the law-of-the-wall for non-dimensional temperature as seen in equation (32) [38]. Here, on the LHS, $T_w$ is the wall temperature, $T_p$ is the temperature of the wall adjacent cell-centroid, $u^{*}$ is calculated as $u^{*}=C_{\mu }^{{\displaystyle \frac{1}{4}}}{k}_p^{{\displaystyle \frac{1}{2}}}$. The first term on the RHS $(T_c^{*}$) determines the convective-conductive heat transfer through equation (33). The equations available to calculate $T_c^{*}$ are implemented based on the height of the viscous ($y^{*}$) and the thermal sublayer ($y_t^{*}$) where$E$ is an empirical and $K$is the von Kármán constant. $P$ is a function that accounts for the influence of roughness and the Prandtl number on heat transfer in the laminar sublayer.

The second and third term on the RHS of equation (32) are included to account for compressibility effects when using pressure-based solvers. $k_p$ and $U_p$ are the turbulent kinetic energy and mean velocity at the wall adjacent cell-centroid.

$$T^{*}={\displaystyle \frac{\left(T_w-T_p\right)\rho c_p{u}^{*}}{\dot{q}}}=T_c^{*}+{\displaystyle \frac{1}{2}}\rho Pr{\displaystyle \frac{C_{\mu }^{\frac{1}{4}}{k}_p^{\frac{1}{2}}}{\dot{q}}}{U}_p^2+{\displaystyle \frac{1}{2}}\rho {\displaystyle \frac{C_{\mu }^{\frac{1}{4}}{k}_p^{\frac{1}{2}}}{\dot{q}}}\left({Pr}_t{U}_p^2+\left(Pr-{Pr}_t\right)U_c^2 \right)$$

(32)

$$T_c^{*}=\left\{\begin{array}{c}\mathrm{Pr}{y}^{*}, y^{*}<y_T^{*}\\ {Pr}_t\left[{\displaystyle \frac{1}{K}}\ln \left(E{y}^{*}\right)+P\right], y^{*}>y_T^{*}\end{array}\right.$$

(33)

Radiation is modelled using the discrete ordinance (DO) model, with four angular divisions and three pixels in both the azimuthal (phi) and polar (theta) directions. The fluid absorptivity is determined through the weighted-sum-of-grey-gasses method (WSGGM), only considering the tri-atomic molecules of CO₂ and H₂O. Radiation scattering is negligible and, therefore, neglected since the radiating medium is gaseous.

The evolution of the reactant and product species through the domain is modelled in the present work using the partially-premixed combustion approach. Turbulence-chemistry interaction is modelled using the non-adiabatic FGM model coupled with the full GRI-Mech 3.0 reaction mechanism [31] for hydrocarbons. The Zimont turbulent flame speed model is implemented for premixed combustion, and the progress variable variance is treated as a transported scalar. To account for mixing and chemical timescales, the turbulent flame speed setting is selected for turbulence-chemistry interactions. Laminar flame speed is estimated using the Metgalchi-Keck law [39] for propane, and the precomputed PDF values are used for hydrogen and fuel mixtures. In ANSYS Fluent only the mass fractions of CO and CO₂ are considered for hydrocarbon combustion with a weight of unity, ensuring that the progress variable increases through the flame [38]. Since CO exists in small quantities relative to CO₂ within a combustor, the mass fraction of CO₂ will bias the progress variable towards 1. This may cause non-physical behaviour since the FGM model will regard a mixture as burnt when the progress variable is close to or equal to 1 (high CO₂ concentrations), causing reaction progress to slow or stop. Therefore, in the standard FGM model, when combustion is only mapped using the progress variable, the estimation of slower evolving species like NO and the post flame burnout of CO is inaccurate [40, 41]. To remedy this, as suggested by Klarmann, et al. [40] and Efimov, et al. [41], transport equations for CO and NO mass fractions are solved and their net reaction rates are corrected based on the transported values.

2.3. Model Development

In the following section the above thermofluid network theories (flow equations, pressure losses, CRNs and heat transfer) will be applied to model the case-study gas turbine combustor. Furthermore, the 3-D CFD model setup of the combustor will also be discussed. Firstly, the three experimental running points simulated will be discussed along with the various propane-hydrogen mixture simulation boundary conditions.

2.3.1. Simulation Cases

Initial validation of the CFD and network models will be performed for propane-only operation. To this end, the experimental results of [22] are used as simulation boundary conditions. 3 engine running points have been chosen across the systems operating curve resulting in 3 simulation cases. Running point MGT3 is considered the primary operating point as it is closest to the design point of the system.

Once the CFD and thermofluid network models are validated using the propane-only experimental data, the effect of hydrogen blending is investigated. Since there is no experimental data available for the hydrogen-blending operating points, the CFD simulations will be used as the ground truth when investigating the network model accuracy. For the three cases (MGT1-MGT3), different hydrogen blending ratios are simulated. The pure hydrogen mass flow rates were selected to achieve outlet combustor temperatures similar to that of pure propane combustion for the three running points. Table 2 below shows the mass flow rates for pure propane and pure hydrogen operation.

The blending ratios were defined based on the nominal fuel mass flow rates for the pure fuel conditions. For example, for the 75-25 ratio (propane-hydrogen), the fuel flow rate for MGT3 is 75% of the nominal mass rate flow of propane (0.0045 kg/s) and 25% of the nominal hydrogen mass flow rate (0.0018 kg/s). For input simplicity in CFD, the total fuel mass flow rate is 0.0033 kg/s with a 0.883 mass fraction of propane and 0.117 mass fraction hydrogen.

Table 3 summarizes the mass flow rates and corresponding species mass fractions for all mixture operating points. This method of blending fuels ensures matched outlet temperatures across cases. Doing so will allow a direct comparison of the CFD data for all fuel compositions to determine the effect of hydrogen addition on combustor performance.

2.3.2. 3-D CFD Model Setup

The domain used in the 3-D CFD simulation is shown in Figure 5. As illustrated, the inlet boundary conditions are the air and fuel mass flow rates and the static temperatures of each stream. The compositions of the fuel and air are specified in the setup of the FGM model as mass fractions of the species in the fuel and oxidiser streams. Therefore, for each boundary condition, a mixture fraction is specified which is 1 for the fuel and 0 for the air (oxidiser). A static pressure outlet boundary condition is specified at the outlet of the combustor.

Figure 5 (bottom) shows an isometric view of the combustor highlighting the segmented volumes corresponding to the combustion zones and outlet transition duct (nozzle). These zones are present in the generated mesh and thus available in Ansys Fluent for post-processing purposes. The segmentation enables the extraction of mass-weighted average quantities from each zone for direct comparison with the network model.

The FGM model is implemented using the nonadiabatic diffusion flamelet formulation with full GRI-Mech 3.0 reaction mechanism and the accompanying thermodynamic and transport properties. The fuel and oxidiser were specified as mass fractions corresponding to each running point and fuel mixture in Table 2 and Table 3 with the respective inlet temperatures. The outlet static pressure of the combustor is operating pressure for the generation of the flamelet and PDF table. This pressure corresponds to the operating value specified in the operating conditions of the simulation and therefore the combustor outlet gauge pressure is 0 Pa for all cases.

As suggested in ANSYS [38], for large chemical mechanisms, the absolute error tolerance in the Flamelet generation process was reduced to 1e-20 and the relative error tolerance to 1e-05. The Automated Grid refinement of the Flamelet table was implemented with a maximum change in value ratio of 0.35. The initial scalar dissipation value was reduced to 0.001. The maximum number of flamelets was increased to 400, and the number of grid points in the table was increased to 64 in the mixture and reaction space which improved simulation stability. In the PDF table parameters, the minimum temperature of 260 K was set well below the lowest temperature provided in the boundary conditions, and the number of species in the PDF table is 20.

To model heat transfer, the solid components of the combustor seen in the isometric view of Figure 5 were included in the mesh for a CHT simulation. The solid-fluid wall interfaces of the mesh are therefore coupled to the energy solution, and the exterior wall of the casing is allowed to radiate to the surroundings. Radiation within the combustor is modeled using the DO method with four angular divisions and three pixels in both the azimuthal (phi) and polar (theta) directions. The absorption coefficient is determined using the WSSG method where the beam length is based on the domain.

Pressure and velocity are coupled, with spatial gradients calculated using the cell-based least squares method. Pressure is discretized using the PRESTO! approach, and all equations are discretized using a second-order upwind formulation, except radiation, which uses first-order upwind. The pseudo-time method is implemented with a global time step to enhance simulation stability.

2.3.3. Mesh Independence Study

The 3-D domain was discretised into a poly-hexcore mesh. Care was taken to ensure that the y-plus values were between 30 and 300 for the non-equilibrium wall treatment. Two mesh refinements were performed and evaluated for a grid-independent solution. In Figure 6, the static temperature and velocity magnitude along the axis of the combustor for the coarse, medium and fine mesh are plotted for MGT3. The results show that minor differences in temperature and velocity predictions are observed moving from the medium to fine mesh.

To further demonstrate that the medium mesh is appropriate and mesh independent, the GCI values for outlet total temperature and inlet total pressure were estimated. The results are shown in Table 4 below.

For both temperature and pressure, low GCI values are estimated with a range of 0.00% to 0.12% indicating a high degree of numerical precision for the medium and fine mesh. The refinement ratio is 1.3 or greater, indicating satisfactory refinement [42]. In addition, despite having significantly less cells, the results of the medium mesh deviate only slightly from the fine mesh. Consequently, the medium mesh was selected for further use in this study due to the lower computational cost and levels of precision and accuracy comparable to the fine mesh.

2.3.4. Thermofluid Network Model Setup

The 1-D network model, shown in Figure 7, has been constructed in Flownex SE. The model was developed with an emphasis on simplicity, while being informed by 3-D CFD results. As such, the network model closely represents the layout introduced by Bragg [43], where the primary combustion zone is assumed to be a WSR and post-flame burnout is modelled using PFRs.

Standard flow elements are used to model heat transfer and all flow passages. Jets (holes) are modelled as restrictors with discharge coefficients; other sub-components are modelled using pipe elements with loss coefficients for frictional and obstruction losses. The recirculation effects within the combustor, discussed later, are approximated by redirecting mass flow through flow resistance elements with prescribed mass flow rates. Discharge and loss coefficients, as well as recirculation mass flow rates, are determined by the correlations in Section 2.2.2. The calculations are performed by embedded MathCad user-defined scripts, which are executed after each iteration and update the corresponding flow element parameters until convergence is achieved.

Flownex SE can execute Python scripts within its C# scripting components through PythonNET. This facilitates the use of Cantera WSR and PFR in Flownex SE. In this network model, a C# script is executed within a compound component. It obtains the reactor mass flow rate, inlet temperature, pressure, and list of reactant mass fractions and sends it to the Python script that executes the Cantera reactor calculations. An exit temperature and a list of product mass fractions is then returned to Flownex SE and written to the outlet of the pipe element used to approximate the combustion zone in the flow network.

The results of the 3-D simulation were used to inform the shape of the network. Shown in Figure 8 are streamline traces initiated from the PZ jets (left) and the SZ jets (right). It is evident that there is recirculation into the PZ from the flow entering through the PZ jets. This flow phenomenon is not present in the SZ. The magnitude of the reculation mass flow is estimated using equation 15 from Lefebvre [26] instead of being informed by the CFD. This is, of course, a crude estimation as the correlation does not account for swirl flow within the flame tube and the recirculation zone upstream of the opposed PZ jets. Furthermore, consistent with the mechanism described in Lefebvre [26], recirculation flow induced by swirl, is approximated in the network model. This is achieved by returning products at the outlet of the PZ to its inlet at a rate determined by equation 15. Mass flow entrainment due to swirl is critial for flame holding in this combustor as entrained hot reaction products mix with incoming unburnt reactants in the primary zone to sustain combustion.

The type of reactor used to model each reaction zone was informed by the average Damköhler number of each combustion zone in the CFD results. The Damköhler number was calculated as a ratio of the combustion zones average Kolmogorov time scale to the rate coefficient of the slowest reaction estimated at the average combustion zone temperature [44]. In Figure 9 are three plots for each running point. The plots mark the Damköhler number of each combustion zone for three trends that represent three fuel compositions, namely: propane-only (100-0), hydrogen-only (0-100) and a 50 % share of each fuel (50-50).

Notably, the Damköhler number is much less than 1 for all zones in each data set. This indicates that each zone is mixed sufficiently such that the reaction rates will be chemically limited. Therefore, for the combustor in this study, the selection of WSRs to approximate the PZ and SZ is justified. The DZ also exhibits chemically limited combustion, and it is modelled as a PFR to account for the effect of an increased residence time and elevated reaction temperature on reaction kinetics.

3. Results and Discussion

This section presents the results of the 3-D numerical and network model. First, 3-D CFD results for combustor operation using propane, hydrogen and propane-hydrogen fuel mixtures are presented. The intent is to develop an understanding for how hydrogen supplementation may augment combustor operation. Secondly, a detailed validation of the network model is carried out for propane only operation. Finally, the network model is validated for operation on fuel mixtures of hydrogen and propane.

3.1. 3-D Numerical Model

The experimental data available to validate the CFD model is limited to the combustor outlet total temperature at various running points as recorded in Ssebabi (2020) for propane-only operation (see Table 1). Although the validation dataset is limited, an agreement with the outlet total temperature is considered sufficient for this study as the CFD results will be used to validate a system-level thermofluid network model. A detailed validation of the internal flow structure and flame dynamics is beyond the scope of this work since the network model’s predictive accuracy at the outlet of the combustor is of primary importance.

The CFD and experimental results for the combustor outlet total temperature are presented in Table 5. The CFD model is shown to predict a combustor outlet total temperature 64 K to 85 K lower than the experimental results of Ssebabi [22]. This represents a 5.9 % (MGT1) to 7.6 % (MGT2) difference which is within the typical error range of approximately 5 % and 7 % found in the CFD validation studies of combustors in Cam, et al. [45] and Benaissa, et al. [46], respectively. The close agreement in the outlet total temperature suggests that the combustion and heat transfer processes are captured reasonably well by the numerical model for its intended application.

Figure 10 is a compound figure that displays the change in key combustor parameters with the increase in hydrogen fuel content for each running point. Figure 10-1a is the PZ equivalence ratio and 1b is the global equivalence ratio. Figure 10-2a and 2b are plots for the change in total inlet air and fuel pressure pressure, respectively, normalised to the data for propane only. Figure 10-2c shows the combustor efficiency as a function of blending ratio.

Seeing as hydrogen has a higher mass energy density compared to propane, it is observed that the equivivalence ratios reduces as the blending ratio is increased. Despite common operating conditions for each fuel mixture, the change in total pressure across the combustor varies with changes in fuel blending ratios as seen in Figure 10-2a. The trends for MGT2 and MGT3 are similar, initially decreasing, then increasing from blending ratio 50-50 to pure hydrogen operation. When compared to propane only, the change in total pressure across the combustor is at most 3.5 % and 2.1 % higher for MGT2 and MGT3, respectively. The change in total pressure for MGT1 is lower for all blending ratios. Expectedly, and due to hydrogens lower density, increasing its share in the fuel mixture increases the fuels total inlet pressure (Figure 10-2b). Noteably, the maximum inlet total fuel pressure is for the 25-75 mixture for all running points, and not for hydrogen only operation.

Data for combustion efficiency presented in Figure 10-2c is estimated by comparing the total outlet enthalpy of the CFD simulations to the values estimated by a well-stirred reactor for the same boundary conditions. Operating point MGT3 is closest to the design point and has the highest combustion effiency of 89 % for propane-only operation. It is followed by MGT2 and MGT1 with values of 79 % and 66 %, respectively.

For the fuel blends, when compared to propane-only operation, there is an increase in efficiency for all operating points. In the data for the first blending ratio 75-25, when compared to pure propane, MGT1 displayed the largest increase of 22 % followed by MGT2 with 14 %, and 4 % for MGT3. In addition, there is a narrow combustion efficiency range for the fuel blends where MGT2 and MGT3 have a similar range of between 89 % and 93 %, and MGT1 is between 88 % and 90 %. For pure hydrogen operation, the combustion effiency of MGT1 and MGT2 is increased to 89 % from 66 % and 79 % and the efficiency of MGT3 is unchanged from propane-only operation at 89 %.

It is clear from the above observation that fuel blending with hydrogen improves combustion efficiency, especially at off-design operating points. This is attributed to the enhanced chemical kinetics of the H₂-O₂ reaction mechanisms. However, for the current combustor, the combustion effiency is not improved when operating close to the design point on pure-hydrogen.

Figure 11 is a compound figure. It presents mass average data (static temperature and species mass flow rate) calculated at the inlet and outlet planes of each combustion zone. The planes at the interfaces of zones are labelled with the abbreviations of both zones. All plots are concerned with running point MGT3 only.

Similar to the work of Tamang & Park (2023) , the CFD data shows that increasing the share of hydrogen reduces the mass fraction of CO₂ (Figure 11-c) and increases NO (Figure 11-b) at the outlet of the combustor. CO₂ decreases with increasing hydrogen content due to the reduced carbon mass in the fuel mixture, as hydrogen blending is optimized to maintain outlet gas temperatures equivalent to pure propane combustion.

There is no fuel bound nitrogen available to form NO, and the formation of Prompt-NO is suppressed due limited availability of C_xH_y radicals in the current combustor under the fuel-lean conditions. Therefore, the Thermal-Zeldovich route is the dominant NO formation pathway as tempertures in the PZ and SZ exceed 1500 °C [47]. This verifies the increased net formation of NO with increasing hydrogen as the effective flame temperature is simultaneously increased due to reduced flame radation [48]. In Figure 11-b for the mass flow rate of NO, the hydrogen only case is shown to produce 1.8 times more NO than the propane only case at the outlet of the combustor, similar margins to those reported in Benaissa, et al. [46]and Tamang & Park [5]. Hydrogen addition causes the mass of NO to peak in the SZ before sharply declining in the DZ which is a shift upstream relative to the results for pure-propane operation. It occurs because hydrogen's enhanced combustion kinetics, driven by increased OH, H, and O radical concentrations, accelerate NO formation through thermal, NNH, and N2O pathways compared to propane combustion [49].

For propane only operation, NO formation in the PZ and SZ occurs at a lower rate when compared to fuel blend and pure hydrogen operation (Figure 11-b) despite similar zone temperatures (Figure 11-a). A fuel rich PZ and SZ, for pure propane firing, increases the presence of C_xH_y radicals necessary for NO reburn, and a high temperature would promote Thermal NO formation [47]. These are competing mechanism and would explain the reduced NO formation rate. As C_xH_y radicals are consumed in the SZ and DZ of the combustor, the NO reburn rate will decrease increasing the dominance of the thermal pathway. However since dilution air entering the DZ, after plane SZ-DZ, lowers fluid temperatures below the threshold for Thermal NO formation, other NO formation pathways like the Prompt, N₂O and NNH will be prominent.

Less CO₂ is formed for blended fuels (Figure 11-c), however, an increased share of hydrogen in the fuel is found to reduce conversion rates of CO to CO₂. This is evident in Figure 11-d, where the outlet CO mass flow rate for mixture fuels is found to be 1.38 (25-75) to 1.7 (50-50) times more than the propane-only case. Despite being an indication of incomplete combustion, the large variation and concentration of fuel intermediates (not shown here) is found to be beneficial for CO conversion in propane-only combustion. At the low temperatures in DZ, these hydrocarbons increase CO conversion to CO₂ primarily through indirect pathways with O and H radicals [28]. With a much larger O an H radical pool in fuel mixture combustion, an increased rate of hydrogen abstraction subverts this effect by forming fewer and lighter hydrocarbons which limit the available CO to CO₂ reaction pathways.

Figure 12 shows the area-averaged wall temperatures for the simulated cases and blending ratios per case. Wall temperatures for propane operation show a monotonic upward trend, with PZ and SZ walls averaging 734-754 K and DZ walls ranging 1031-1038 K across all operating points. Both hydrogen mixtures and pure hydrogen reduce zonal wall temperatures compared to propane. The most significant reduction occurs in the PZ of the 50-50 mixture for MGT1 (normalized temperature of 0.71), while the smallest reduction is in the SZ of the 75-25 and 25-75 mixtures for MGT3 (0.884).

The DZ wall is highlighted by the rectangles in Figure 13. Static temperature contours reveal that propane-only operation exhibits active combustion along the DZ wall due to fuel accumulation in this region. Hydrogen introduction increases fuel velocity through the flame tube centre, preventing wall-bound fuel accumulation and combustion. This change, combined with reduced flame radiation from increased water vapor content in hydrogen combustion products, successfully lowers PZ and SZ wall temperatures under blended firing configurations [48, 50].

3.2. Thermofluid Network Model Validation Using CFD Results for Propane-Only Firing

In this section the network model is validated by benchmarking it against CFD results for propane-only operation using the running points found in Table 1. Validation is achieved through a comparison of the mass flowrate distribution, outlet total temperature, change in absolute total pressure and the mass fractions of major and minor species at the combustor outlet. Subsequently, the network model’s ability to predict combustor performance on mixed fuels is determined by comparison with the CFD data presented in the previous section.

The CFD results for outlet total temperature (OTT) and absolute pressure (OTP) is shown in the left-hand plot of Figure 14 for each running point. The error of the network model’s estimates for OTP and OTT, relative to the CFD, is shown in the right-hand plot. The network model predicts higher outlet temperatures than CFD but maintains consistent trends, showing temperature increases from MGT1 to MGT3 with less than 0.5% relative error. Outlet temperatures differ by an average of 5.69 K ± 0.29 K across all operating points. Inlet pressure predictions show excellent agreement with CFD, with relative errors ranging from 0.14% (MGT1) to 0.26% (MGT3).

Figure 15 presents the mass flow rates for each combustor flow path from both CFD and network model simulations, along with the corresponding total pressure drops per flow path. The total pressure drop across each component and the resulting mass flow rate is determined by the loss models presented in Section 2.2.1. The network model underestimates total combustor pressure drop by 3-4% compared to CFD. While cooling holes and pre-diffuser show the largest deviations (6% and 8% respectively), they contribute only 6.3% of total pressure drop. The major pressure drop components (snout, SZ jets, and PZ jets) show better agreement with maximum deviations of 3%, 3%, and 6% respectively. Despite the difference in the change in total pressure across each sub-component, the mass flow distribution predicted by the network model closely agrees with the CFD results.

Averaged inner wall temperatures (left) and zone temperatures (right) for each combustion zone, as estimated by the CFD, are shown in the top row of compound Figure 16. The bottom row shows the errors of the network model relative to the CFD.

The average combustion zone temperature estimated by the network model displays a similar trend to the CFD. Zone temperatures increase from the primary to the secondary zone and subsequently decreases in the dilution zone. This behaviour is expected as the primary zone jets supply air to complete combustion in the secondary zone increasing the average zone temperature, and the dilution zone jets (DZ) then quench the gas entering the dilution zone. Despite similar trends to the CFD data, the network model estimates lower average zone temperatures for all running points. The network model underpredicts the PZ temperature with a percentage error range of -9.6 % to -18.3 %, and a temperature difference of between -109 K to – 238 K. The SZ has an error range between -11.8 % (-197.5 K) and 1.2 % (19.2 K). The DZ is well predicted with an error of -3.4 % or -36.4 K

There is good agreement between the CFD and network model for the inner wall temperatures of the PZ and SZ for MGT1 and MGT2, with errors between -2.5 % and 7.3 %. For MGT3, the SZ wall temperature correlates well with an error 6.4 %, however the PZ wall temperature is under predicted by 15.5 %. Large errors are observed for the DZ wall temperature of all running points with an error range between –16.2 % and – 21.13 %.

The current network model is configured to estimate the convection heat transfer between a combustion zone, represented as a single flow element, and the solid node of the conduction heat transfer element. As such the difference between the average combustion zone temperature and the temperature of the solid node is used to estimate convection heat transfer. In contrast, considering a FVM CFD simulation, convective heat transfer is determined between the wall and the neighboring cell(s) only. Therefore, CFD simulations can account for localized changes in the flow field and combustion, and their effects on gas temperatures and heat transfer. This explains why, despite the low errors in the estimates for the DZ temperatures, there is a large error in the wall temperature estimates.

The network model has demonstrated satisfactory agreement with the CFD model in predicting the global operating conditions, internal aerodynamics, and heat transfer within a can combustor. The distinguishing feature of the network model in this study is the coupling of the aerodynamics with the chemical reactors. The following section of the validation is performed to show that by including chemical reactor models, the network model can predict changes in major and minor combustion species, and the pollutant NO across a system operating curve. Doing so positions the network model as a tool to improve combustion efficiency and pollutant emissions.

The top row of Figure 17 shows plots of CFD data only for the major species mass fraction (left), and the minor and pollutant species mass fraction (right) for each running point measured as a mass average across the combustor outlet plane. The bottom row are the corresponding plots for the error of the network model relative to the CFD (left) and the normalized mass fractions (right).

For the major species, CO₂ and H₂O, there is good agreement between the network model and the CFD. CO₂ is underpredicted within a range of -1.1 % to -2.4 %. The error for H₂O mass fraction is negligible and less than 0.03 %. A low degree of error for the major species indicates that similar degrees of combustion completeness have been predicted.

There is a considerable difference in the results for NO mass fraction seen in the right-hand plot of Figure 17. Despite having a similar trend, a mass fraction of 8.2 (MGT3) to 10.7 (MGT1) times higher than the CFD is estimated by the network model. This difference is attributed to primary and secondary combustion zones being approximated as well-stirred reactors which assume combustion proceeds uniformly. As a result, the elevated temperature throughout the reactor increases the thermal NO formation rate. In 3-D CFD, combustion models like the FGM used in this study account for the influence of local turbulence-chemistry interactions on combustion rates and may result in reduced NO formation rates when compared to a well-stirred reactor.

3.3. Comparison Between Thermofluid Network and CFD Model Results for Blended Firing

In the previous section, the network model was validated against CFD results for propane-only operation and demonstrated the network model’s capability to predict combustor performance using a fuel it was designed for. In this section, the predictive capability of the network model will be assessed. It will examine the model’s capacity to extrapolate combustor performance for fuel blends of propane (C₃H₈) and hydrogen (H₂), by validating once again to CFD results. The validation will be executed for each running point and fuel composition discussed in Section 2.3.1. As in Section 3.2, this validation will focus on the global performance of the combustor.

With data for all running points and fuel mixtures, Figure 18 is a compound figure displaying the following: labelled 1a is the total outlet temperature for the CFD simulation and 1b is the relative error of the network model. Similarly, 2a and 2b is a plot of the major species mass fractions for CO₂ and H₂O at the combustor outlet with the corresponding relative error for each plot shown in 3a and 3b, respectively.

For the outlet total temperature, a close agreement is observed between the network model and the results of the CFD simulation for all running points and fuel mixtures. The error relative to the CFD is within ± 0.55 %. This result is expected as Figure 9 showed that the Damköhler number is much less than 1 for each combustion zone and fuel mixture. This indicates that each combustion zone proceeds as a chemically limited system similar to the reactors employed to model the combustor.

For all fuel mixtures, as shown in Figure 18- 2a and 2b, the major species mass fractions (CO₂ and H₂O) at the outlet of the combustor as estimated by the network model corresponds well with the CFD results. For CO₂, the error is within the range of -2.4 % to 2.1 %. The H₂O mass fraction correlates well and it has an error range of -0.6 % to 1.3 %. This close correlation indicates similar degrees of combustion completeness and heat release, verifying the agreement in total outlet temperature.

In top row of Figure 19, the results for the minor species mass fraction CO (left) and NO (right) are presented as determined at the outlet of the combustor. In the bottom row are plots of the data normalized to the corresponding CFD results. The data labelled with the suffix F is from the network model, and the data without is from the CFD model.

It is observed that the CO mass fraction estimated by the network model is of a similar magnitude and follows the trend seen in the CFD data, decreasing with an increased share of H₂. The NO mass fraction is overestimated by the network model. Nonetheless, a trend of increasing NO mass fraction with an increase in H₂ is captured by the network model for all running points.

Notably, with an increase in H₂, the network model’s normalised mass fraction for CO is monotonically decreasing. It is also observed that for propane-only operation the model overestimates CO, and with the addition of H₂ CO is then underestimated. For NO, the error reduces initially and plateaus at the 50-50 mixture.

Under an assumption of an evenly distributed species composition in the chemical reactors, the fuel, O and H radical pool and intermediate hydrocarbons are mixed at a molecular level and reactions can proceed rapidly. Therefore, the reactions that are kinetically favoured at the reactor conditions will proceed. As such the intermediate hydrocarbons are rapidly oxidised and are found to have a concentration of near zero within all the network’s reactors. Therefore, in the DZ of the network model, there are insufficient intermediate hydrocarbons present to promote the low temperature conversion of CO to CO₂ observed in the CFD results for propane only operation. The result is a reduced conversion rate of CO to CO₂ for propane-only operation for the network model when compared to the CFD results. Similarly, unlike the observation in the CFD data, propane-only operation predicted by the network model is shown to produce more CO when compared to blended fuel operation.

With a low concentration of fuel intermediates, the rate of H₂-O₂ combustion is not limited by the kinetically favoured H abstraction reactions. Therefore, the radicals produced by the H₂-O₂ system are free to increase the CO conversion rate. This is evident in the normalised CO mass fraction for blended fuel firing that monotonically reduces from 4.2 (100-0) to a value of 0.25 for the 25-75 blending ratio.

Figure 20- 1a and 1b are plots for the CFD and network model wall temperatures respectively. In a similar order and presentation, the data for the averaged zone temperature is presented in Figure 20- 2a to 1b.

The stark increase in the PZ wall temperature in Figure 20-1b is reflected in the PZ temperature shown in Figure 20- 2b. Despite predicting similar SZ temperatures, the network model estimates for the SZ wall temperature are 1.3 to 1.4 times larger than the CFD for mixtures fuels. In contrast, the DZ zone and wall temperatures differ proportionately and increases to 1.07 and 1.08 times the CFD values for hydrogen only operation, respectively.

Similar to the results for propane only in the previous section, the considerable difference between the CFD and network model is due to the localised wall effects that a 3D FVM CFD simulation can resolve. This difference is intensified by the introduction of Hydrogen as combustion proceeds through the core of the combustor due to the fuel’s high inlet velocity. The flow along the walls is cooler as the intensity of the wall bound combustion is reduced and subsequently reduces the convective heat transfer rate. In addition, introducing hydrogen reduces the radiation from the flame which decreases the incident radiation [48,50].

The effect of the well-stirred assumption is evident in the network model’s PZ temperature estimate. From propane to hydrogen only operation, the average temperature of the PZ increases steadily from 1046 K to 1661 K. Subsequently this results in the elevated PZ wall temperature that is two times greater than the CFD value for hydrogen only operation.

3.4. Comparison Between CFD and Thermofluid Network model Computing Resources

The principal benefit of using a network model is the significant reduction in time. Table 6 presents the hardware used for each type of simulation and the average compute time for a single simulation to achieve convergence. These results highlight the benefits, specifically for system-wide simulations that the thermofluid network modelling approach offers.

4. Conclusions

This study details the development of a thermofluid network model employed to predict the performance of a MGT can combustor operating on pure and fuel blends of propane and hydrogen. It presents the treatment of sub-component pressure losses and heat transfer, and describes the first implementation of well-stirred and plug-flow reactors into Flownex SE to model combustion.

The CFD data showed that blended and pure hydrogen operation increases combustion efficiency and reduces the flame-tube wall temperatures. However, it increases pollutant emissions when compared to propane-only operation with 70 % and 80 % more CO and NO, respectively. The network model accurately predicts bulk flow distribution and pressure changes throughout the combustor. In addition, major species mass fraction and total temperature at the outlet of the combustor are predicted well for blended and pure fuels.

Nevertheless, the network model struggles to predict pollutant species (CO and NO) and flame-tube wall temperatures due to the simplifying assumption of uniform combustion zones. This highlights the need for increased fidelity in the reactor network to predict temperature distributions that better represent the 3D case. Despite this, the network model is shown to accurately predict the overall aerodynamic and thermodynamic changes in the combustor and does so 420x faster than the CFD model. This makes it highly suitable for rapid and accurate design space exploration using both pure and blended fuels.

Future work will involve implementing models to account for the effect of mixing rates on combustion to improve the prediction of pollutant emissions and, the influence of increasing the size of the reactor network on the accuracy of predicted flame-tube wall temperatures will be investigated. Moreover, the successful development of the combustor network model in Flownex SE renders it immediately suitable for integration into a full system-MGT model, enabling investigation into the effect of blended fuel firing on overall engine performance.