# Chemical timescale analysis of the Partially Stirred Reactor model for a hydrogen-fuelled scramjet

Alessandro Piscopo, Ward de Paepe, Alessandro Parente, Salvatore Iavarone

► To cite this version:

Alessandro Piscopo, Ward de Paepe, Alessandro Parente, Salvatore Iavarone. Chemical timescale analysis of the Partially Stirred Reactor model for a hydrogen-fuelled scramjet. Results in engineering, 2024, 23, pp.102834. 10.1016/j.rineng.2024.102834. hal-05224661

HAL Id: hal-05224661

https://hal.science/hal-05224661v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research paper

# Chemical timescale analysis of the Partially Stirred Reactor model for a hydrogen-fuelled scramjet

Alessandro Piscopo<sup>a,b,c</sup>, Ward De Paepe<sup>c</sup>, Alessandro Parente<sup>a,b,d,*</sup>, Salvatore Iavarone<sup>a,b,e</sup>

## ARTICLE INFO

### Keywords:

Supersonic combustion

Hydrogen

Scramjet

Reynolds-averaged Navier-Stokes

Partially Stirred Reactor

Chemical timescales

## ABSTRACT

Scramjet engines offer a promising option for powering hypersonic vehicles, although the unique characteristics of supersonic combustion present numerous challenges. Notably, turbulence-chemistry interactions are significant in supersonic flows, demanding specific combustion formulations in numerical simulations due to their dissimilarity from traditional combustion systems. This work evaluates the performance of the Partially Stirred Reactor (PaSR) closure of the mean chemical source term by performing Reynolds-Averaged Navier-Stokes (RANS) simulations of the German Aerospace Centre (DLR) scramjet. The PaSR approach requires defining a chemical timescale, and this study assesses several formulations. Among them, the Fastest Formation Rate (FFR) does not provide stable combustion, while the Slowest Formation Rate (SFR) achieve fairly good predictions. Moreover, the sensitivity to species selection in calculating the chemical timescale is investigated, as well as the impact of two PaSR residence time formulations. A newly proposed Batchelor-based mixing timescale performs well with both the SFR and FFR formulations and guarantees a coupling between chemical and mixing times. It also accounts for species diffusivity, which can have important effects on hydrogen-based fuel combustion.

## 1. Introduction

Supersonic flight remains an important research topic in the aero-space community, given its potential for several applications like space exploration, commercial air transport, and military strategy. The related industries have made sustained efforts to develop economical solutions featuring advanced vehicle designs and propulsion systems. Among the latter, air-breathing engines, and particularly scramjets, have gained interest and relevance mainly thanks to the fact that these engines use atmospheric oxygen for combustion, eliminating the need to carry on-board an oxidizer and thereby enhancing fuel efficiency and aircraft payload compared to conventional rocket engines. Moreover, a scramjet (a term that stands for supersonic combustion ramjet) generates thrust without any moving parts, resulting in a simplified layout and a lighter engine, and avoids the strong flow deceleration from supersonic to sub-sonic regime, which is a peculiar feature of a ramjet engine and a source

of performance decline at high Mach numbers ($M > 5$) [1,2]. Therefore, scramjets are ideal for supersonic/hypersonic flights due to their structural simplicity, high-speed impulse, and long-range efficiency. Despite the appealing features of scramjet engines, the implementation of supersonic combustion presents several challenges involving complex physical phenomena such as shock waves, compressibility effects, and ignition [1]. Moreover, efficient supersonic combustion is difficult to achieve due to the short residence times of the fuel in the high-speed flow and the consequent poor air-fuel mixing, for which suitable flame holders (such as cavities, struts, ramps, and back-facing steps) in scramjets have been investigated to provide stable combustion and boost combustion efficiency [3]. Research on scramjet engines is necessary to further improve their combustion efficiency, reduce their fuel consumption and operational costs, and ultimately unlock their potential to dramatically decrease travel times [4] and related carbon footprint of long-distance flights, and to power reusable hypersonic vehicles for

<footnote><sup>a</sup> Université Libre de Bruxelles, École Polytechnique de Bruxelles, Aero-Thermo-Mechanics Laboratory, Avenue Franklin Roosevelt 50, Brussels, 1050, Belgium</footnote>

<footnote><sup>b</sup> Université Libre de Bruxelles and Vrije Universiteit Brussel, Brussels Institute for Thermal-Fluid Systems and Clean Energy (BRITE), Avenue Franklin Roosevelt 50, Brussels, 1050, Belgium</footnote>

<footnote><sup>c</sup> University of Mons, Thermal Engineering & Combustion Unit, Rue de l'Epargne 53, Mons, 7000, Belgium</footnote>

<footnote><sup>d</sup> WEL Research Institute, Avenue Pasteur 6, Wavre, 1300, Belgium</footnote>

<footnote><sup>e</sup> Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, Gif-sur-Yvette, 91190, France</footnote>

<footnote>* Corresponding author at: Université Libre de Bruxelles, École Polytechnique de Bruxelles, Aero-Thermo-Mechanics Laboratory, Avenue Franklin Roosevelt 50, Brussels, 1050, Belgium.

E-mail address: Alessandro.Parente@ulb.be (A. Parente).</footnote>
sustainable space missions. To explore and optimise the performance of scramjet engines without the expense and challenges associated with extensive experimental testing, researchers resort to numerical simulations, particularly Computational Fluid Dynamics (CFD).

Previous LES studies demonstrated their ability to capture unique experimental trends, such as bi-modal velocity and temperature profiles, for different injection configurations and fuels [5–8]. Berglund et al. [5] emphasised the role of the recirculation regions and the limitations of applying flamelet models to supersonic conditions. While flamelet approaches with tabulated chemistry have shown good performances, they fail to reproduce the experimental trends at conditions outside the tabulated chemical manifold, which should then be extended [9]. This represents a critical modelling issue when replicating the scramjet performance over the whole flight history is necessary, along with predicting potential emissions of pollutants such as NOx in the scope of civil transportation, for which there is a scarcity of experimental data. Genin et al. [6] employed LES for a strut-based engine and highlighted the complexity of the mixing process and the need to model it accurately. Potturi et al. conducted different numerical studies on the DLR combustor [7,8], performing both RANS and hybrid LES/RANS simulations of the cold and reactive flows. The reactive flow was also analysed using different finite-rate combustion models, such as laminar chemistry, Partially Stirred Reactor (PaSR), and Scale Similarity (SS) [8]. The cold flow profiles obtained during the RANS simulations better matched with the experimental profiles, while the hybrid LES/RANS showed improved predictions for the reacting flow, as it includes the effects of local cooling and straining of the eddy structures. Still, the reaction rates were underestimated in the region downstream of the wedge base, providing lower temperature estimations. The RANS temperature profiles showed an important overprediction of the reaction region width. The later study [8] confirmed that the mixing process in the DLR combustor is dominated by large structures with important turbulence intensity.

Other LES investigations have stressed the importance of finite-rate chemistry effects and turbulence-chemistry interactions (TCI) in supersonic regimes [10,11]. Sun et al. [12] investigated the flame-holding characteristics of a cavity-based supersonic combustor by means of hybrid LES/RANS simulations with an improved blending function. They concluded that the shear layer was crucial for the complex flame-holding cycle involving unsteadiness and recirculation regions, emphasised by the intermittent combustion occurring in the cavity. The performance of finite-rate chemistry models was also studied by Xiang et al. [13], who used the PaSR and Unsteady PaSR (U-PaSR) approaches in RANS simulations of the German Aerospace Centre (DLR) and the NASA SCHOLAR combustors [14]. PaSR is a finite-rate, reactor-based combustion model that allows for detailed chemistry while maintaining a reasonable computational cost. It was shown that taking turbulence-chemistry interactions into account is a requisite for the accurate representation of the supersonic combustion physical processes. By considering the compressibility modification of Ingenito and Bruno [15], the computational cost was slightly increased and little improvements were observed, mainly in the velocity predictions. Fureby et al. [16] performed LES of the DLR scramjet [17] with the PaSR model, highlighting the role of the shock waves and their interactions with boundary layers and reactive zones. Indeed, the thickening of the boundary layer implies higher heat transfer rates and could lead to separation [18].

Although LES with finite-rate chemistry approaches is a powerful tool for supersonic combustion modelling, it requires significant computational resources. This limits the use of LES for tasks where several numerical simulations are necessary, such as parametric design exploration, optimisation, uncertainty quantification, and control. RANS simulations are suitable for addressing the challenges of supersonic combustion and designing scramjets in a more cost-effective way. While flamelet approaches may fail to accurately reproduce the experimental data in RANS context [19], the PaSR model has shown promising performances in different supersonic combustion applications [20,21]. Nevertheless, the number of RANS studies about supersonic combustion with

finite-rate chemistry remains limited, though it has shown adequate outcomes with simplified models [22,23], leaving margin for improvement. Evaluating and improving the PaSR model for RANS simulations could lead to more reliable and efficient designs of scramjets achieving stable and efficient combustion, high thrust, and potentially minimal pollutant emissions, contributing to advancements in aerospace propulsion. In a previous work [24], the authors assessed the performance of the PaSR combustion closure in RANS simulations of two selected test cases, namely, the Burrows and Kurkov's cavity-based combustor [25] and the DLR strut-based scramjet [17]. The impact of different turbulence models and formulations of the mixing timescale on the PaSR velocity, temperature, and species predictions was evaluated. Building upon this prior research, this work considers various sub-models for the chemical timescale, aiming to select the most appropriate formulation for representing the supersonic process. The DLR test case was retained as the strut-based fuel injection system has significant advantages over other injection systems since it creates a low-speed zone, increasing the residence time of the injected fuel, which helps flame stabilization. In addition to assessing the performance of different chemical timescale formulations, this paper also focuses on highlighting the sensitivity of the PaSR predictions to the species retained in the chemical timescale calculation and the reactor residence time estimation. Finally, an approach that guarantees a coupling between chemical and mixing timescale formulations is proposed and validated.

## 2. Partially Stirred Reactor (PaSR) model

In the realm of reactor-based models, the chemical source terms, which are returned to the CFD solver, can be expressed as the collective contribution from all potential composition spaces:

$$ \bar{\omega}_i = \int_{\Psi} P(\psi) \dot{\omega}_i d\psi. \qquad (1) $$

In Eq. (1), the i-th species chemical source term $\bar{\omega}_i$ is then formulated as an integral over the composition space, where $P(\psi)$ represents the joint scalar Probability Density Function (PDF) over the domain $\Psi$, and $\psi = [T, Y_i]^T$ denotes the composition space including temperature and species mass fractions [20]. Under certain conditions characterised by high Reynolds and Damköhler numbers, the PDF $P(\psi)$ takes on a bi-modal shape [26]. This PDF can then be decomposed into fine structures, indicated by (.), and their immediate surroundings, indicated by (.) (and assuming an infinitely thin flame front). The mean quantities for a given scalar $\phi$ of $\psi$ are computed as a weighted average of the corresponding values in fine structures and surroundings:

$$ \bar{\phi} = \gamma \phi^* + (1 - \gamma) \phi^0, \qquad (2) $$

where $\phi^*$ and $\phi^0$ are the fine structures and the surroundings scalar values, respectively. The term $\gamma$ is the cell reacting volume fraction and is determined by the interplay between a mixing timescale $\tau_m$ and a chemical timescale $\tau_c$ [27] as:

$$ \gamma = \frac{\tau_c}{\tau_c + \tau_m}. \qquad (3) $$

The Damköhler number, $Da = \tau_m/\tau_c$, being the ratio between mixing and the chemical timescales, indicates the importance of turbulence-chemistry interactions. It is assumed that most exothermic reactions occur within the fine structures [28,29], due to favourable mixing conditions, resulting in $\dot{\omega}_i(\psi^0) \sim 0$, so that $\bar{\omega}_i \sim \gamma \dot{\omega}_i(\psi^*)$ with $\psi^* = [T^*, Y_i^*]^T$ being the conditions within the fine structures. In the Partially Stirred Reactor (PaSR) model, the averaged chemical reaction rate for a given species $\bar{\omega}_i$ results from a mass exchange between the fine structures and its surroundings, and it is computed as:

$$ \bar{\omega}_i = \gamma \frac{\bar{\rho}(Y_i^* - Y_i^0)}{\tau^*}, \qquad (4) $$
where $\bar{\rho}$ is the Reynolds-averaged density, $Y_i^*$ and $Y_i^0$ are the i-th species mass fractions in the fine structures and in the surroundings, respectively, and $\tau^*$ represents the residence time in the fine structures. The fine structures are treated as ideal reactors, and two commonly used formulations, namely Perfectly Stirred Reactor (PSR) and Plug-Flow Reactor (PFR), have been employed in previous studies [30,31]. This work adopts a PFR approach for computational efficiency [32]. The fine structures are considered to evolve from an initial state $\psi^0 \sim \bar{\psi}$ for a duration equivalent to the residence time $\tau^*$. Integrating the PFR equations provides estimates for the sub-grid quantities $Y_i^*$. The residence time $\tau^*$ is calculated as the minimum of the mixing and chemical timescales, i.e., $\tau^* = \min(\tau_m, \tau_c)$ [33,34]. This formulation assumes that the residence time is constrained and governed by the fastest process between chemistry and turbulence.

Two different formulations for the chemical timescale $\tau_c$, already available in the literature [30,34–36], namely, the Slowest Formation Rate (SFR) and the Fastest Formation Rate (FFR), are assessed in this work, and one novel formulation, i.e., the Geometric mean of the Formation Rates (GFR), based on the previous two, is introduced. The SFR approach, which has already been validated on hydrogen-enriched non-premixed flames [35], considers that the slowest reacting species drive the reactions in the fine structures and computes the chemical timescale as:

$$ \tau_c^{SFR} = \max_{i=[1,N_S-N_D]} \left( \frac{Y_i}{|\dot{\omega}_i|} \right), \quad (5) $$

where $\dot{\omega}_i$ is the net formation/consumption rate of the i-th chemical species among $N_S$ species, $Y_i$ is the i-th species mass fraction, $N_D$ is the number of dormant species that are characterised by an absolute formation/consumption rate $|\dot{\omega}_i|$ lower than a given threshold $\epsilon = 10^{-8}$.

The second approach considers the Fastest Formation Rate (FFR) as the one driving the combustion process. In [34], it led to an accurate estimation of the chemical source terms for all the species in a turbulent, non-premixed syngas-air jet flame. The corresponding timescale is derived as follows:

$$ \tau_c^{FFR} = \min_{i=[1,N_S-N_D]} \left( \frac{Y_i}{|\dot{\omega}_i|} \right). \quad (6) $$

The latest formulation, named Geometric mean of the Formation Rates (GFR), considers the geometric mean of the first two approaches, providing timescale estimates between those of the SFR and FFR approaches, which correspond to the chemical time of the fastest species among the slowest ones. The GFR formulation addresses the limitations of SFR and FFR, which consider the slowest or the fastest species as the leading scale, respectively, and might provide unrepresentative values for most species and reactions involved. Thus, GFR alleviates this issue by estimating intermediate values possibly more representative of the chemical process, expressed as:

$$ \tau_c^{GFR} = \sqrt{\tau_c^{SFR} \cdot \tau_c^{FFR}}. \quad (7) $$

Several mixing time formulations can be used in the PaSR framework, and even though this work focuses mainly on the definition of chemical timescale, three different formulations are briefly discussed here. The first formulation, the integral-based formulation ($\tau_{m,I}$), estimates the mixing time as a fraction of the integral timescale ($\tau_I$) related to the large eddy break-up time:

$$ \tau_{m,I} = C_{mix} \tau_I = C_{mix} \frac{k}{\epsilon}, \quad (8) $$

where $k$ is the turbulent kinetic energy and $\epsilon$ its dissipation rate, and $C_{mix}$ is a constant parameter, which can be case-dependent and must be tuned accordingly [37]. The Kolmogorov timescale ($\tau_K$) is the smallest in the energy cascade model, and its associated mixing formulation ($\tau_{m,K}$) is defined as:

$$ \tau_{m,K} = \tau_K = \sqrt{\frac{\nu}{\epsilon}}, \quad (9) $$

<table><caption>**Table 1**
Inlet air and fuel conditions.</caption><thead><tr><th></th><th>T [K]</th><th>M [J]</th><th>Y<sub>H<sub>2</sub></sub></th><th>Y<sub>O<sub>2</sub></sub></th><th>Y<sub>H<sub>2</sub>O</sub></th><th>Y<sub>N<sub>2</sub></sub></th><th>p [Pa]</th></tr></thead><tbody><tr><td>Air inlet</td><td>340</td><td>2</td><td>0</td><td>0.232</td><td>0.032</td><td>0.736</td><td>101325</td></tr><tr><td>Fuel inlet</td><td>250</td><td>1</td><td>1</td><td>/</td><td>/</td><td>/</td><td>101325</td></tr></tbody></table>

where $\nu$ is the kinematic viscosity and $\epsilon$ is the turbulent kinetic energy dissipation rate. By taking the square root of the integral and the Kolmogorov scales, the geometric mean formulation can be defined:

$$ \tau_{m,G} = \sqrt{\tau_I \cdot \tau_K}. \quad (10) $$

An advantage of this formulation is avoiding the $C_{mix}$ coefficient tuning while providing similar estimates to the integral-based formulation [24]. Finally, a coupling between $\tau_c$ and $\tau_m$ is proposed in this paper based on the scale proposed by Batchelor [38], which depends on the molecular diffusivity $D$. The laminar diffusivity emphasises the effect of an individual species on the mixing timescale, whereas the turbulent diffusivity affects all species uniformly. The mixing timescale is then defined as proposed in [39], in analogy with the Kolmogorov approach:

$$ \tau_{m,B} = \sqrt{\frac{D_{i,m}}{\epsilon}}. \quad (11) $$

The approach thus relies on the same species for the determination of the relevant chemical and mixing timescales. The diffusion coefficient $D_{i,m}$, which is the diffusivity in the gas mixture of the species retained for the $\tau_c$ computation, indicated by the subscript $i$, is calculated as follows:

$$ D_{i,m} = \frac{1 - X_i}{\sum_{j,j \neq i} \left( \frac{x_j}{D_{ij}} \right)}. \quad (12) $$

In Eq. (12), $X$ is the species mole fraction, $D_{ij}$ is the laminar diffusivity of the retained species $i$ in the $j$-th species, and the summation is done over all the species in the mixture. Relying on a mixing scale as defined in Equation (11) can be important when molecular hydrogen or radical species are selected by the chemical timescale formulation. It must be noted that the proposed Batchelor-based formulation for the mixing timescale can only be coupled with either the SFR or the FFR formulation since the proposed GFR approach does not retain a chemical timescale linked to a specific species.

### 3. DLR scramjet

The DLR scramjet test engine experimental campaign was presented by Waidmann et al. [17]. The hydrogen fuel is injected parallel to the air stream through 15 1-mm diameter holes, located 2.4 mm apart at the base of a wedge with a half angle of 6°, as highlighted in Fig. 1a. The upper wall presents a 3° angle divergence to compensate for the boundary layer growth from $x = 58$ mm. The total length of the combustor is 340 mm, with a width of 40 mm and an inlet height of 50 mm. The strut injector helps in solving the well-known ignition problems in scramjet engines. The flame holder length is 32 mm, and its height is 6 mm [5]. A sketch of the DLR scramjet is shown in Fig. 1b with the measurement locations retaken. The air and hydrogen inlet conditions are summarised in Table 1. To reduce the computational cost, a smaller domain was kept, with only one injection hole and periodic boundary conditions applied in the spanwise $z$ direction. At the walls, no-slip boundary conditions were imposed. In addition, 12 boundary layers were added to the different walls with a first height of 0.02 mm, to keep the $y^+$ value in a reasonable range. The walls are assumed to be adiabatic due to the short time available for heat exchange. A flux-based pressure far-field boundary condition was chosen for the outlet. The available experimental data include the velocity profiles collected via Laser Doppler Velocimetry (LDV) and the static temperature profiles measured by Coherent Stokes
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) DLR geometry with the 15 fuel injectors and (b) sketch of the DLR scramjet combustion chamber adapted from [19].</figcaption></figure>

Raman Spectroscopy (CARS) at different axial locations. In addition, Particle Image Velocimetry (PIV) was used for flow visualization [17]. Cold flow pressure profiles are also available. These experimental data were used to validate the results obtained by simulations.

## 4. Numerical setup

This work presents RANS simulations of the DLR scramjet performed on the commercial software ANSYS Fluent R2022 R1. The DLR mesh was filled with the poly-hexcore method, and the result of a grid dependency study led to a final mesh composed of 950k cells. The chemical source reaction rate term was closed by implementing the PaSR model via a User Defined Function (UDF). The Laminar Finite Rate (LFR) model was used to compare the results with the different PaSR sub-models and assess the importance of turbulence-chemistry interactions. This model simplifies the treatment of turbulent fluctuations by straightforwardly approximating the source term. It assumes that the mean chemical source term depends only on mean temperature, density, and species mass fractions as follows:

$$ \bar{\omega}_k = \dot{\omega}_k(\bar{\rho}, \bar{Y}_1, \dots, \bar{Y}_{ns}, \bar{T}) = M_{w,k} \sum_{r=1}^{N_R} \hat{R}_{k,r}(\bar{T}, \bar{Y}_k), \quad (13) $$

where $M_{w,k}$ is the molecular weight of the k-th species, and $\hat{R}_{k,r}(\bar{T}, \bar{Y}_k)$ is the molar net production rate of species k from reaction r. This approach might be suitable for specific applications, such as supersonic combustion and atmospheric boundary layers. However, it is not advisable for most combustion scenarios, where turbulence fluctuations play a significant role. Nevertheless, in situations characterised by low Damköhler numbers, indicating minimal turbulence fluctuations, this model's performance could rival that of more sophisticated models, as demonstrated by prior research [40].

The Li kinetic mechanism, involving 9 species and 21 reactions [41], was used and a comparison was made with the Zetterval (Z25) scheme, including 11 species and 25 reactions [42]. The results of this comparison are shown in Appendix A. At first, four species were retained for the computation of $\tau_c$: H$_2$, O$_2$, OH, and H$_2$O. The sensitivity to the species selection was also assessed by considering all the species (except N$_2$) for the calculation of $\tau_c$. Following Piscopo et al. [24], the geometric mean formulation of the mixing timescale was selected. To further assess the relation between the chemical and mixing timescales and the reaction rate, the Kolmogorov formulation was employed with the FFR formulation. The Reynolds Stress Model (RSM) was selected as the turbulence model, with the Baseline (BSL) stress and Generalised $k-\omega$ (GEKO) options [43–45], also including compressibility effects. A density-based solver with an implicit formulation, characterised by a broader stability range with respect to the explicit formulation, was used. The flow field was initialised with the full multi-grid initialisation, and the second-order upwind scheme was employed for spatial discretization, as well as for the turbulent kinetic energy and the dissipation rate. Turbulent Prandtl and Schmidt numbers were fixed at 0.85 and 0.7, respectively.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Numerical density gradient contour (top) and experimental shadow image [19] (bottom) of the non-reactive case.</figcaption></figure>

## 5. Results and discussion

The comparison between the simulation results and the experimental data is first presented for the non-reactive case. Then, the performances of different assumptions for chemistry and its timescales in the PaSR model for RANS are assessed.

### 5.1. Non-reactive flow

Fig. 2 shows the simulated and experimental Schlieren images of the non-reactive flow in the DLR combustor. The tip of the wedge is the starting point of two shock waves, which are reflected at the lower and upper walls. With its 3° angle, the upper wall reflects the shock wave further downstream on the subsonic H$_2$ jet than the lower wall (1a in Fig. 2). In the experiment, the shock reflection occurs at approximately $x = 105$ mm, while in the CFD, it is estimated to be around $x = 102.5$ mm. The experimental value at which the intersection of the shock reflected on the lower wall and the jet (1b) is located also differs from the CFD one by 2.5 mm. The jet acts as a compression ramp, making the shock wave intersect with the recompression shock near $x = 130$ mm on the upper wall (2). The discrepancy between the experiment and the simulation is found in the intensity of the two shocks generated at the base of the wedge (near $x = 65$ mm) that are strongly underestimated in the simulation, resulting in the absence of the subsequent shock. The discrepancies likely arise from the limitations of the turbulence model in accurately capturing these complex flow interactions. The origin of these shocks could be from the non-uniformity of the wedge base, according to [19].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of cold pressure profiles with the experiments at the walls (a) and along the centreline (b) y = 25 mm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of cold velocity profiles with the experiments at four different axial locations: (a) 78 mm, (b) 125 mm, (c) 157 mm and (d) 233 mm.</figcaption></figure>

The cold flow pressure profiles are shown in Fig. 3, at the lower and upper wall (Fig. 3a), and along the centreline (Fig. 3b). At the lower wall, the CFD results are reasonably in agreement with the experiments, while the upper wall measurements are not available and only the numerical profile is shown for comparison with the shockwaves. At the latter location, it can be observed that the first pressure increase is slightly underestimated and the intensity of the second shockwave is not caught properly.

In Fig. 4, the cold flow velocity profiles at four different axial locations are shown. The profiles from the simulations are compared with the experimental data, and a reasonable agreement can be observed. At the first two locations, $x = 78$ mm (Fig. 4a) and $x = 125$ mm (Fig. 4b), the numerical profile captures the trend of the experimental points. At $x = 157$ mm (Fig. 4c) and $x = 233$ mm (Fig. 4d), a slight shift appears between the simulated and experimental velocity profiles, which is noticeable in the minimum speed location (2-3 mm downstream in the simulations with respect to the experiments), also appearing in the RANS of Oevermann et al. [19] and in the LES of Berglund et al. [5].

## 5.2. Reactive flow

Fig. 5 shows the experimental shadow image and the numerical density gradient contour in the reactive flow case. In both instances, the release of combustion heat expands the low-density area within the designated region. The comparison between the numerical and experimental cases demonstrates a robust agreement, with the distinct shocks portrayed in the figures displaying reasonable similarities. The distance between the reflection points of the CFD and the experiment is now minimal. At the wedge's tip, oblique shocks are formed and are subsequently reflected by both the upper and lower walls following the strut's interaction with the boundary layer. This interaction thickens the boundary layer, which is a phenomenon capable of inducing separation. Additionally, robust expansion fans, commonly referred to as Prandtl-Meyer expansion fans, accompanied by recompression shocks, redirect the flow towards the centreline. The thickening and separation of the boundary layer lead to a marked increase in static temperature and pressure, subsequently driving a rise in the heat transfer rate. Moreover, the density contour visually represents the significant density fluctuations inherent
<figure><img src="image_6.png" /><figcaption>Fig. 5. Numerical density gradient contour (top) and experimental shadow image [19] (bottom) of the reactive case.</figcaption></figure>

in the supersonic combustion process, directly influencing temperature and pressure dynamics. As the analysis progresses along the engine, the intensity of the shocks demonstrates a noticeable decline, ultimately dissipating before reaching the outlet.

The velocity profiles of the reactive case obtained by the PaSR model with different chemical timescale formulations and the LFR model are shown in Fig. 6. The simulations employed the Li mechanism [41]. The FFR formulation leads to higher velocities than all the other definitions, which estimate similar trends. While the LFR and SFR approaches result in almost matching velocity profiles, the GFR one leads to a higher minimum speed in Fig. 6b and a lower one in Fig. 6c.

The corresponding temperature profiles are presented in Fig. 7, where the bimodal trend observed at the initial axial location ($x = 78$ mm) remains elusive across all PaSR submodels. In the work of Li et al. [46], a bimodal trend appears, with modified version of the SST $k-\omega$ and $k-\epsilon$ models affecting the penetration depth of the hydrogen jet. Still, the values of the local minima and maxima are not predicted correctly. The SFR formulation outperforms both the GFR and LFR at $x = 125$ mm and $x = 233$ mm. At the last axial location, the different models get closer to each other. This could indicate the weakening of turbulence-chemistry interactions when getting closer to the outlet. Nevertheless, the importance of considering turbulence-chemistry interactions appears clearly in Fig. 7a and Fig. 7b, where the LFR predicts important temperature peaks with respect to the PaSR models. The GFR model predicts reduced temperature values at the three axial locations.

Flame extinction occurs in the FFR case, justifying the very different velocity profiles obtained in Fig. 6. To understand the cause of the flame extinction, one has to focus on the different variables of Equation (4), particularly the cell reactive volume fraction (or fine structures volume fraction) $\gamma$ and the residence time $\tau^*$ within the fine structures, which are both functions of chemical and mixing timescales.

Fig. 8 shows the values along the centreline of mixing (Fig. 8a) and chemical (Fig. 8b) timescales obtained by the various chemical timescale definitions and the geometric mean mixing time. The values shown for the FFR approach are taken before flame extinction occurs, as the simulation was initialised from the converged LFR simulation. It can be observed that from $x = 75$ mm to $x = 80$ mm the $\tau_c$ values are very low, almost of the order $\mathcal{O}(10^{-8})$. Thus, three orders of magnitude separate the mixing and the chemical time. The much lower $\tau_c$ impacted the fine structures residence time $\tau^*$, selected as the minimum between $\tau_c$ and $\tau_m$. A decrease in the residence time has two effects: on one hand, it increases the reaction rate (as it appears in the denominator of Eq. (4)) but, on the other hand, it affects the difference between the mass fractions in the fine structures and the ones in the surrounding (which are the initial ones in the ideal reactor) as a small residence time does not

allow the species to react. Concurrently, the parameter $\gamma$ was also affected by the lower $\tau_c$ value and kept decreasing, making the reaction rates smaller, until flame extinction, where $\tau_c$ was artificially clipped to 1 for numerical stability. Two different zones can be considered for the evolution of the chemical time for the other two approaches: from $x = 80$ mm to $x = 140$ mm, the $\tau_c$ value increases from its minimal value having an order of magnitude $\mathcal{O}(10^{-6})$ up to $\mathcal{O}(10^{-4})$ for the GFR, and from $\mathcal{O}(10^{-5})$ to $\mathcal{O}(10^{-2})$ for the SFR. This first section corresponds to the start and core of the flame. From $x = 140$ mm, the value stays approximately constant at $\mathcal{O}(10^{-4})$ for GFR and $\mathcal{O}(10^{-2})$ for SFR. The mixing time values are $\mathcal{O}(10^{-5})$, and the profiles are similar for all the formulations before flame extinction occurs for the FFR approach.

### 5.2.1. Sensitivity to species selection

As described in Section 2, calculating the chemical timescale with the different approaches depends on the species for which the formation/consumption rate is computed. In this section, the sensitivity to the species selection is assessed. The simulations employing the geometric mean mixing timescale and the three different chemical timescale definitions are considered. The chemical timescale is calculated from the formation/consumption rates of all the species (except $N_2$).

The simulation results providing a stable flame are shown in Fig. 9. While the results with SFR and FFR remain identical to the previous ones, the GFR approach provides flame extinction. This is explained by the fact that considering all the species, some of them fast reacting, brings down the FFR timescale and, consequently, the final $\tau_c$ estimate in GFR. Indeed, as per the definition of the geometric mean, the $\tau_c$ value of the GFR is closer to the FFR estimate than the SFR estimate. When this occurs for GFR, the values of $\gamma$ and $\tau^*$ reduce and bring down the reaction rates for the above-mentioned reasons.

### 5.2.2. Effect of the residence time definition

The previous modifications affected the estimates of the chemical timescale $\tau_c$ and impacted the cell reactive volume fraction $\gamma$ and the residence time $\tau^*$. According to Equation (4), the reaction rate depends on $\gamma$, the species mass fractions, and the residence time. Modifying the residence time also affects the mass fraction in the fine structures: the longer the residence time, the closer the fine structure mass fraction is to its equilibrium value. Fig. 10 shows the impact on the temperature profiles when switching from $\tau^* = \min(\tau_c, \tau_m)$ to $\tau^* = \tau_m$ (indicated as $\tau_m^*$). The latter is the standard formulation of the PaSR employed in the literature [47]. Only the SFR formulation provides a stable flame, with the peak temperature depending on the residence time formulation employed. A lower peak is predicted when using the mixing time as the residence time, which is more significant at $x = 78$ mm (Fig. 10a). These results can be explained by Fig. 10d, which illustrates the mass fraction variation $\Delta Y = Y_{H_2O}^* - Y_{H_2O}^0$ divided by the residence time calculated with the two different assumptions (with both the SFR and GFR formulations). The reaction rate is proportional to this ratio, and the latter is higher when using $\tau^* = \min(\tau_c, \tau_m)$, thus resulting in higher temperature peaks. From $x = 75$ mm to $x = 130$ mm, the difference in the ratio ranges between one to two orders of magnitude, as the baseline approach results in a maximum ratio of the order $\mathcal{O}(10^5)$ and the modified one is of the order $\mathcal{O}(10^3)$. The profiles obtained with the two residence times are similar except for the peak temperatures. Thus, it seems that the induction zone (i.e., $x = 78$ mm), is mostly sensitive to the residence time formulation. From the comparison of the GFR formulations, it can be seen that the base case provides a stable flame and exhibits a higher ratio of the order $\mathcal{O}(10^4)$ than the one obtained with the modified residence time $\tau_m^*$, which provides extinction.

### 5.2.3. Effect of the Damköhler number

The previous sections revealed that the FFR approach could never predict a stabilised flame using the geometric mean formulation for the mixing timescale. This mixing time formulation was chosen based on a previous work where the SFR approach was employed to calculate
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Experimental and simulated velocity profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, (c) 157 mm, and along the centreline (d) $y = 25$ mm, obtained by RANS simulation with the LFR and the various PaSR sub-models. The LES results from [5] are also reported.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained with the LFR model and the PaSR model with the SFR and GFR formulations. The LES results from [5] are also reported.</figcaption></figure>

the chemical timescale [24]. A hypothesis was made in Section 5.2 that $\tau_{m,G} \gg \tau_c$, minimising the value of the reactive volume fraction $\gamma$, and thus the reaction rate, leading to flame extinction. To verify this hypothesis, the FFR formulation was coupled with the Kolmogorov formulation

for the mixing timescale, which provides the lowest values for such a timescale [24]. Fig. 11 shows the temperature profiles obtained using the FFR approach and the Kolmogorov mixing timescale $\tau_{m,K}$. It appears that the FFR formulation can now predict a stable flame, also in the sim-
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the PaSR (a) mixing timescale and (b) chemical timescale obtained with $τ_{m,G}$ and the three different $τ_c$ formulations along the centreline.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained by the PaSR model with the SFR formulation considering 8 species (i.e., all the species of the kinetic mechanism except N₂) for the calculation of $τ_c$.</figcaption></figure>

ulations considering the other residence time definition ($τ^* = τ_m$) and those retaining four and eight species for the $τ_c$ calculation. Considering all the species (except N₂) for the $τ_c$ calculation determines reduced peak temperatures and reaction rates compared to the other cases.

Fig. 12 shows the mixing (Fig. 12a) and chemical (Fig. 12b) time estimates along the centreline y = 25 mm. It can be observed that the mixing times are of the order of magnitude $O(10^{-6})$ while the chemical times vary from $O(10^{-7})$ to almost $O(10^{-4})$. The most impacting parameter in the FFR formulation appears to be the species selection for the $τ_c$ calculation, which also impacts the residence time $τ^*$, in turn impacting the mass fraction exchange term $ΔY_i = Y_i^* - Y_i^0$ of Equation (4). On the other hand, changing the residence time definition only slightly changes the timescale estimations. While the effect of the species selection on $τ_m$ is noticeable but limited, the impact on $τ_c$ is significant as it results in a very different profile along the centreline with respect to the baseline case ($τ_{m,K}$) and the case with the modified residence time ($τ^* = τ_{m,K}$). While the SFR approach mainly depends on the main species, the FFR is more sensitive to radicals and minor species selection (as is the GFR

formulation, although in a limited way) when coupled with the Kolmogorov formulation. Species selection emerges as the most influential parameter for both FFR and GFR formulations. In the FFR approach, considering exhausted species (i.e., species resulting in low $τ_c$ values) may compromise numerical stability. This impact is mitigated in the GFR formulation by incorporating stable species. Fig. 13 depicts the Damköhler values in the domain for the case employing $τ_{m,G}$ and $τ_c^{SFR}$ and the one using $τ_{m,K}$ and $τ_c^{FFR}$. The contours reveal that the values for both cases fall within the range of Damköhler numbers approximately of the order of $O(0.1-10)$, which promotes a stable flame. However, distinct distributions appear. The geometric mean with the SFR formulation exhibits two narrow zones with $Da = 4.5$ within the flame. In contrast, the Kolmogorov with the FFR formulation displays a U-shaped region with $Da = 14.5$, encompassing the flame ignition region.

### 5.2.4. Diffusivity-based mixing timescale

The previous sections emphasised the impact of different mixing and the chemical timescale estimates on the flame stability. This sec-
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained with the assumption $\tau^* = \tau_m$ for the PaSR model with the SFR formulation. (d) Effect of the residence time definition on the ratio between the H₂O mass fraction exchange and the residence time for the SFR and GFR formulations.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained with the FFR chemical timescale $\tau_c^{FFR}$, and the Kolmogorov mixing timescale $\tau_{m,K}$. The profiles obtained by the simulations retaining all the species for the $\tau_c^{FFR}$ calculation and the two residence time definitions are also shown.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparison of (a) mixing timescales and (b) chemical timescales, obtained with the Kolmogorov and the FFR formulations, along the centreline. The impact of the residence time definition and the species selection is also shown.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Damköhler number contour obtained with (a) $τ_{m,G}$ and the SFR formulation and (b) with $τ_{m,K}$ and the FFR formulation considering 4 species for the calculation of $τ_c$.</figcaption></figure>

tion presents the results obtained by the PaSR model employing the Batchelor-based mixing time formulation, which introduces a coupling between the mixing and chemistry timescales. The temperature profiles are shown in Fig. 14. Both the SFR and FFR formulations predict a stable flame, also when considering all the species to calculate $τ_c$, although with different reaction rates, as it can be observed in Fig. 14a where the peak temperatures given by the FFR formulation are the lowest. As previously mentioned, the GFR formulation cannot be considered with this mixing timescale. The flame appears thinner with the FFR formulation when including all the species, which is confirmed by the profile at $x = 125$ mm (Fig. 14b) due to the inclusion of minor species. At $x = 233$ mm (Fig. 14c), all the profiles get similar toward the outlet.

Fig. 15 shows the mixing and chemical time estimates. By comparing this figure with Fig. 8, the following observations can be made: (1) the mixing time formulation has a significant impact on the chemical estimates, (2) although the range of chemical times covered by $τ_c^{SFR}$ is lower than before, SFR still predicts a stable flame, and (3) $τ_c^{FFR}$ trend is similar to the one obtained in Fig. 8a but the corresponding $τ_m$ is an order of magnitude smaller and this provides a diverse outcome, namely, a stable flame.

The competition between the diffusive and the turbulent mixing can be analysed by looking at the mixing time estimates, shown in Fig. 15a and the mass diffusion coefficient $D_m$, shown in Fig. 16. In Fig. 16, the species retained for the calculation of $D_m$ (and $τ_c$) among the preselected four ones, i.e., H₂, O₂, OH, and H₂O, is also indicated. When SFR is used, H₂ and H₂O are retained as slowest reacting species in the flame region (up to about $x = 125$ mm), with H₂ prevailing mostly near the injection as it is the species with the highest mass fraction, although already reacting. Then, H₂O produced from the combustion process is picked as species for the calculation of $τ_c$. In the post-flame region, H₂ is totally consumed, the reactions stop, and SFR selects the species with the highest concentration, i.e., O₂. On the other hand, when FFR is used, OH is mostly retained in the flame region, being either the species with

the lowest concentration or the highest formation rate in this region. A rapid consumption of O₂ in the region close to the injection makes this species retained at around $x = 80$ mm. In the post-flame region, H₂ is retained as FFR picks the species with the lowest concentration in that region. It can be observed from Fig. 16 that the diffusivity values are of the same order of magnitude ~ $O(10^{-4})$ with both chemical timescale formulations. In both cases, the highest value $O(10^{-3})$ is reached when the species retained is hydrogen. The discrepancy between the diffusivity estimates becomes larger beyond $x = 115$ mm. By comparing Fig. 16 with Fig. 15a, the importance of diffusive and turbulent mixing processes can be understood (as per Equation (11)). It appears that, until $x = 120$ mm, the $τ_m$ profiles obtained with SFR and FFR (4 species) follow the trends of the $D_m$ profiles. After that, the diffusivity decreases while the mixing time estimates increase, meaning that the turbulent kinetic energy dissipation strongly increases (by viscous heating), decreasing the turbulent kinetic energy. Once the chemical reactions are completed, it becomes challenging to maintain a high level of turbulence, explaining the decay in the turbulent kinetic energy.

Fig. 17 shows the Damköhler values obtained with the Batchelor-based mixing timescale approach and the two chemical time formulations. With SFR (Fig. 17a), the obtained contour shape is quite different from the one obtained when SFR was coupled with the geometrical mean mixing timescale $τ_{m,G}$. Moreover, the values are one order of magnitude smaller. With FFR, the maximum value is closer to that obtained when FFR was coupled with the Kolmogorov mixing timescale $τ_{m,K}$ (i.e., 12.3 against 14). However, the region where this maximum is attained is smaller than before, while the shape remains similar. Overall, the Batchelor-based mixing sub-model achieves robustness and ensures a more direct coupling between the characteristic timescales of the PaSR model.

### 5.2.5. *Quantitative comparison*

The performances of different assumptions for the characteristic timescales in the PaSR model for RANS are assessed in the reactive
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained with both SFR and FFR chemical timescales and the Batchelor-based mixing timescale $\tau_{m,B}$. The profiles obtained by the simulations retaining all the species for $\tau_c^{SFR}$ and $\tau_c^{FFR}$ calculation are also shown.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison of (a) mixing timescales and (b) chemical timescales, obtained with the Batchelor-based mixing time approach $\tau_{m,B}$ and both the SFR and FFR formulations, along the centreline. The impact of the species selection is also shown.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Mass diffusion coefficient $D_m$ retained by the SFR and FFR formulations, among 4 species, along the centreline.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Damköhler number Da contour obtained (a) with $τ_{m,B}$ and SFR and (b) with $τ_{m,B}$ and FFR, considering 4 species for the calculation of $τ_c$ and the Li mechanism.</figcaption></figure>

<caption>**Table 2**</caption>

Performance of the different combustion models and sub-models as provided by the averaged relative error [%] (Equation (14)) calculated for the temperatures measured on the centreline. Two axial locations are considered: $x = 125$ mm, and $x = 233$ mm.

<table><thead><tr><th>τ<sub>c</sub></th><th>Species</th><th>Residence time</th><th>τ<sub>m</sub></th><th>Flame</th><th>Error [%]</th></tr></thead><tbody><tr><td>LFR</td><td>/</td><td>/</td><td>/</td><td>Stable</td><td>9.77</td></tr><tr><td>SFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Stable</td><td>6.24</td></tr><tr><td>SFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Stable</td><td>6.56</td></tr><tr><td>SFR</td><td>4</td><td>τ<sub>m</sub><sup>*</sup></td><td>τ<sub>m,G</sub></td><td>Stable</td><td>7.03</td></tr><tr><td>SFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,B</sub></td><td>Stable</td><td>9.91</td></tr><tr><td>SFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,B</sub></td><td>Stable</td><td>6.56</td></tr><tr><td>GFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Stable</td><td>17.65</td></tr><tr><td>GFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Extinguished</td><td></td></tr><tr><td>GFR</td><td>4</td><td>τ<sub>m</sub><sup>*</sup></td><td>τ<sub>m,G</sub></td><td>Extinguished</td><td></td></tr><tr><td>FFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Extinguished</td><td></td></tr><tr><td>FFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,G</sub></td><td>Extinguished</td><td></td></tr><tr><td>FFR</td><td>4</td><td>τ<sub>m</sub><sup>*</sup></td><td>τ<sub>m,G</sub></td><td>Extinguished</td><td></td></tr><tr><td>FFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,K</sub></td><td>Stable</td><td>8.24</td></tr><tr><td>FFR</td><td>4</td><td>τ<sub>m</sub><sup>*</sup></td><td>τ<sub>m,K</sub></td><td>Stable</td><td>6.1</td></tr><tr><td>FFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,K</sub></td><td>Stable</td><td>10.62</td></tr><tr><td>FFR</td><td>4</td><td>τ<sub>min</sub></td><td>τ<sub>m,B</sub></td><td>Stable</td><td>4.64</td></tr><tr><td>FFR</td><td>8 (∀ sp.)</td><td>τ<sub>min</sub></td><td>τ<sub>m,B</sub></td><td>Stable</td><td>15.53</td></tr></tbody></table>

case by calculating the averaged relative error metric, as proposed in [48]:

$$ \left| \frac{E}{y_e} \right| = \frac{1}{N} \sum_{i=1}^{N} \left| \frac{y_m - y_e}{y_e} \right|, \qquad (14) $$

where E is the error, $y_m$ is the model prediction, $y_e$ is the experimental value, and N is the number of experimental points. Equation (14) was applied to the temperatures measured on the centreline at the two last axial locations, namely $x = 125$ mm and $x = 233$ mm. The coordinate $x = 78$ mm was not considered as the error calculation is strongly impacted by the fact that the bimodal shape observed in the experiments is not captured by the various models. Table 2 summarises the averaged relative errors obtained from the various combustion models and sub-models investigated in this work. Errors from cases with an extinguished flame are not reported. The SFR formulation provides similar error values ranging between 6 and 10%. The highest error with the SFR approach is obtained with the Batchelor-based mixing time considering 4 species. When coupling the Batchelor-based timescale with the FFR formulation, the lowest error is attained with a value of 4.64%. Conversely, when including all the species in the $τ_c$ calculation with the former mixing time, the error decreases to 6.56% for SFR and rises to 15.53% for FFR. The GFR formulation has an error value of 17% when it provides a stable flame. Finally, the coupling of the FFR approach with $τ_{m,K}$ provides an error estimated between 6.1% and 10.62%. Coupling the chemical and mixing timescales via the Batchelor-based approach ensures robustness, since a stable flame is obtained with all the chemical timescale formulations, without compromising the accuracy, or even improving it, as observed for the FFR sub-model.

## 6. Conclusions

In this work, a sensitivity analysis of the Partially Stirred Reactor (PaSR) combustion model was conducted by considering several sub-modelling options in RANS simulations of a supersonic combustion chamber burning hydrogen. Three different chemical timescale formulations, i.e., Slowest Formation Rate (SFR), Geometric mean of the Formation Rates (GFR), and Fastest Formation Rate (FFR), were assessed while maintaining the same sub-model for the mixing timescale. Moreover, the impact of the species selected for the chemical timescale calculation and the residence time formulation on the PaSR predictions was also assessed. While stable flames were always obtained when resorting to the SFR formulation, the GFR and FFR approaches provided extinguished flames when the relative difference between chemical and mixing timescales was significant, strongly reducing the reaction rates. SFR provided better results than the laminar finite-rate (LFR) chemistry model, emphasising the importance of turbulence-chemistry interactions, and showed its robustness in predicting the temperature profiles fairly well, independently of the species retained to calculate the chemical timescale and the residence time formulation considered. When predicting a stable flame, the GFR approach did not provide better results than SFR. Also, GFR was more sensitive to the selection of radical species and the residence time formulation. Finally, the FFR formulation predicted a stable flame only when coupled with the Kolmogorov mixing timescale, which ensured Damköhler numbers of the same orders of magnitude as those obtained by SFR with the initial mixing timescale sub-model. These results emphasised the need for a coherent coupling between chemical and mixing timescale formulations, which was then implemented via a Batchelor-based mixing timescale approach. This formulation ensured stable flames for all the chemical timescale sub-models considered. Moreover, it provided the best results when combined with the FFR sub-model. This study outlines two viable strategies for the supersonic combustion case under investigation: (1) coupling the geometric mean formulation of the mixing timescale with the SFR approach, and (2) pairing the Kolmogorov mixing timescale formulation with the FFR approach. A third strategy was proposed via an approach that binds together mixing and chemical timescales formulations more effectively: a Batchelor-based sub-model that calculates the mixing timescale from the diffusion coefficient of the species selected by the chemical timescale approach, ensuring a robust coupling with both FFR and SFR. Future works will focus on Uncertainty Quantification studies to improve the fidelity of the PaSR submodels and validate the framework for LES simulations.

### CRediT authorship contribution statement

**Alessandro Piscopo:** Writing – original draft, Software, Methodology, Investigation, Conceptualization. **Ward De Paepe:** Supervision, Funding acquisition. **Alessandro Parente:** Writing – review & editing, Supervision, Methodology, Conceptualization, Funding acquisition. **Sal-**
<figure><img src="image_19.png" /><figcaption>Fig. A.1. Experimental and simulated temperature profiles at three different axial locations, i.e., (a) 78 mm, (b) 125 mm, and (c) 233 mm, obtained with the LFR model and the PaSR model with the GFR and SFR formulations when using the Z25 mechanism [42].</figcaption></figure>

**vatore Iavarone:** Writing – review & editing, Software, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

A. Piscopo and S. Iavarone acknowledge the financial support of the Fonds National de la Recherche Scientifique (FRS-FNRS). The authors thank Bruce Crawford and Didier Bessette for providing valuable insights that significantly contributed to the implementation of the Partially Stirred Reactor model in ANSYS Fluent.

## Appendix A. Impact of the kinetic mechanism

The sensitivity of the simulation results to the choice of the kinetic mechanism is presented and discussed in this appendix. The results obtained by the simulations employing the Zetterval mechanism [42] are therefore presented. Fig. A.1 compares the outcomes of the simulations featuring the different PaSR submodels and the LFR model to assess the performances of this mechanism with respect to the Li mechanism. It can be observed that the LFR peak temperature is strongly decreased with respect to the Li mechanism at $x = 78$ mm. The SFR formulation estimates a slightly lower peak temperature at that location, while the GFR predicts a higher value than what was obtained with the Li mechanism in Fig. 7a. At $x = 125$ mm, a wider reaction region can be observed for all the models with an increased peak temperature. Again, at $x =$

233 mm, all the models get closer to each other, with the LFR being the closest to the experimental peak temperature, followed by SFR and GFR, respectively. Using a different mechanism does not prevent flame extinction when the FFR formulation is considered.

## References

[1] V.A. Sabelnikov, V.V. Vlasenko, Combustion in Supersonic Flows and Scramjet Combustion Simulation, Springer, Singapore, Singapore, 2018, pp. 585–660.

[2] J.S. MacKay, R.J. Weber, An Analysis of Ramjet Engines Using Supersonic Combustion, Sep. 1958.

[3] L. Sam, P. Idithsaj, P.P. Nair, A. Suryan, V. Narayanan, Prospects for scramjet engines in reusable launch applications: a review, Int. J. Hydrog. Energy 48 (92) (2023) 36094–36111, https://doi.org/10.1016/j.ijhydene.2023.05.341, https://www.sciencedirect.com/science/article/pii/S0360319923027593.

[4] M. Berglund, K. Fedina, C. Fureby, J. Tegnér, V. Sabelnikov, Finite rate chemistry large-eddy simulation of self-ignition in supersonic combustion ramjet, AIAA J. 48 (2010) 540–550, https://doi.org/10.2514/1.43746.

[5] M. Berglund, N. Wikström, C. Fureby, Numerical simulation of scramjet combustion, Tech. Rep., Swedish Defence Research Agency - FOI, 2005.

[6] F. Génin, S. Menon, Simulation of turbulent mixing behind a strut injector in supersonic flow, AIAA J. 48 (3) (2010) 526–539, https://doi.org/10.2514/1.43647.

[7] A. Potturi, J.R. Edwards, LES/RANS simulation of a supersonic combustion experiment, in: 50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2012.

[8] A. Potturi, J.R. Edwards, Investigation of subgrid closure models for finite-rate scramjet combustion, in: 43rd Fluid Dynamics Conference, 2013.

[9] A. Saghafian, V.E. Terrapon, H. Pitsch, An efficient flamelet-based combustion model for compressible flows, Combust. Flame 162 (3) (2015) 652–667, https://doi.org/10.1016/j.combustflame.2014.08.007.

[10] A. Mura, A. Techer, G. Lehnasch, Analysis of high-speed combustion regimes of hydrogen jet in supersonic vitiated airstream, Combust. Flame 239 (2022) 111552, https://doi.org/10.1016/j.combustflame.2021.111552.

[11] W.Y. Kun Wu, Peng Zhang, X. Fan, Numerical investigation on flame stabilization in dlr hydrogen supersonic combustor with strut injection, Combust. Sci. Technol. 189 (12) (2017) 2154–2179, https://doi.org/10.1080/00102202.2017.1365847.

[12] M.-B. Sun, Z.-G. Wang, J.-H. Liang, H. Geng, Flame characteristics in supersonic combustor with hydrogen injection upstream of cavity flameholder, J. Propuls. Power 24 (4) (2008) 688–696, https://doi.org/10.2514/1.34970.

[13] Z. Xiang, S. Yang, S. Xie, J. Li, H. Ren, Turbulence-chemistry interaction models with finite-rate chemistry and compressibility correction for simulation of supersonic
turbulent combustion, Eng. Appl. Comput. Fluid Mech. 14 (1) (2020) 1546-1561, https://doi.org/10.1080/19942060.2020.1842248.

[14] A.D. Cutler, P.M. Danehy, R.R. Springer, S. O'Byrne, D.P. Capriotti, R. DeLoach, Coherent anti-Stokes Raman spectroscopic thermometry in a supersonic combustor, AIAA J. 41 (12) (2003) 2451-2459, https://doi.org/10.2514/2.6844.

[15] A. Ingenito, C. Bruno, Physics and regimes of supersonic combustion, AIAA J. 48 (3) (2010) 515-525, https://doi.org/10.2514/1.43652.

[16] C. Fureby, E. Fedina, J. Tegnér, A computational study of supersonic combustion behind a wedge-shaped flameholder, Shock Waves 24 (2013) 41-50, https://doi.org/10.1007/s00193-013-0459-2.

[17] W. Waidmann, F. Alff, U. Brummund, M. Böhm, W. Clauss, M. Oschwald, Experimental investigation of the combustion process in a supersonic combustion ramjet, Tech. Rep., DLR Institute for Chemical Propulsion and Engineering, Hardthausen, Germany, 1994.

[18] J. Quinlan, J.C. McDaniel, T.G. Drozda, G. Lacaze, J.C. Oefelein, A priori analysis of flamelet-based modeling for a dual-mode scramjet combustor, https://arc.aiaa.org/doi/pdf/10.2514/6.2014-3743, https://doi.org/10.2514/6.2014-3743, https://arc.aiaa.org/doi/abs/10.2514/6.2014-3743.

[19] M. Oevermann, Numerical investigation of turbulent hydrogen combustion in a scramjet using flamelet modeling, Aerosp. Sci. Technol. 4 (2000) 463-480.

[20] Y. Moule, V. Sabel'nikov, A. Mura, Modelling of self-ignition processes in supersonic non premixed coflowing jets based on a PaSR approach, in: 17th AIAA International Space Planes and Hypersonic Systems and Technologies Conference, American Institute of Aeronautics and Astronautics, 2011.

[21] W. Zhang, X. Zheng, P. Yan, R. Xue, Numerical simulation of scramjet engine with compressibility corrections on k-ε model based on openfoam, Int. Commun. Heat Mass Transf. 158 (2024) 107893, https://doi.org/10.1016/j.icheatmasstransfer.2024.107893, https://www.sciencedirect.com/science/article/pii/S0735193324006559.

[22] K.A. Verma, K.M. Pandey, M. Ray, K.K. Sharma, The numerical investigation of combustion performance of scramjet combustor with variation in angle of attack, Results Eng. 15 (2022) 100507, https://doi.org/10.1016/j.rineng.2022.100507, https://www.sciencedirect.com/science/article/pii/S2590123022001773.

[23] N. Bordoloi, K.M. Pandey, M. Ray, K. Kumar Sharma, Computational investigation to interpret the reacting flow conditions for two different fuels in the scramjet combustor, Results Eng. 15 (2022) 100568, https://doi.org/10.1016/j.rineng.2022.100568, https://www.sciencedirect.com/science/article/pii/S2590123022002389.

[24] A. Piscopo, S. Iavarone, M. Savarese, M. Riis, B. Crawford, D. Bessette, S. Orsino, A. Coussement, W. De Paepe, A. Parente, Mixing time scale analysis of the partially stirred reactor model for high-speed turbulent combustion of hydrogen in vitiated air, Acta Astronaut. 218 (2024) 70-89, https://doi.org/10.1016/j.actaastro.2024.02.009.

[25] M.C. Burrows, A.P. Kurkov, Analytical and experimental study of supersonic combustion of hydrogen in a vitiated airstream, Tech. Rep. NASA-TM-X-2828, NASA, 1973.

[26] K.N.C. Bray, J.B. Moss, A unified statistical model of the premixed turbulent flame, Acta Astronaut. 4 (3) (1977) 291-319, https://doi.org/10.1016/0094-5765(77)90053-4.

[27] J. Chomiak, A. Karlsson, Flame liftoff in diesel sprays, Symp., Int., Combust. 26 (2) (1996) 2557-2564, https://doi.org/10.1016/S0082-0784(96)80088-9.

[28] J. Chomiak, A possible propagation mechanism of turbulent flames at high Reynolds numbers, Combust. Flame 15 (3) (1970) 319-321, https://doi.org/10.1016/0010-2180(70)90014-3.

[29] M. Tanahashi, M. Sato, M. Shimura, T. Miyauchi, DNS and combined laser diagnostics of turbulent combustion, J. Therm. Sci. Technol. 3 (3) (2008) 391-409, https://doi.org/10.1299/jtst.3.391.

[30] Z. Li, A. Cuoci, A. Sadiki, A. Parente, Comprehensive numerical study of the adelaide jet in hot-coflow burner by means of rans and detailed chemistry, Energy 139 (2017) 555-570, https://doi.org/10.1016/j.energy.2017.07.132.

[31] M.T. Lewandowski, I.S. Ertesvåg, Analysis of the eddy dissipation concept formulation for MILD combustion modelling, Fuel 224 (2018) 687-700, https://doi.org/10.1016/j.fuel.2018.03.110.

[32] A. Péquin, M.J. Evans, A. Chinnici, P.R. Medwell, A. Parente, The reactor-based perspective on finite-rate chemistry in turbulent reacting flows: a review from traditional to low-emission combustion, Appl. Energy Combust. Sci. 16 (2023) 100201, https://doi.org/10.1016/j.jaecs.2023.100201.

[33] S. Iavarone, A. Péquin, Z.X. Chen, N.A.K. Doan, N. Swaminathan, A. Parente, An a priori assessment of the partially stirred reactor (pasr) model for mild combustion, Proc. Combust. Inst. 38 (4) (2021) 5403-5414, https://doi.org/10.1016/j.proci.2020.06.234.

[34] A. Péquin, S. Iavarone, R. Malpica Galassi, A. Parente, The partially stirred reactor model for combustion closure in large eddy simulations: physical principles, submodels for the cell reacting fraction, and open challenges, Phys. Fluids 34 (5) (2022) 055122, https://doi.org/10.1063/5.0090970.

[35] S. Iavarone, M. Cafiero, M. Ferrarotti, F. Contino, A. Parente, A multiscale combustion model formulation for nox predictions in hydrogen-enriched jet flames, Int. J. Hydrog. Energy 44 (41) (2019) 23436-23457, https://doi.org/10.1016/j.ijhydene.2019.07.019.

[36] M. Evans, C. Petre, P. Medwell, A. Parente, Generalisation of the eddy-dissipation concept for jet flames with low turbulence and low Damköhler number, Proc. Combust. Inst. 37 (4) (2019) 4497-4505, https://doi.org/10.1016/j.proci.2018.06.017, https://www.sciencedirect.com/science/article/pii/S1540748918302001.

[37] M. Ferrarotti, Z. Li, A. Parente, On the role of mixing models in the simulation of mild combustion using finite-rate chemistry combustion models, Proc. Combust. Inst. 37 (2019) 4531-4538, https://doi.org/10.1016/j.proci.2018.07.043.

[38] G.K. Batchelor, Small-scale variation of convected quantities like temperature in turbulent fluid Part 1. General discussion and the case of small conductivity, J. Fluid Mech. 5 (1) (1959) 113-133, https://doi.org/10.1017/S002211205900009X.

[39] P. McMurtry, Turbulence - ME 7960, Course outline.

[40] T. Poinsot, D. Vynante, Theoretical and Numerical Combustion, RT Edwards, Inc., 2005.

[41] J. Li, Z. Zhao, A. Kazakov, F.L. Dryer, An updated comprehensive kinetic model of hydrogen combustion, Int. J. Chem. Kinet. 36 (2004) 566-575, https://doi.org/10.1002/kin.20026.

[42] A. Vincent-Randonnier, V. Sabelnikov, A. Ristori, N. Zettervall, C. Fureby, An experimental and computational study of hydrogen-air combustion in the lapcat ii supersonic combustor, Proc. Combust. Inst. 37 (3) (2019) 3703-3711, https://doi.org/10.1016/j.proci.2018.05.127.

[43] B.E. Launder, Second-moment closure: present... and future?, Int. J. Heat Fluid Flow 10 (4) (1989) 282-300, https://doi.org/10.1016/0142-727X(89)90017-9.

[44] M. Gibson, B. Launder, Ground effects on pressure fluctuations in the atmospheric boundary layer, J. Fluid Mech. 86 (1978) 491-511, https://doi.org/10.1017/S0022112078001251.

[45] B. Launder, G. Reece, W. Rodi, Progress in the development of a Reynolds stress turbulence closure, J. Fluid Mech. 68 (1975) 537-566, https://doi.org/10.1017/S0022112075000184.

[46] Z. Li, J. Wang, X. Li, Application of hydrogen mechanisms in combustion simulation of dlr scramjet combustor and their effect on combustion performance, Fuel 349 (2023) 128659, https://doi.org/10.1016/j.fuel.2023.128659, https://www.sciencedirect.com/science/article/pii/S0016236123012723.

[47] R. Amaduzzi, A. Bertolino, A. Özden, R.M. Galassi, A. Parente, Impact of scalar mixing uncertainty on the predictions of reactor-based closures: application to a lifted methane/air jet flame, Proc. Combust. Inst. (2022), https://doi.org/10.1016/j.proci.2022.06.028.

[48] W.L. Oberkampf, M.F. Barone, Measures of agreement between computation and experiment: validation metrics, J. Comput. Phys. 217 (1) (2006) 5-36, https://doi.org/10.1016/j.jcp.2006.03.037.