# Influences of steps in hybrid rocket engines:Simulation and validation on simplified geometries

Christopher Glaser, jouke hijlkema, Jean-Yves Lestrade, Jérôme Anthoine

► To cite this version:

Christopher Glaser, jouke hijlkema, Jean-Yves Lestrade, Jérôme Anthoine. Influences of steps in hybrid rocket engines:Simulation and validation on simplified geometries. Acta Astronautica, 2023, 208, 10.1016/j.actaastro.2023.03.037. hal-04062608

HAL Id: hal-04062608

https://hal.science/hal-04062608

Submitted on 7 Apr 2023

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
<figure><img src="image_2.png" /></figure>

Research paper

# Influences of steps in hybrid rocket engines:
Simulation and validation on simplified geometries

C. Glaser *, J. Hijlkema, J.-Y. Lestrade, J. Anthoine

ONERA/DMPE, Université de Toulouse, F-31410 Mauzac, France

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

## ARTICLE INFO

**Keywords:**

Regression rate

Backward facing step

Forward facing step

Simplified geometries

Arrhenius law

## ABSTRACT

In this study, we investigate the effect of single steps inside the fuel grains of typical hybrid rocket motors numerically. Steps can be a promising candidate to increase the regression rate of hybrid engines. We show that simulations on simplified geometries yield surprisingly close numerical regression rate profiles when compared to experimental data. Moreover, there is no necessity to calculate the progression of the fuel port at different time steps, as a single simulation on the average diameter can match the experiments sufficiently well. This observation allows for larger parametric studies on optimal positioning of the steps while limiting computational resources. The potential to increase the average regression rate through steps can be predicted numerically, and the simulations are following the same trend as the experiments. The augmented mixing and heat transfer inside the turbulent boundary layer can be predicted. It is concluded that backward facing steps increase the regression rate proportional to the step height up to 10 mm, and are preferably utilized upstream. For forward facing steps, the regression rate augmentation is saturated starting with 5 mm step height, and the position should be further downstream.

## 1. Introduction

The unique features of Hybrid Rocket Engines (HREs) such as safety and simplicity make them promising candidates for low cost space access. The downsides like low regression rates, slower transients, possibly high fuel slivers and low maturity, however, need to be addressed. The field of regression rate enhancement is broad and includes a variety of methods and concepts. In a previous article [1], we grouped them into (a) chemical adjustments to the fuel, (b) advanced injection methods and (c) novel combustion chamber designs. Steps can be counted to combustion chamber designs that increase mixing and turbulence.

In this article, the effect of steps inside HREs is investigated by a parametric study using the multi-physics solver CEDRE of ONERA by simplifying the HRE geometry with academical computational domains. In evaluating several step configurations such as Forward Facing Step (FFS) and Backward Facing Step (BFS) with different heights, an important height-to-length relationship for the recirculation zone can be derived. Although such values, even for reactive flow, can be found in literature for BFS and FFS, an extensive case study for an HRE with steps (and its unique interdependence of fuel injection from the fuel surface, mixing, heat transfer and regression rate) is lacking to the knowledge of the authors. However, these relationships are of

high importance when it comes to distributing multiple steps along a predefined profile or geometry as we proposed in Ref. [2].

First, we give a brief overview into turbulence enhancing devices for HREs in the literature. In a next step, we elaborate the numerical set-up and the boundary conditions used for the simulations. Once the foundation for the parametric study is laid, we validate the simplified geometries approach with experiments performed at ONERA. Finally, we discuss the results of the parametric study that follows. This article serves as direct continuation of the preliminary work presented in Ref. [3].

## 2. Diaphragms and steps in the literature

The concept of steps or diaphragms in the fuel port to promote turbulence and mixing is not new. ONERA employed a diaphragm in the LEX (Lithergol Experimental) rocket in the 1960s [4], making it the first diaphragm in a hybrid rocket motor with flight heritage. Gany and Timnat [5] observed average regression rate increases of 50% when employing a diaphragm in the middle of their polyester fuel grain already in 1972. Grosse [6] conducted studies on diaphragm shape and position in paraffin HREs (a work that has been continued at University of Padua [7]). The common results show regression

<footnote>* Corresponding author.
E-mail address: christopher.glaser@onera.fr (C. Glaser).</footnote>
## Nomenclature

## Abbreviations

<table><tr><td>O/F</td><td>Oxidizer-to-fuel ratio</td></tr><tr><td>BFS</td><td>Backward Facing Step</td></tr><tr><td>FFS</td><td>Forward Facing Step</td></tr><tr><td>GSI</td><td>Gas Surface Interaction</td></tr><tr><td>HDPE</td><td>High-Density Polyethylene</td></tr><tr><td>HRE</td><td>Hybrid Rocket Engine</td></tr><tr><td>HYCAT</td><td>Hybrid with Catalyzer</td></tr><tr><td>LEX</td><td>Lithergol Experimental</td></tr><tr><td>ONERA</td><td>The French Aerospace Lab</td></tr><tr><td>RZL</td><td>Recirculation Zone Length</td></tr><tr><td>SST</td><td>Shear Stress Transport</td></tr></table>

## Symbols

<table><tr><td>δt<sub>end</sub></td><td>Relative error end time</td></tr><tr><td>δt<sub>start</sub></td><td>Relative error start time</td></tr><tr><td>δt<sub>b</sub></td><td>Relative error burn time</td></tr><tr><td>δ<sub>cons</sub></td><td>Conservative relative error estimation</td></tr><tr><td>δ<sub>cut</sub></td><td>Relative error cutting</td></tr><tr><td>δ<sub>scan</sub></td><td>Relative error scanning</td></tr><tr><td>δ<sub>total</sub></td><td>Total relative error</td></tr><tr><td>ω̇</td><td>Production rate</td></tr><tr><td>η</td><td>Radial coordinate</td></tr><tr><td>λ<sub>g</sub></td><td>Thermal conductivity gas</td></tr><tr><td>λ<sub>s</sub></td><td>Thermal conductivity solid</td></tr><tr><td>μ</td><td>Dynamic viscosity</td></tr><tr><td>ρ</td><td>Gas density</td></tr><tr><td>ρ<sub>s</sub></td><td>Solid fuel density</td></tr><tr><td>A</td><td>Area</td></tr><tr><td>A<sub>pyr</sub></td><td>Arrhenius pre-exponent</td></tr><tr><td>a</td><td>Marxman empirical constant</td></tr><tr><td>C<sub>norm</sub></td><td>Normalization coefficient</td></tr><tr><td>c<sub>p,fuel</sub></td><td>Heat capacity fuel</td></tr><tr><td>D<sub>0</sub></td><td>Initial diameter</td></tr><tr><td>D<sub>2</sub></td><td>Diameter after the step</td></tr><tr><td>D<sub>avg</sub></td><td>Average diameter</td></tr><tr><td>D<sub>avg,flat</sub></td><td>Flat average diameter for mesh</td></tr><tr><td>D<sub>avg,var</sub></td><td>Varying average diameter for mesh</td></tr><tr><td>D<sub>f</sub></td><td>Final diameter</td></tr><tr><td>D<sub>f,cut</sub></td><td>False diameter due to cutting</td></tr><tr><td>D<sub>m</sub></td><td>Diffusion coefficient</td></tr><tr><td>E<sub>a</sub></td><td>Activation energy</td></tr><tr><td>G<sub>fuel</sub></td><td>Fuel mass flux</td></tr><tr><td>G<sub>ox</sub></td><td>Oxidizer mass flux</td></tr><tr><td>ΔH<sup>°</sup><sub>f,fuel</sub></td><td>Formation enthalpy</td></tr><tr><td>h</td><td>Step height</td></tr><tr><td>h<sub>i</sub></td><td>Enthalpy species i</td></tr><tr><td>h<sub>s</sub></td><td>Enthalpy solid</td></tr><tr><td>I<sub>turb</sub></td><td>Turbulent intensity</td></tr><tr><td>L<sub>g</sub></td><td>Fuel grain length</td></tr><tr><td>ṁ<sub>ox</sub></td><td>Oxidizer mass flow</td></tr></table>

rate increases downstream the diaphragms by up to 90% for lab-scale motors and combustion efficiency increases from 76% to 95%. Zhang et al. [8] pointed out that the regression rate after the diaphragm peaks consistently at the point of the re-attachment of the flow. Kumar and

<table><tr><td>Δm<sub>f</sub></td><td>Fuel mass loss</td></tr><tr><td>N</td><td>Total number of species</td></tr><tr><td>n</td><td>Marxman empirical constant</td></tr><tr><td>n<sub>pixel</sub></td><td>Number of miscalculated pixels</td></tr><tr><td>P<sub>c</sub></td><td>Chamber pressure</td></tr><tr><td>Q<sub>rad</sub></td><td>Radiation</td></tr><tr><td>R</td><td>Gas constant</td></tr><tr><td>Re<sub>h</sub></td><td>Step Reynolds number</td></tr><tr><td>Δr</td><td>Average burnt thickness</td></tr><tr><td>Δr<sub>cut</sub></td><td>Average burnt thickness obtained from cutting</td></tr><tr><td>res<sub>scan</sub></td><td>Scanner resolution</td></tr><tr><td>ṁ</td><td>Regression rate</td></tr><tr><td>ṁ<sub>norm</sub></td><td>Normalized regression rate</td></tr><tr><td>ṁ<sub>true</sub></td><td>True regression rate</td></tr><tr><td>T</td><td>Temperature</td></tr><tr><td>t<sub>b</sub></td><td>Burn time</td></tr><tr><td>U</td><td>Free stream velocity</td></tr><tr><td>x/h</td><td>Length-to-height ratio</td></tr><tr><td>x/l</td><td>Step position</td></tr><tr><td>y</td><td>Mass fraction</td></tr></table>

Kumar [9] (numerically) and Dinesh and Kumar [10] (experimentally) took advantage of this effect: both teams placed multiple diaphragms in the fuel grain to further extend the effect of diaphragms. They evaluated that the diaphragms are best to be placed after the recirculation zone of the preceding diaphragm. Quantitatively speaking, for Kumar and Kumar [9] the ideal spacing was 8–10 times the height of the diaphragm in numerical simulations, which corresponds to the observed recirculation zone length (RZL). Interestingly, Dinesh and Kumar observed the optimal spacing to be about 2.7 times of the RZL.

Concerning the use of steps in HREs, first records of a hybrid rocket motor with a backward facing step can be found in Korting et al. [11] in 1987. It was shown that a BFS at the port entrance doubled the local regression rate after the step. Noteworthy, the peak of the regression rate enhancement appeared at the point of the re-attachment of the flow after the step, as has been pointed out for diaphragms also. Lee et al. [12] proposed the use of a backward facing step instead of a diaphragm to avoid the pressure oscillations often induced by diaphragms. The stepped geometry was achieved by combining two fuel grains with different inner diameters. With a 7.5 mm step, the regression rate increased by about 50% downstream the backward facing step. The motor in the study of Kamps et al. [13] uses also a step in the middle of the fuel grain to promote mixing. Sakashi et al. [14] used multiple steps with alternating inner diameters. Their so-called concave–convex design effectively resembles a sequence of BFS and FFS with a maximal height difference of 9 mm. In doing so, the regression rate could be increased experimentally by up to 100%. Kumar and Joshi [15] used 5 cylindrical fuel grains with either 25 mm or 35 mm inner diameter. Combining them, a BFS/FFS fuel port could be created that increased regression rates between 30–50%. Contrary to the multi-diaphragm research, for steps, the optimal spacing due to the RZL has not been exploited.

Previous investigations of diaphragms and other turbulence enhancing devices in HREs all traced the regression rate enhancing effect back to the recirculation zone. Moreover, it was concluded throughout literature that in order to deploy multiple devices (steps, diaphragms) in a single motor, the ideal spacing should be done so that the new obstacle does not intervene with the previous recirculation zone. Hence, the height-to length ratio is of upmost importance for distributing multiple steps along an optimized profile (or any profile for that matter).
## 3. Numerical methods

The simulations are carried out with the multi-physics solver CEDRE of ONERA. It is used in research as well as industry for energetics and propulsion applications and validated on a large database [16]. CEDRE has been successfully employed for numerical simulations and validations on hybrid engines in the past [17,18]. This section provides a brief overview of the underlying numerical models and assumptions of the hybrid rocket simulations. As this study is focused on the application of the model rather than the numerical implementation, the concepts are briefly described. For a detailed and thorough understanding it is highly recommended to refer to the work of Durand et al. [17,19] who implemented the HRE regression rate model into CEDRE.

The base-line of the simulation is based on the HYCAT (Hybrid with CATalyzer) facilities of ONERA. The HYCAT family uses hydrogen peroxide (H₂O₂) and High-Density Polyethylene (HDPE) as propellants. The hydrogen peroxide is decomposed with a catalyst into hot water steam and gaseous oxygen. For this reason, only the gas phase has to be considered for the simulations, decreasing also the numerical complexity. As for the pyrolysis products of HDPE, only gaseous ethylene (C₂H₄) is assumed. Together with the decomposed hydrogen peroxide this leads to the following chemical species in the simulations: H₂O, CO, CO₂, O₂ and C₂H₄. The heat capacity of the species are calculated as an eighth-order temperature polynomial. The Sutherland law is used for the dynamic viscosity and the thermal conductivity is expressed by the Eucken approach. Using a constant Schmidt number equal to unity, the kinematic viscosity is obtained [17]. Given the significant difference of fuel regression rate (mm/s) and oxidizer flow velocity (m/s), the flow field is considered as quasi-steady.

The simulations are considered 2D-axisymmetric. It has been shown in three-dimensional numerical simulations on the transient behavior of HREs that the flow-field displays non-symmetrical behavior, especially in the start-up phase [20]. Nonetheless, the assumption of a 2D-axisymmetric computational domain is still leading to a reasonable accuracy when compared to the experimental regression rate (4.5% on average when using a moving mesh in transient simulations [21]). Even when the flow field is disturbed by a turbulence inducing device (diaphragm), the assumption of axisymmetry in numerical simulations match experimental regression rates well [8]. Gong et al. [22] observed little variation of the experimental regression rate profiles in circumferential direction, a similar trend to what we also obtained while investigating our post-firing grains. Moreover, in our study, the oxidizer is axially injected, without swirl or complex injection patterns. DiMartino et al. [23] pointed out that the use of more complex injection patterns (such as radial injectors) are preferably investigated using three-dimensional simulations, however, simple axial injection supports the assumption of 2D-axisymmetry, as has been also stated by Durand et al. [17]. Since the goal of this research is to enable larger parametric studies while limiting computational cost and using steady state simulations on an averaged fuel port diameter, the assumption of axisymmetry can be justified. Moreover, when comparing the trend of our axisymmetric regression rate profiles with three-dimensional simulations such as in Ref. [20] at elevated burn-times, the trend of the regression rate profiles of both axisymmetric and three-dimensional simulations are resembling progressively with increasing burn time. We will see in Section 4 that our simulations are able to match the our experimental regression rates reasonably well.

The flow is described by the compressible Reynolds-Averaged Navier-Stokes equations given that the Reynolds number in typical HYCAT tests is in the order of 10⁵. The Reynolds stress tensor is modeled by the Boussinesq assumption. The Reynolds analogy and the postulation of a turbulent Prandtl number of 0.9 is used to model the turbulent heat flux. As for the turbulent mass flux, a constant turbulent Schmidt number of 0.9 and the Reynolds analogy allow to calculate the turbulent mass flow. With the generalized gradient diffusion hypothesis, scalar diffusion is established [17,19]. The set

of equations is closed with the Shear Stress Transport (SST) model of Menter [24], combining the benefits of the k-ε model in the core flow and the k-ω model near the walls. Due to the fuel injection at the fuel surface, the correction of Wilcox [25] for the specific dissipation at the fuel wall is considered. The fuel inlet is modeled with a Gas Surface Interaction (GSI) model using an Arrhenius-type description to calculate the regression rate based on the fuel surface temperature as is described in Section 3.2.

### 3.1. Chemical reaction

The combustion inside an HRE is described by a turbulent diffusion flame. The two-step kinetic chemical reaction denotes according to Westbrook and Dryer to [26]:

$$C_2H_4 + 2O_2 \rightarrow 2CO + 2H_2O \quad (1)$$

$$CO + 0.5O_2 \leftrightarrow CO_2 \quad (2)$$

The detailed rates of production and destruction can be found in Ref. [17].

### 3.2. Gas-surface interaction model

The GSI model at the fuel surface used in CEDRE by Durand et al. [17,19] is based on an energy and mass balance at the fuel surface. It couples the fuel surface temperature to the regression rate, allowing to calculate the regression rate rather than fixing it at the fuel surface. According to mass conservation at the fuel surface, the mass flux of the fuel ($G_{fuel}$) can be written as:

$$G_{\text{fuel}} = \rho_s \cdot \dot{r}, \quad (3)$$

with $ρ_s$ being the solid fuel density and $\dot{r}$ the regression rate. The production rates $\dot{\omega}_i$ of each species i are considered source terms and yield the mass balance:

$$y_i \cdot G_{\text{fuel}} - \rho D_{m,i} \frac{\partial y_i}{\partial \eta} = \dot{\omega}_i, \quad (4)$$

with $D_m$ being the diffusion coefficient, y the mass fraction and η the radial coordinate. We consider the sum of the mass source terms $\dot{\omega}_i$ to be equal to the total fuel mass flux and ethylene as the sole pyrolysis product:

$$\sum_{i=1}^{N} \dot{\omega}_i = G_{\text{fuel}}. \quad (5)$$

Moreover, the energy balance (gas phase left, solid phase right) denotes to:

$$\lambda_g \frac{\partial T}{\partial \eta} + Q_{\text{rad}} - \sum_{i=1}^{N} \dot{\omega}_i \cdot h_i = \lambda_s \left( \frac{\partial T}{\partial \eta} \right)_s - G_{\text{fuel}} \cdot h_s. \quad (6)$$

$\lambda_{g/s}$: Thermal conductivity gas/solid

$T$: Temperature

$Q_{\text{rad}}$: Net radiation

$h_{i/s}$: Enthalpy species i/solid

$N$: Total number of species

Assuming a quasi-steady heat flux through the solid fuel:

$$\lambda_s \left( \frac{\partial T}{\partial \eta} \right)_s = G_{\text{fuel}} \cdot c_{p,\text{fuel}} \cdot \Delta T, \quad (7)$$

and defining the solid fuel enthalpy ($h_s$) as:

$$h_s = c_{p,\text{fuel}} \cdot \Delta T + \Delta H_{f,\text{fuel}}^{\circ}, \quad (8)$$
<table><caption>**Table 1**
Properties of HDPE [17].</caption><thead><tr><th>A<sub>pyr</sub> [kg/m<sup>2</sup>]</th><th>E<sub>a</sub> [mol]</th><th>ΔH<sup>°</sup><sub>f,fuel</sub> [J/kg]</th></tr></thead><tbody><tr><td>3.5*10<sup>6</sup></td><td>125.6*10<sup>3</sup></td><td>-6.365*10<sup>5</sup> ± 1.14*10<sup>4</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Computational domain for backward facing step case.</figcaption></figure>

Eq. (6) simplifies (by neglecting radiative heat transfer) to [17]:

$$ \lambda_g \frac{\partial T}{\partial \eta} = \sum_{i=1}^{N} \dot{\omega}_i \cdot h_i - G_{\text{fuel}} \cdot \Delta H_{f, \text{fuel}}^{\circ}, \quad (9) $$

with $c_{p,fuel}$ and $\Delta H_{f,fuel}^{\circ}$ being the heat capacity and formation enthalpy respectively. Finally, the production rate of ethylene ($C_2H_4$) – hence the fuel mass flux $G_{fuel}$ – is given by the Arrhenius law:

$$ G_{\text{fuel}} = A_{\text{pyr}} \cdot \exp\left(-\frac{E_a}{R \cdot T}\right). \quad (10) $$

The values for the pre-exponent $A_{pyr}$, the activation energy $E_a$ and the formation enthalpy $\Delta H_{f,fuel}^{\circ}$ for HDPE are listed in Table 1.

### 3.3. Computational domain

The computational domain in this article is based on the HYCAT facilities at ONERA as described in Section 3. In order to decrease the number of influences on the results of the fuel surface and to decrease computational costs, the HYCAT family is approximated by a 2D-axisymmetrical simplified domain, as depicted in Fig. 1 for the BFS case. In Section 4, we will show that approximating the full motor geometry with an simplified approach greatly decreases the computational effort while yielding sufficient details for a parametric study. The open source software GMSH [27] is used to create the meshes. The total fuel grain length can vary between 240–500 mm. The inlet fixes the oxidizer mass flux ($G_{ox}$), the inlet temperature and the turbulent intensity ($I_{turb}$). The fuel surface is modeled according to the GSI presented in Section 3.2 and the outlet is a pressure outlet which fixes the pressure to a representative chamber pressure. The mesh is chosen to be Cartesian because it lessens the number of cells (as compared to triangles) and allows to refine the mesh size close to the fuel surface to approximate the boundary layer region. The minimal cell size at the fuel surface is calculated to ensure an $y+ < 1$ (around 1 µm) for all possible total cell counts. Euler’s implicit method with a fixed time step of $10^{-6}$ is used and the spatial discretization is second order accurate [17].

The inner diameter before the step is fixed to 25 mm as this represents the minimum diameter of fuel grains used in HYCAT. The diameter after the step is dependent on the step height ($h$). For a

maximum step size of 10 mm this leads to a diameter after the step of 45 mm. Changing from BFS to FFS domain is conveniently achieved by switching the inlet and outlet conditions. This signifies that for the FFS cases, the diameter after the step is always fixed to 25 mm and the diameter before the step changes according to the step height. Translated to a real world application this means that a cylindrical fuel grain with two different inner diameters is simply turned around.

### 3.4. Mesh convergence

In order to prove mesh convergence, six different meshes with increasing cell count (from roughly 15 000 to over 400 000) are created for a representative BFS domain with a step height of 2 mm. All meshes ensure a $y+ < 1$ at the fuel surface. The results are presented in Table 2. Mesh 1 (14 902 cells) serves as starting point and the cell count is approximately doubled with each mesh. After each individual simulation is converged at a fixed time step of $10^{-6}$ after around 120 000 iterations, the average values of surface pressure, surface temperature, total fuel mass flow and RZL are gathered for 10 000 iterations with 100 data points at the same physical time. Mesh 1 presents the base case, the values for meshes 2–6 give the relative change compared to the preceding simulation. For example, for mesh 4, the total fuel mass flow increased by 0.16% as compared to the fuel mass flow of mesh 3. The boundary conditions of the mesh convergence study are loosely based on the conditions inside a typical HYCAT firing and can be found in Table 3.

Looking at Table 2, it becomes obvious that the average values of pressure and temperature converge fast and could be considered mesh independent starting from mesh 1. For the fuel mass flow, however, the values converge between mesh 2 and 3. The largest changes in values can be seen for the RZL. After the first refinement to mesh 2, the length of the recirculation zone increased by almost 20%. This is due to the fact that the coarse mesh 1 is not able to properly capture the eddies that are separating after the step. After mesh 3, the RZL finds a more stable value. However, in order to correctly judge the mesh independence for the RZL, it is necessary to also consider the shape of the recirculation zone. Fig. 2 depicts the appearance of the recirculation zone for Meshes 1, 4 and 5. Here, the RZL for case 1 is not at all developed and too short. For mesh 4, however, the length and the height of the RZL have a more consistent appearance. Nonetheless, a smaller second recirculation zone inside the bigger zone is observable, which could have a non-negligible effect on the simulations. Moreover, the contour of the flame zone is less consistent as compared to mesh 5. Mesh 6 only changes the RZL by under 2%, but uses double the cells. For this reason, mesh 5 with 208 307 cells is considered as the most suitable mesh for the numerical test campaign in terms of mesh independence and computational cost.

Having set the level of mesh refinement, it is possible to perform an analysis of the accumulation of numerical errors as proposed for reactive gas dynamics in Ref. [28]. In our simulations, the ratio of actual timesteps to the number of theoretically allowable timesteps ensuring a maximum error of 5% is $180*10^3$ in the worst case geometry (smaller cell count or big mean cell ratio) and $763*10^3$ for the best case scenario (high cell count or small mean cell ratio). Thus, the reliability of the simulations in terms of accumulation of numerical errors [28] is given in all cases.

## 4. Validation of simplified geometries with experiments

Here, we elaborate on the assumptions made when using simplified geometries to represent the flow-field inside the hybrid engine. Therefore, the two types of geometries investigated are presented and the results obtained with these meshes compared to experiments to prove their validity.
<table><caption>**Table 2**
Mesh convergence study.</caption><thead><tr><th>Mesh</th><th>Cells</th><th>Pressure [Pa]</th><th>Temperature [K]</th><th>Fuel mass flow [g/s]</th><th>RZL [mm]</th></tr></thead><tbody><tr><td>1</td><td>14902</td><td>3602121</td><td>924.74</td><td>5.978</td><td>8.56</td></tr><tr><td>2</td><td>26199</td><td>+0.03%</td><td>+0.14%</td><td>+2.11%</td><td>+19.77%</td></tr><tr><td>3</td><td>56263</td><td>+0.03%</td><td>+0.03%</td><td>+0.56%</td><td>+2.91%</td></tr><tr><td>4</td><td>96683</td><td>± 0%</td><td>+0.01%</td><td>+0.16%</td><td>+3.77%</td></tr><tr><td>5</td><td>208307</td><td>+0.01%</td><td>+0.05%</td><td>+0.90%</td><td>-2.73%</td></tr><tr><td>6</td><td>412097</td><td>± 0%</td><td>+0.02%</td><td>+0.60%</td><td>+1.87%</td></tr></tbody></table>

<table><caption>**Table 3**
Boundary conditions of mesh convergence.</caption><thead><tr><th>T<sub>inlet</sub></th><th>P<sub>outlet</sub></th><th>G<sub>ox</sub></th><th>y<sub>H2O</sub></th><th>y<sub>O2</sub></th><th>I<sub>turb</sub></th></tr></thead><tbody><tr><td>1000 K</td><td>36 bar</td><td>200 kg/m<sup>3</sup></td><td>0.59</td><td>0.41</td><td>1%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Recirculation zone (numerical) after the step for different meshes.</figcaption></figure>

## 4.1. Definition of the simplified geometries

Instead of simulating the motor at different instances of time throughout the burn, the arithmetic mean between the fuel port diameter at the beginning ($D_0$) and the end of the burn ($D_f$) is chosen to limit the number of simulation per test case to one:

$$D_{\text{avg}}(x) = \frac{D_0(x) + D_f(x)}{2} \qquad (11)$$

This choice is justified, as it has been shown analytically that the average mass flux inside an HRE is best to be approximated by the average diameter throughout the burn [29]. Moreover, numerical simulations by other authors show satisfactory results of the average diameter approach even when compared to results obtained with considerably more time consuming and complex simulations with fuel surface shape changes across different time steps [30].

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The two types of meshes.</figcaption></figure>

The process to obtain the average diameter for the simplified geometries to be compared to experiments is as follows: after the burn, the fuel grain is cut along its longitudinal axis and the final fuel port diameter digitized using an office scanner and image processing. For details on this methods the reader is referred to Ref. [2]. With the help of the final diameter $D_f(x)$, the known initial diameter $D_0(x)$ and the burn time $t_b$ the local regression rate $\dot{r}(x)$ can be calculated as:

$$\dot{r}(x) = \frac{D_f(x) - D_0(x)}{2t_b} \qquad (12)$$

Moreover, the average diameter $D_{\text{avg}}$ for the preparation of the meshes is obtained in two different ways in order to evaluate the better fit with the experiments:

1. For the calculation of $D_{\text{avg, var}}(x)$, the true final diameter obtained from the cut fuel grains is used, thus forming a varying profile along the fuel port axis $x$.

2. The final fuel port profile is averaged along the profile leading to a flat average diameter $\bar{D}_{\text{avg, flat}}$.

In Fig. 3, these two ways to estimate the average diameter are depicted and the corresponding mesh fuel surface shown for the final diameter after the test of HYCAT 46.

The most prominent difference between both mesh designs is the missing injector effect between 0 m and roughly 0.05 m. In the experiments, this area regresses faster as the remaining fuel because of the formation of a recirculation zone and increased mixing [31] – the very effect we want to exploit by introducing BFS and FFS inside the fuel grain. In the varying mesh design, we model this cavity, whereas for the flat design this area vanishes due to the averaging process. Apart from the missing zone of injector effects, the two meshes are relatively
<table><caption>**Table 4**
Summary of experimental data.</caption><thead><tr><th>Test</th><th><i>D</i><sub>0</sub> [mm]</th><th><i>t</i><sub>b</sub> [s]</th><th><i>m</i>'<sub>ox</sub> [kg/s]</th><th>O/F [-]</th><th><i>G</i>'<sub>ox,avg</sub> [kg/m<sup>2</sup>s]</th><th><i>r</i>' [mm/s]</th><th>Δ<i>m</i>'<sub>f</sub> [g]</th><th><i>P</i>'<sub>c</sub> [bar]</th></tr></thead><tbody><tr><td>H46 (Ref.)</td><td>25</td><td>9.41</td><td>0.303</td><td>16.49</td><td>371.23</td><td>0.77</td><td>172.7</td><td>20.4</td></tr><tr><td>H47 (Ref.)</td><td>40</td><td>9.28</td><td>0.329</td><td>26.10</td><td>220.11</td><td>0.39</td><td>116.8</td><td>20.0</td></tr><tr><td>H48 (Ref.)</td><td>25</td><td>9.65</td><td>0.337</td><td>7.96</td><td>388.02</td><td>0.86</td><td>408.9</td><td>26.1</td></tr><tr><td>H49 (Ref.)</td><td>40</td><td>9.65</td><td>0.348</td><td>12.58</td><td>228.02</td><td>0.42</td><td>266.7</td><td>24.1</td></tr><tr><td>H50 (FFS)</td><td>40/25</td><td>9.60</td><td>0.343</td><td>10.19</td><td>299.37</td><td>0.59</td><td>322.8</td><td>25.4</td></tr><tr><td>H52 (BFS)</td><td>25/40</td><td>9.66</td><td>0.346</td><td>8.84</td><td>289.59</td><td>0.67</td><td>378.2</td><td>28.3</td></tr></tbody></table>

similar. Nonetheless, the flat surface mesh could theoretically also be obtained with the mass loss method using the fuel mass loss ($\Delta m_f$), fuel density ($\rho_f$) and the total grain length ($L_g$) [32]:

$$D_f = \sqrt{D_0^2 + \frac{4\Delta m_f}{\pi \rho_f L_g}}, \quad (13)$$

which would eradicate the necessity to cut the fuel grains. For this reason it is of interest to compare the results of both meshing approaches with the experiments.

## 4.2. Experimental data and error estimation

A total of 6 tests are being used for validation purposes, H46–H49 as reference cases without step and H50 (FFS) and H52 (BFS) as stepped cases. A stepped case consists of two fuel grains with different initial port diameter that are connected together with a pin. For the FFS case of H50, for example, the initial port diameter in the first grain is 40 mm, whereas the second grain has a diameter of 25 mm, effectively forming a step of 7.5 mm. The same applies for the BFS of H52, however, for H52 the initial port diameter changes from 25 to 40 mm.

Table 4 lists the experimental data of the tests. We define the average oxidizer mass flux $G_{ox,avg}$ as:

$$G_{ox,avg} = \frac{16\dot{m}_{ox}}{\pi(D_0 + D_f)^2}, \quad (14)$$

with $D_f$ being the final diameter estimated with the mass loss method following the recommendation of Ref. [29]. In the stepped cases, the initial diameter is calculated with the arithmetic mean between the two grain halves, thus 32.5 mm. The average regression rate $\dot{r}$ can be calculated with the final diameter ($D_f$) and the burn time $t_b$ as:

$$\dot{r} = \frac{D_f - D_0}{2t_b}. \quad (15)$$

In order to be able to validate the numerical simulations with the experimental results, it is necessary to estimate the relative error of the regression rate for tests H46–H49. Following Eq. (15), the main sources of error for the regression rate are identified as:

1. Uncertainty in burn time $t_b$.

2. After the tests, the fuel grain is cut along its longitudinal axis to obtain the final diameter and, thus, the regression rate profiles. Errors during the cutting process of the post-firing grain need to be considered.

3. After the grains are cut, they are digitized using an office scanner and the final diameter calculated using image post processing. Resolution of the scanner and errors during the image processing have an impact on the obtained regression rate values

The error in the estimation of burn time arises from the method used to calculate the burn duration. In this work, the burn time is derived from the chamber pressure curve during the tests based on the approach of Durand et al. [33]. Using the quadratic module of chamber pressure time derivative:

$$\left| \frac{dP}{dt} \right|^2, \quad (16)$$

the start and the end of the start up and shut-off phases can be identified by the peaks of the derivative. Hence, the error in the burn time

estimation can be calculated using the uncertainty of the start ($\delta t_{start}$) and end phase ($\delta t_{end}$) using the propagation of errors as [33]:

$$\delta t_b = \sqrt{\delta t_{start}^2 + \delta t_{end}^2}. \quad (17)$$

The values for $\delta t_{start}$ and $\delta t_{end}$ are defined as difference between the first occurrence of a maximum peak and a later second peak with 1% of the value of the preceding peak [33].

For the error induced by cutting the grain in half and the image processing, the average burnt thickness ($\Delta r$) has to be calculated first for each case. Using the mass loss method given in Eq. (13) and the initial diameter $D_0$, the burnt thickness can be calculated:

$$\Delta r = \frac{D_f - D_0}{2}. \quad (18)$$

Next, the relative error due to cutting needs to be estimated. Because of possible parallel displacement of the cutting axis, the falsely measured final diameter after the cut ($D_{f,cut}$) can differ from the true final port diameter. The method proposed in [2] allows to estimate the false final diameter ($D_{f,cut}$) based on the cutting error of each test. In the following we assume a displacement of the saw of 1 mm. Hence, the burnt thickness obtained after the cut ($\Delta r_{cut}$) can be described as:

$$\Delta r_{cut} = \frac{D_{f,cut} - D_0}{2}. \quad (19)$$

Finally the relative error of the cutting method denotes to:

$$\delta_{cut} = \frac{|\Delta r_{cut} - \Delta r|}{\Delta r}. \quad (20)$$

As for the assessment of the error during image processing when digitizing the fuel grain halves, the resolution of the scan ($\text{res}_{scan}$) has to be taken into account:

$$\delta_{scan} = \frac{\text{res}_{scan} \cdot n_{pixel}}{\Delta r}, \quad (21)$$

with $n_{pixel}$ being the number of miscalculated pixels during image processing. The scan resolution $\text{res}_{scan}$ in all scans amounts to 0.044 mm/pixel, the number of falsely calculated pixel is set as a best guess to 5. This leads in all tests to an absolute error of 0.21 mm.

Summarizing all relative errors ($\delta$), the total error on the regression rate can be estimated as:

$$\delta_{total} = \sqrt{\delta t_b^2 + \delta_{cut}^2 + \delta_{scan}^2}. \quad (22)$$

The errors obtained for the tests H46-H49 are listed in Table 5. The cases of H46 and H48 with a higher initial mass flux show estimated total errors of around 3% while the cases of H47 and H49 with smaller initial mass fluxes display errors of up to 6%. To reflect a more conservative estimation ($\delta_{cons}$) which includes also errors induced by the assumption of axisymmetrical fuel consumption and other uncertainties that are hard to quantify, the errors for cases with small initial diameter (H46 & H48) are set to 5%, while cases of H47 and H49 are fixed to 10%. For the stepped cases H50 and H52 where the grain consists of different diameters in the first half and the second half, the error estimation is separated between the different halves. For example, in the case of H50 (FFS), the grain in the first half has an initial diameter of 40 mm (same as for H47 & H49) and the second half 25 mm (as H46 & H48). Thus, the conservative error estimation for H50 is 10% in the first section and 5% in the second section. For the BFS of H52, it is inverse. For all cases, the estimated error is indicated as a gray box in the respective plots in the next sections.
<table><caption>**Table 5**
Error estimation for regression rate profiles.</caption><thead><tr><th>Test</th><th>δt<sub>b</sub> [%]</th><th>δ<sub>cut</sub> [%]</th><th>δ<sub>scan</sub> [%]</th><th>δ<sub>total</sub> [%]</th><th>δ<sub>cons</sub> [%]</th></tr></thead><tbody><tr><td>H46</td><td>1.12</td><td>0.44</td><td>2.93</td><td>3.2</td><td>5.0</td></tr><tr><td>H47</td><td>0.82</td><td>0.60</td><td>5.89</td><td>6.0</td><td>10.0</td></tr><tr><td>H48</td><td>0.65</td><td>0.29</td><td>2.56</td><td>2.7</td><td>5.0</td></tr><tr><td>H49</td><td>0.77</td><td>0.51</td><td>5.20</td><td>5.3</td><td>10.0</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Validation of tests HYCAT46 & 48.</figcaption></figure>

### 4.3. Validation of cylindrical cases

Fig. 4 compares the experimental regression rate profiles (and the estimated error) of H46 and H48 with the numerical data obtained both with the flat profile and the varying mesh as explained in Section 4.1. To illustrate the length difference of H46 (250 mm) and H48 (500 mm) they are depicted in scale. Apart from the difference in length, both tests have the same initial fuel port diameter of 25 mm and therefore roughly the same average oxidizer mass flux. For this reason they are grouped here together.

For H46, the numerical regression rate profiles match the experimental profile surprisingly well, keeping in mind the assumptions and simplifications made for the meshes. Only the area of the higher regression rate due to the injector effect at the beginning of the grain (until around 0,075 m) is not covered by our numerical profiles. This is because the injector is not modeled in the simplified geometries in order to simplify the preparations and calculations. A decision that can be justified given the fact that the numerical profiles of H46 match the experiments well. Moreover, no considerable difference can be observed between the numerical flat and the numerical varying profile.

Turning the attention to H48, overall, the numerical fit is still satisfactory, even though less clear as for the H46 case. However, when comparing the experimental profiles of H46 and H48, the mismatch between 0.1 m and 0.25 m can be understood. H48 exhibits an area of higher regression rate directly after the injector effects (at approximately $x = 0.075$ m). H46 does not display this area. We suppose that instabilities (as can be noted when looking at the noisy experimental profile of H48 compared to H46) lead to this local increase in regression rate which is not covered by the numerical simulations. Nonetheless, the numerical profiles catch up to the experimental profile once it becomes less noisy, starting at around $x = 0.25$ m. Supposing the experimental profile of H48 would follow the trend of the test of H46, the numerical profiles, again, match sufficiently close. We will see when looking at H47, H49, H50 and H52 that this second peak after the injector effects is unique to H48. This indicates that the numerical profiles of H48 can be considered a satisfactory fit. Again, no discernible difference between the flat and the varying numerical profile can be observed.

In Fig. 5, the tests of H47 & H49 are compared to the simulations because they were both tested with the same oxidizer mass flux and

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Validation of tests HYCAT47 & 49.</figcaption></figure>

initial diameter of 40 mm. Looking at H47, the numerical and experimental trend are identical, however, the numerical profiles both overestimate the experimental data considerably (around 30%). To understand this behavior, it is necessary to take the absolute value of the regression rate into account. The tests of H47 and H49 were conducted with the initial diameter of 40 mm leading to relatively low mass flux levels of about 220 kg/m²s. Given the dependency of the regression rate on the mass flux, the absolute value of the regression rate for H47 & H49 (0.39–0.42 mm/s) is around half of the value of H46 & H48 (0.77–0.86 mm/s). Therefore, these tests are more susceptible to errors in cutting and overall post processing (see also Table 5). More importantly, the simulations are carried out on a single calculation with an average diameter – an assumption that becomes prone to relative errors at low mass fluxes given the small changes in diameter throughout the burn. Nonetheless, apart from a constant overestimation of the regression rate, the numerical profiles follow the same trend as the experiments.

Looking at H49, the numerical and experimental regression rates are close, with little deviation from the experimental data. Again, the varying numerical profile and the flat profile are almost identical, proving that using flat profiles on the average diameter is a valid approach allowing even simpler meshing approaches.

Summarizing, the flat cases without steps can be assumed validated, as has been shown in four tests. This allows using simplified geometries for the parametric study. Moreover, it has been shown that flat meshes and varying surface meshes show almost identical results. Thus, from here on we employ only flat surface meshes to maximize the collected data while limiting computational cost.

### 4.4. Validation of stepped cases

The approach to validate the stepped cases is as follows: first the flat numerical profile is compared to the experimental case (including error) at the correct oxidizer mass flow (0.343 kg/s for H50, red dotted line). The curves are depicted in Fig. 6. From here on, the step position will be at 0 mm axial distance. Clearly, after the step, the peak of the regression rate is underestimated and the numerical curve is increasing with further distance, whereas the experimental profile becomes horizontal. We suspect that the GSI model described in Section 3.2 overestimates the impact of the oxidizer mass flow on the regression rate when the flow is disturbed by obstacles that form a recirculation zone and, thus, increase the mixing. This could be due to mechanisms and interrelations in the determination of the regression rate through the GSI that are showing in these areas of increased mixing and turbulence. Moreover, radiation could play a non-negligible role at the zone directly after the step. Future investigations need to prove these assumptions.

To follow the suspicion that the impact of the mass flow after the step is falsely estimated in our GSI, we conducted additional simulations on lower mass flow levels (from the initial 0.343 kg/s to 0.2 kg/s
<figure><img src="image_9.png" /><figcaption>Fig. 6. Validation of test HYCAT50 (forward facing step).</figcaption></figure>

<table><caption>**Table 6**
Normalization coefficients for forward facing step case.</caption><thead><tr><th>Case</th><th>C<sub>norm</sub> [s/mm]</th></tr></thead><tbody><tr><td>Experimental</td><td>2.7</td></tr><tr><td>ṁ = 0.343 kg/s</td><td>2.3</td></tr><tr><td>ṁ = 0.2 kg/s</td><td>3.6</td></tr><tr><td>ṁ = 0.1 kg/s</td><td>5.9</td></tr></tbody></table>

<table><caption>**Table 7**
Normalization coefficients for backward facing step case.</caption><thead><tr><th>Case</th><th>C<sub>norm</sub> [s/mm]</th></tr></thead><tbody><tr><td>Experimental</td><td>1.2</td></tr><tr><td>ṁ = 0.346 kg/s</td><td>1.3</td></tr><tr><td>ṁ = 0.2 kg/s</td><td>2.1</td></tr><tr><td>ṁ = 0.1 kg/s</td><td>3.6</td></tr></tbody></table>

and 0.1 kg/s) which all lie in the bandwidth of historical HYCAT data [17]. To compare the regression rate profiles that are at different mass flows, we normalize the regression rate to the non-dimensional value of one just before the step. This normalization can be described by:

$$ \dot{r}_{\text{norm}} = \dot{r}_{\text{true}} \cdot C_{\text{norm}}, \qquad (23) $$

with $\dot{r}_{\text{norm}}$ being the regression rate after the normalization, $\dot{r}_{\text{true}}$ the real regression rate and $C_{\text{norm}}$ the coefficient used for normalization. The different normalization coefficients can be found in Table 6.

Looking at Fig. 6, it is clearly visible that the regression rates before the step are identical in shape after the normalization, justifying this approach. After the step, the initial suspicion of the mass flow dependence becomes obvious: with decreasing mass flow, hence, mass flux, the profiles predict the initial peak better. Moreover, at the 0.1 kg/s case (blue line with crosses), the numerical profile becomes also horizontal, as it is the case for the experimental curve. The reason why the position of the peak in the numerical simulations is more pronounced than the experimental is the fact that the step in the experiment moves with time as it is consumed also laterally. In our academic cases, we do not consider the progression of the step over time. Nonetheless, as it can be seen in Fig. 6, the simplified geometries at lower mass flow allow a satisfactory representation of the trends of the regression rate profiles with FFS, allowing them to be used in the parametric study.

Fig. 7 shows the same validation approach for H52 (BFS). The regression rates are normalized to the non-dimensional value of one before the step and the coefficients used for the normalization can be found in Table 7. Again, the numerical profiles before the step are almost identical in shape, showing the approach to normalize the regression rates from different mass flows is applicable for qualitative investigations.

Directly after the step, the 0.346 kg/s and 0.2 kg/s simulations underestimate the regression rate and follow a sloped profile after. The profile of 0.1 kg/s is able to correctly follow the trend directly

<figure><img src="image_10.png" /><figcaption>Fig. 7. Validation of test HYCAT52 (backward facing step).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Temperature contours of axisymmetric hybrid rocket simulation.</figcaption></figure>

after the step with some smaller underestimations and a steeper profile. Nonetheless, for qualitative analyses such as the parametric study of this article, the numerical profile at low mass flows yields suitable results.

Summarizing, we have shown the validity of the numerical model and the simplified geometries for a parametric study. For flat cases without steps, the profiles largely follow the experiments qualitatively and quantitatively. In the cases with steps (H50 & H52), the model loses some of its quantitative prediction, however, at lower mass flows the regression rates obtained with the simplified geometries are still well suited for a parametric study and will be applied from here on. Further effects such as the role of the turbulence and combustion models – as well as considering radiation – are to be investigated in the future to regain some quantitative prediction for stepped cases.

## 5. Results and discussion of the parametric study

In this section, the test matrix and the decision process for the boundary conditions is explained, before we discuss the outcome of the simulations. The results to be exploited are: (a) temporal and spatially averaged regression rates after the step, (b) temperature and regression rate profiles in axial direction at the fuel surface and (c) height-to-length relationship for the RZLs.

Before we continue, we present in short the basic characteristics of the flow over a BFS and FFS respectively. Fig. 8 shows the temperature fields and streamlines for a representative flow without step, a BFS with 5 mm height and a FFS with the same height. For the BFS case, the flame zone detaches after the step and re-attaches after the recirculation zone. The flow before the step is unaltered. For the FFS case, however, the step already lifts up the flame zone the closer it gets to the step and a first recirculation zone is formed before the step. After the step, the flow re-attaches considerably faster than for the BFS case. This is due to the fact that in the BFS case, the diameter increases and therefore the flow velocity decreases for a constant mass flow. In the case of the FFS, on the other hand, the diameter decreases and the flow accelerates. This fundamental difference has to be kept in mind when discussing the flow phenomena in the succeeding sections.
<table><caption>**Table 8**
Constant boundary conditions.</caption><thead><tr><th>T<sub>inlet</sub></th><th>P<sub>outlet</sub></th><th>Y<sub>H2O,inlet</sub></th><th>Y<sub>O2,inlet</sub></th><th>I<sub>turb</sub></th></tr></thead><tbody><tr><td>900 K</td><td>36 bar</td><td>0.59</td><td>0.41</td><td>11%</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Definition of reference case.</figcaption></figure>

## 5.1. Test matrix

The simulations will be carried out on the flat domain displayed in Fig. 1. The boundary conditions are being held close to the HYCAT data at lower mass flows as explained in Section 4. The grain length is set to 240 mm to limit the number of cells. Using only a single step at this stage allows to single out the effects of the steps. The interactions between multiple steps will be considered in future work.

The following major investigations are carried out:

### 1. Height study:

The impact of different step heights $h$.

### 2. Position study:

The effect of the axial position of the step ($x/l$).

### 3. Mass flow study:

The importance of the oxidizer mass flux $G_{ox} = \frac{\dot{m}_{ox}}{A}$ given by the oxidizer mass flow $\dot{m}_{ox}$.

These studies are carried out both for the BFS and FFS cases. Between the three different parametric campaigns, all other values are being held constant. The boundary conditions (based on typical HYCAT findings) that are similar throughout all calculations are listed in Table 8, the different inlet conditions, step heights and positions of the test matrix can be found in Table 9. The results of all studies are collected and their impact on the average regression rate, recirculation zone and temperature and regression rate profiles assessed in the respective subsections.

## 5.2. Definition of the reference cases

Keeping in mind that the diameter of the fuel port increases after the step, we need to consider a different reference case after the step than before for proper comparison. For example, with a mass flow of 100 g/s in a 5 mm FFS case, the reference $G_{ox}$ before the step yields 104 kg/m²s. After the step, it increases to around 203 kg/m²s due to the smaller diameter. Given the strong dependence of the regression rate on the mass flux, the reference needs to compare a different mass flux before and after the step. This is why the reference case is composed of two sub-references as displayed in Fig. 9 with a FFS as example.

## 5.3. Backward facing step cases

In the following sections, the average regression rate behind the step, the temperature and regression rate profiles of the BFS cases are analyzed separately.

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Regression rate after step for backward facing step.</figcaption></figure>

### 5.3.1. Regression rate behind step for backward facing step

Fig. 10 depicts the spatially and temporal averaged regression rate values for the BFS cases limited to behind the step. It allows to evaluate the impact of all parametric studies in one single image. The points show the reference cases and the crosses indicate the parametric cases. The label S2 i.e., corresponds to a step height of 2 mm, the label S5 end refers to the case of 5 mm step height at the position $x/l = 0.75$. Additionally, the Marxman [34] fit is given as:

$$ \dot{r} = a \cdot G_{ox}^{n} \qquad (24) $$

where $a$ and $n$ are parameters used to best fit the data points. For the Marxman fit, only the data points of the mass flow study are considered. This allows to extrapolate approximated values even outside of ranges that have been simulated. The mass flux $G_{ox}$ is calculated using the diameter after the step ($D_2$):

$$ G_{ox} = \frac{4\dot{m}_{ox}}{\pi D_2^2} \qquad (25) $$

Comparing the Marxman fits of the reference cases and the BFS cases, the most prominent observation is that the average regression rate of the BFS cases is inferior at mass fluxes above approximately 100 kg/m²s. Below this value, the average regression rate behind the step is higher than the reference cases without step. Moreover, the cases S1, S2, and S5 end are below the reference cases; whereas the cases S10 and S5 front are above the reference cases. This observation illustrates the existence of two competing effects: directly after the BFS, the flame zone is further from the fuel surface (see also Fig. 8), therefore decreasing the fuel surface temperature. On the other hand, the step introduces a recirculation zone that promotes mixing and heat transfer. If the recirculation zone is too small (as it is the case for 1 mm and 2 mm for higher mass fluxes), the regression rate after the steps never recovers from the flame zone being further from the fuel surface directly after the step. This is also why the case of S5 front, where the step is positioned at the beginning of the fuel port rather than the middle, is superior to the reference cases. The turbulence induced by the step has more time and space to compensate the decrease of the surface temperature directly after the step. In order to discuss this effect we introduce the temperature and regression rate surface profiles of the different cases in the next part.

### 5.3.2. Surface profiles for backward facing step

Fig. 11 shows the temperature distribution along the fuel port surface. Moreover, the gray boxes indicate the length of the recirculation zone induced by the step. The green dotted line specifies the profile of the respective reference cases. Keeping in mind that the diameter of the fuel port increases after the step, we need to consider a different reference case after the step than before as explained in Section 5.2. The step position is at 0 m distance.

With the help of these surface profiles, the aforementioned competing effects of recirculation zone and flame distance from the surface
<table><caption>**Table 9**
Values for parametric studies.</caption><thead><tr><th>Study</th><th>h [mm]</th><th>x/l [-]</th><th>m&#x0307;<sub>ox</sub> [g/s]</th><th>Labels</th></tr></thead><tbody><tr><td>Height</td><td>[1, 2, 5, 10]</td><td>0.5</td><td>100</td><td>[S1, S2, S10]</td></tr><tr><td>Position</td><td>5</td><td>[0.25, 0.5, 0.75]</td><td>100</td><td>[S5 front, S5 end]</td></tr><tr><td>Mass</td><td>5</td><td>0.5</td><td>[50, 100, 200]</td><td>None</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Temperature of the backward facing step height study.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Regression rate profile for different step positions. The step height is constant at 5 mm.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Regression rate of the backward facing step height study.</figcaption></figure>

become clearer. In all cases, the surface temperature directly after the step drops significantly, given that the flame is further from the surface. However, the surface temperature increases exactly until the re-attachment point (after the recirculation zone). From there on, the heat transfer and the mixing is increased, clearly visible in steeper temperature profiles as compared to the reference case. With the help of Fig. 11 it can be explained why the cases for S1 and S2 are inferior to the reference cases. The increased heat transfer induced by the step is too small to compensate the loss of surface temperature directly after the step. Starting with S5, the augmented heat transfer is considerably higher than the loss of temperature immediately after the step and can in total increase the surface temperature.

Fig. 12 displays the regression rate profiles of the same cases as Fig. 11.

The regression rate profiles follow the exact same trend as the temperature profiles. This is not surprising given that the surface temperature was modeled in the GSI (recall also Section 3.2) to be the

main influence of the regression rate. Knowing this, we will continue from here on to only consult the regression rate profiles for the sake of efficiency. The regression rate profiles confirm the competing effects: directly after the step, the regression rate descends well below the reference value. At the re-attachment point, the regression rate shows a peak (an effect also observed by Refs. [5,8,11]) and after a smaller rebound continues with a steeper profile than the reference cases (with the exception of the 10 mm case).

The effect of the sloped profile can be further illustrated with the results of the step position study in Fig. 13. Here, the step height is constant at 5 mm but the position moved from 25% (front) to 75% (end) of the total grain length. It becomes apparent that a BFS at 75% of the grain length diminishes the effect of the step.

Contrarily, for a BFS at 25%, the sloped regression rate profile considerably increases the regression rate for over the half of the total fuel grain. This explains also why the spatially averaged regression rate of the S5 front case in Fig. 10 is superior to the other S5 cases. Keeping in mind that the GSI model of CEDRE overestimates the mass flux dependence after the step (therefore showing a steeper profile) compared to the experiments. For the experiments (recall Fig. 7), the profile after the step is rather horizontal than sloped. However, even with an hypothetical horizontal profile in S5 front the case would be benefiting from a position upstream.

### 5.4. Forward facing step cases

Following the same procedure as for the BFS cases, this section discusses the outcome of the FFS simulations separated into average regression rate and local regression rate profiles.

#### 5.4.1. Regression rate for forward facing step

Fig. 14 reports the average regression rate behind the FFS in red and compares them to the reference cases without a step. The prominent observation is that – contrary to the BFS cases – all FFS cases exhibit a higher regression rate than the reference cases. However, the difference declines with increasing mass flux. Nonetheless, recalling Fig. 6, we observed the numerical regression rate profiles to be following a steep profile whereas for the experiment, the profile was rather horizontal. This leads to the averaged numerical regression rates being in all cases higher than the reference case. The reason why most of the FFS cases
<figure><img src="image_17.png" /><figcaption>Fig. 14. Regression rate after step for forward facing step.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Detail A of regression rate after forward facing step.</figcaption></figure>

are clustered at the same $G_{ox}$ value is that for the FFS domain, the diameter after the step is fixed to 25 mm, whereas the diameter before the step changes depending on the step height. That is why in Fig. 15, a more detailed view of the area labeled *Detail A* is provided.

Here, it is to be noted that the S1 case has about the same regression rate as the reference cases. Starting with S2, the difference grows. Interestingly, there seems to be a negligible difference between the S5 and S10 cases. A trend opposite to the BFS cases where the S10 BFS was by far superior to all other cases.

Moreover, for the BFS cases, the regression rate was highest when the step position is in the front. For the FFS, on the other hand, the step position at the end exhibits the highest regression rate. This observation will also be further investigated with the help of the surface profiles in the next section.

### 5.4.2. Surface profiles forward facing step

Before taking the surface profiles into account, it is reminded that we will concentrate solely on the regression rate profiles as explained in Section 5.3. Likewise to the BFS surface profiles, Fig. 16 compares the regression rate to the reference cases while indicating also the RZLs. Due to the diameter being smaller after the step (leading to increased $G_{ox}$) the reference cases before and after the step are different in order to refer to the proper mass flux as we elaborated in Section 5.2.

The first observation is that contrary to the BFS cases, the influence of the FFS is already noticeable before the step. This is due to the step ‘pushing’ the flame further away from the surface the closer it gets to the step. As soon as the recirculation zone before the step is reached, the regression rate drops drastically. This leads to an overall smaller regression rate before the step as compared to the reference cases. For the 10 mm case, this effect is most visible.

After the step, the regression rates show a steep increase with a distinct peak exactly at the re-attachment point. For the BFS cases, this peak was visible, but less pronounced. Moreover, for the BFS cases,

<figure><img src="image_19.png" /><figcaption>Fig. 16. Regression rate of the forward facing step height study.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. Regression rate profile for different step positions. The step height is constant at 5 mm.</figcaption></figure>

the recirculation zone is considerably larger than for the FFS cases. This translates to the FFS regression rate peak being reached earlier than for the BFS cases. After the peak, the regression rate follows approximately the same angle as the reference case. Interestingly, the difference between the 5 mm and 10 mm step is almost negligible after the step, whereas the negative effects before the step itself are considerably higher. This leads to the conclusion that the FFS need not be higher than 5 mm. Another fundamental differences to the findings of the BFS cases, where the 10 mm case is to be preferred.

If we turn now to the step position study that is detailed in Fig. 17, we can observe another discrepancy to the BFS cases. For the FFS cases, the step position further downstream increases the regression rate the most. At the 25% position, the peak is not only smaller but also less pronounced. Additionally, the decrease in regression rate before the step is noticeable already before the recirculation zone starts. Given that the turbulent boundary layer and its flame zone is less developed at the 25% position, we can assume that the disturbance of the flow at the earlier position is stronger.

### 5.5. Height-to-length relationship

In this section, we put the length of the recirculation zone (x) in relation to the height of the step (h) that is inducing the recirculation zone. This yields the height-to-length relationship $x/h$. In order to compare the different relations amongst each other, we introduce the step Reynolds number $Re_h$, following a common approach in the
<figure><img src="image_21.png" /><figcaption>Fig. 18. Length-to-height ratio of recirculation zones for backward facing step.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. Length-to-height ratio of recirculation zones before forward facing step.</figcaption></figure>

literature [35,36]:

$$Re_h = \frac{U \cdot \rho \cdot h}{\mu}, \qquad (26)$$

with U being the free stream velocity, ρ the flow density and μ the dynamic viscosity. All values are calculated at the domain inlet. Fig. 18 displays the RZL for the BFS cases. Additionally, a power law fit is plotted. The values of the BFS recirculation zones lie between 4.9 and 6.9 and their evolution is coherent with the trends that are found in (non-HRE) literature of BFS flows [37,38]: with larger Reh, the x/h ratio increases and shows saturating behavior. This can serve as one possible explanation why the BFS cases benefit from lower mass flows and, thus, lower Reh values. The local regression rate for BFS cases has shown to be inferior to the reference cases within the recirculation zone. It is typically after the recirculation zones that the regression rate overtakes the reference cases (refer to Fig. 12). Hence, with smaller recirculation zone relative to the step height, the BFS cases become more competitive. Given that the x/h decreases with decreasing mass flow, this could be an explanation for the trend in the BFS regression rates.

In Fig. 19, the x/h ratio before the FFS is plotted. The trend also is well in line with what is reported in (non-HRE) FFS literature [36]: with increasing mass flow and Reynolds number, the recirculation zone ratio before the FFS step decreases. Also the absolute values of the RZL ratio before the step are close to the experimental cases of Awasthi et al. [36].

Unexpectedly, the data for the x/h ratio after the FFS is inverse to what is found in (non-HRE) literature [36,39,40] as can be seen in Fig. 20.

Literature reports the x/h after the FFS to increase with the Reynolds number. However, for our FFS cases, the ratio decreases with increasing

<figure><img src="image_23.png" /><figcaption>Fig. 20. Length-to-height ratio of recirculation after forward facing step.</figcaption></figure>

Reynolds number. An explanation are the differences between the cases in the literature and our HRE related set-up. For example, in the FFS literature, usually the flow is not axisymmetric but on flat surfaces. The HRE FFS cases have an axisymmetric geometry, where the flow significantly accelerates after the FFS due to the smaller diameter. According to the Marxman law [34], the regression rate increases with increasing mass flux after the FFS. This leads to increased blowing from the surface, especially at the end of the recirculation zone where a pronounced peak of regression rate was observable. All these alterations to the typically planar literature case serve as explanation for the unexpected discrepancy of the (non-HRE) literature and our results.

## 6. Conclusion and outlook

In this work, we have conducted a parametric study on hybrid rocket flow-fields with single steps. We have shown that the use of simplified geometries yields satisfactory results as compared to the experimental validation cases. Moreover, a single simulation at the average between the initial and the final fuel port diameter is sufficient to predict the qualitative and – with minor limitations – the quantitative regression rate profiles. This allows to minimize the computational effort for large parametric studies and motivates simplifications to the usually complex computational domains and geometries employed in the field of hybrid rocket numerical simulations. Additionally, these calculations do not need to be conducted at different time steps.

The following main observations from the parametric study are to be recalled:

1. Both backward- and forward facing step have shown to be more effective with decreasing oxidizer mass flux.

2. In all cases, the peak in regression rate was observed to be at the re-attachment point.

3. For backward facing steps, the higher the step, the better the regression rate enhancement; for forward facing steps, the augmentation has shown to be saturated for 5 mm step height.

4. For backward facing step cases, the regression rate increases more after the recirculation zone than within, that is why a step at the beginning of the fuel port yielded better results than at the other positions downstream. This trend is considerably less pronounced for forward facing steps, and might even be inverse.

This work concentrated on single steps only. It is of importance to confirm the results for multiple steps and validate if the interaction of several steps is benefiting the regression rate increase or maybe even diminishing it. Preliminary experimental results on the same motor with multiple steps have shown that multiple steps increase the regression rate more than single steps, therefore hinting at considerable
regression rate increments possible with multi-stepped geometries (as it was also the case for Sakashi et al. [14] and Kumar and Joshi [15]).

Another concern is the evolution of the steps with time. If the steps dissolve rapidly as the burn continues, the turbulence enhancing effects would also get lost. However, in the experimental data (especially visible for the forward facing step of test H50) it has been shown that the steps diminish slower than expected as both surfaces – before and after the step – regress and the initial diameter differences before and after the steps are still clearly observable after our 10 s tests. Interestingly, for the forward facing step case, a noteworthy lateral progression of the steps is visible.

In conclusion, the results of the parametric study will serve as a starting point to develop an algorithm to distribute steps along an idealized profile as we proposed in Ref. [2]. We believe that backward- and forward facing steps can be a promising addition to hybrid rocket engines, especially when they are used to minimize manufacturing effort for complex geometries approximated with steps, while also increasing the performance.

## 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.

## Acknowledgments

The authors would like to acknowledge the support of Quentin Levard for his support with GMSH and CEDRE and we are thankful for the technical support of Jean-Charles Durand, Sébastien Grison, Cédric Clarmont and Christophe Corato during the experimental test campaign.

## Funding

The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 860956. It is part of the ASCenSIon project, an Innovative Training Network (ITN) to advance space access capabilities (https://ascension-itn.eu/).

## References

[1] C. Glaser, J. Hijlkema, J. Anthoine, Evaluation of regression rate enhancing concepts and techniques for hybrid rocket engines, Aerotec. Missili Spazio (2022) http://dx.doi.org/10.1007/s42496-022-00119-4.

[2] C. Glaser, R. Gelain, A.E.M. Bertoldi, Q. Levard, J. Hijlkema, J.Y. Lestrade, P. Hendrick, J. Anthoine, Experimental regression rate profiles of stepped fuel grains in hybrid rocket engines, Acta Astronaut. 204 (2023) 186–198, http://dx.doi.org/10.1016/j.actaastro.2022.12.045.

[3] C. Glaser, J. Hijlkema, J. Anthoine, Approximating idealized hybrid rocket fuel port geometries using steps, in: 8th Edition of the Space Propulsion Conference, Portugal, 2022.

[4] P. Duban, La fusée sonde LEX, Extrait de L'Aéronaut.ique Et L'Astronaut. (2) (1968).

[5] A. Gany, Y. Timnat, Parametric study of a hybrid rocket motor, Israel J. Technol. 10 (1972) 85–96.

[6] M. Grosse, Effect of a diaphragm on performance and fuel regression of a laboratory scale hybrid rocket motor using nitrous oxide and paraffin, in: 45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, American Institute of Aeronautics and Astronautics, Reston, Virginia, 2009, http://dx.doi.org/10.2514/6.2009-5113.

[7] N. Bellomo, M. Lazzarin, F. Barato, A. Bettella, D. Pavarin, M. Grosse, Investigation of effect of diaphragms on the efficiency of hybrid rockets, J. Propuls. Power 30 (1) (2014) 175–185, http://dx.doi.org/10.2514/1.b34908.

[8] S. Zhang, F. Hu, D. Wang, N. Patrick Okolo, W. Zhang, Numerical simulations on unsteady operation processes of N2O/HTPB hybrid rocket motor with/without diaphragm, Acta Astronaut. 136 (2017) 115–124, http://dx.doi.org/10.1016/j.actaastro.2017.03.005.

[9] C.P. Kumar, A. Kumar, Effect of diaphragms on regression rate in hybrid rocket motors, J. Propuls. Power 29 (3) (2013) 559–572, http://dx.doi.org/10.2514/1.b34671.

[10] M. Dinesh, R. Kumar, Utility of multiprotrusion as the performance enhancer in hybrid rocket motor, J. Propuls. Power 35 (5) (2019) 1005–1017, http://dx.doi.org/10.2514/1.B37491.

[11] P.A.O.G. Korting, H.F.R. Schöyer, Y.M. Timnat, Advanced hybrid rocket motor experiments, Acta Astronaut. 15 (2) (1987) 97–104, http://dx.doi.org/10.1016/0094-5765(87)90009-9.

[12] J. Lee, S. Rhee, J. Kim, H. Moon, Combustion instability mechanism in hybrid rocket motors with diaphragm, J. Propuls. Power (2021) 1–10, http://dx.doi.org/10.2514/1.B38492.

[13] L. Kamps, K. Sakurai, Y. Saito, H. Nagata, Comprehensive data reduction for N2O/HDPE hybrid rocket motor performance evaluation, Aerospace 6 (4) (2019) http://dx.doi.org/10.3390/aerospace6040045.

[14] H. Sakashi, Y. Saburo, H. Kousuke, S. Takashi, Effectiveness of concave-convex surface grain for hybrid rocket combustion, in: 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2012, http://dx.doi.org/10.2514/6.2012-4107.

[15] M. Kumar, P. Joshi, Regression rate study of cylindrical stepped fuel grain of hybrid rocket, Mater. Today: Proc. 4 (8) (2017) 8208–8218, http://dx.doi.org/10.1016/j.matpr.2017.07.163.

[16] A. Refloch, B. Courbet, A. Murrone, P. Villedieu, C. Laurent, P. Gilbank, J. Troyes, L. Tessé, G. Chaineray, J. Dargaud, E. Quémerais, F. Vuillot, CEDRE software, AerospaceLab J. (2011).

[17] J.-É. Durand, F. Raynaud, J.-Y. Lestrade, J. Anthoine, Turbulence modeling effects on fuel regression rate in hybrid rocket numerical simulations, J. Propuls. Power 35 (6) (2019) 1127–1142, http://dx.doi.org/10.2514/1.b37528.

[18] J. Messineo, J.-Y. Lestrade, J. Hijlkema, J. Anthoine, Hydrodynamic instabilities modeling in hybrid rocket engines, J. Spacecr. Rockets 56 (5) (2019) 1371–1382, http://dx.doi.org/10.2514/1.a34425.

[19] J.-E. Durand, Développement et Validation Expérimentale D'une Modélisation Numérique Pour La Simulation D'un Moteur-Fusée Hybride (Ph.D. thesis), Université de Toulouse, 2019, URL: https://hal.archives-ouvertes.fr/tel-02919827.

[20] A.G. Kushnirenko, L.I. Stamov, V.V. Tyurenkova, M.N. Smirnova, E.V. Mikhailchenko, Three-dimensional numerical modeling of a rocket engine with solid fuel, Acta Astronaut. 181 (2021) 544–551, http://dx.doi.org/10.1016/j.actaastro.2021.01.028.

[21] X. Meng, H. Tian, X. Chen, X. Jiang, P. Wang, T. Wei, C. Guobiao, Numerical simulation of combustion surface regression based on Butterworth filter in hybrid rocket motor, Acta Astronaut. (2022) http://dx.doi.org/10.1016/j.actaastro.2022.11.003.

[22] L. Gong, X. Chen, O. Musa, H. Yang, C. Zhou, Numerical and experimental investigation of the effect of geometry on combustion characteristics of solid-fuel ramjet, Acta Astronaut. 141 (2017) 110–122, http://dx.doi.org/10.1016/j.actaastro.2017.09.027.

[23] G.D. DiMartino, P. Malgieri, C. Carmicino, R. Savino, A simplified computational fluid-dynamic approach to the oxidizer injector design in hybrid rockets, Acta Astronaut. 129 (2016) 8–21, http://dx.doi.org/10.1016/j.actaastro.2016.08.026.

[24] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J. 32 (8) (1994) 1598–1605, http://dx.doi.org/10.2514/3.12149.

[25] D.C. Wilcox, Reassessment of the scale-determining equation for advanced turbulence models, AIAA J. 26 (11) (1988) 1299–1310, http://dx.doi.org/10.2514/3.10041.

[26] C.K. Westbrook, F.L. Dryer, Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames, Combust. Sci. Technol. 27 (1–2) (1981) 31–43, http://dx.doi.org/10.1080/00102208108946970.

[27] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, Internat. J. Numer. Methods Engrg. 79 (11) (2009) 1309–1331, http://dx.doi.org/10.1002/nme.2579.

[28] N.N. Smirnov, V.B. Betelin, V.F. Nikitin, L.I. Stamov, D.I. Altoukhov, Accumulation of errors in numerical simulations of chemically reacting gas dynamics, Acta Astronaut. 117 (2015) 338–355, http://dx.doi.org/10.1016/j.actaastro.2015.08.013.

[29] M.A. Karabeyoglu, B.J. Cantwell, G. Zilliac, Development of scalable space-time averaged regression rate expressions for hybrid rockets, J. Propuls. Power 23 (4) (2007) 737–747, http://dx.doi.org/10.2514/1.19226.

[30] M. Migliorino, G. Gubernari, D. Bianchi, F. Nasuti, D. Cardillo, F. Battista, Numerical Simulations of Fuel Shape Change in Paraffin-Oxygen Hybrid Rocket Engines, 2022, http://dx.doi.org/10.2514/6.2022-3563.

[31] C. Carmicino, A.R. Sorge, Influence of a conical axial injector on hybrid rocket performance, J. Propuls. Power 22 (5) (2006) 984–995, http://dx.doi.org/10.2514/1.19528.

[32] A. Karabeyoglu, G. Zilliac, B.J. Cantwell, S. DeZilwa, P. Castellucci, Scale-up tests of high regression rate paraffin-based hybrid rocket fuels, J. Propuls. Power 20 (6) (2004) 1037–1045, http://dx.doi.org/10.2514/1.3340.

[33] J.-E. Durand, J.-Y. Lestrade, J. Anthoine, Restitution methodology for space and time dependent solid-fuel port diameter evolution in hybrid rocket engines, Aerosp. Sci. Technol. 110 (2021) 106497, http://dx.doi.org/10.1016/j.ast.2021.106497.
[34] G. Marxman, M. Gilbert, Turbulent boundary layer combustion in the hybrid rocket, in: Symposium (International) on Combustion, vol. 9, (no. 1) Elsevier BV, 1963, pp. 371–383, http://dx.doi.org/10.1016/s0082-0784(63)80046-6.

[35] L. Chen, K. Asai, T. Nonomura, G. Xi, T. Liu, A review of backward-facing step (BFS) flow mechanisms, heat transfer and control, Therm. Sci. Eng. Prog. 6 (2018) 194–216, http://dx.doi.org/10.1016/j.tsep.2018.04.004.

[36] M. Awasthi, W.J. Devenport, S.A.L. Glegg, J.B. Forest, Pressure fluctuations produced by forward steps immersed in a turbulent boundary layer, J. Fluid Mech. 756 (2014) 384–421, http://dx.doi.org/10.1017/jfm.2014.405.

[37] P. Nadge, R. Govardhan, High Reynolds number flow over a backward-facing step: Structure of the mean separation bubble, Exp. Fluids 55 (1) (2014) http://dx.doi.org/10.1007/s00348-013-1657-5.

[38] R.W. Pitz, J.W. Daily, Combustion in a turbulent mixing layer formed at a rearward-facing step, AIAA J. 21 (11) (1983) 1565–1570, http://dx.doi.org/10.2514/3.8290.

[39] M. Sherry, D. Lo Jacono, J. Sheridan, An experimental investigation of the recirculation zone formed downstream of a forward facing step, J. Wind Eng. Ind. Aerodyn. 98 (12) (2010) 888–894, http://dx.doi.org/10.1016/j.jweia.2010.09.003.

[40] J.F. Largeau, V. Moriniere, Wall pressure fluctuations and topology in separated flows over a forward-facing step, Exp. Fluids 42 (1) (2006) 21, http://dx.doi.org/10.1007/s00348-006-0215-9.