Numerical Investigation of Combustion in a Gaseous Bipropellant Rocket Engine

1 Introduction

Liquid rocket propulsion remains a focal point of aerospace research, driven by the ongoing transition toward high-performance, reusable, and cost-effective transportation systems [1,2,3]. Despite its technological maturity, the drive for optimized performance and the adoption of “green” or high-efficiency propellant combinations maintain a high level of scientific interest [4,5]. In this context, the Oxygen–Methane (Ox/CH4) combination has emerged as a primary candidate for next-generation engines, offering a greater specific impulse compared to kerosene, reduced coking tendencies, and compatibility with In-Situ Resource Utilization (ISRU) strategies [6,7]. However, the internal environment of a liquid rocket engine (LRE) involves extreme conditions where high-pressure turbulent flows, complex chemical kinetics, and intense heat release are deeply coupled [8,9,10]. Modeling these phenomena remains one of the most significant challenges in the field. A vast body of literature has addressed both experimental and numerical aspects of liquid bipropellant propulsion, with a specific focus on Ox/CH4 systems due to their widespread industrial application [11,12,13,14]. Significant efforts have been dedicated to characterizing flame stabilization, injection dynamics, and wall heat loads, using both high-fidelity Large Eddy Simulations (LES) and more computationally efficient Reynolds-Averaged Navier-Stokes (RANS) approaches [15,16,17,18].

Despite the wealth of existing studies, the accurate prediction of reactive flows in the proximity of the injector face, where density gradients are steepest and mixing scales are most critical, remains a subject of active debate [19,20,21]. While full-scale testing provides the final validation, there is a persistent need for experimental and numerical data on scaled, cost-effective, laboratory-class engines [22]. These sub-scale applications are fundamental for isolating specific phenomenologies, such as the interaction between injector geometries and combustion efficiency, which are often obscured in complex, full-scale systems. In this context, a comprehensive research line involving both experimental and numerical activities has been initiated. Leveraging a consolidated expertise in hybrid and monopropellant systems [23], the laboratory is now extending its capabilities to liquid bipropellant systems, building upon a solid background of internal flow modeling and experimental characterization of reactive systems [24].

The present work introduces a numerical framework aimed at the analysis of a 200N-class O2/CH4 bipropellant rocket engine platform. The primary objectives are to perform a systematic validation of the modeling strategy against a well-known literature benchmark [25] and to employ this validated framework to investigate the impact of different injection patterns on the combustion field. By comparing several 3D configurations, the study aims to quantify how specific injectors geometric features influence propellant mixing and flame distribution, providing a predictive tool to guide the upcoming experimental test campaigns on the laboratory’s scaled hardware. The paper is organized as follows: Section 2 introduces the experimental facility and the engine design; Section 3 describes the numerical framework; Results are reported in Section 4, and final conclusions are summarized in Section 5.

2 Experimental Set-up and Operating Conditions

Thruster Configuration and Operating Conditions

The numerical investigation focuses on a 200N-class Oxygen–Methane bipropellant rocket engine platform designed at the Aerospace Propulsion Laboratory [26,27]. A sectional view of an exmaple of the engine assembly is shown in Fig. 1.

As seen in the assembly, the platform is composed of a combustion chamber and a converging-diverging nozzle designed as a single graphite component. This component is enclosed by a steel case with a closing ring for protection and integration with the injector part and the test bench. The injector assembly consists of three primary components, as illustrated in Fig. 1: the injector head, which collects the fuel; the injector plate, designed for fuel distribution; and the central oxidizer injector. The first two components (head and plate) were designed to be modular, ensuring adaptability with a wide range of central oxidizer injector configurations. Pressure transducer ports are located near the injector face and upstream of the nozzle to monitor chamber pressure levels and gradients, while ignition port is positioned close to the injector to ensure reliable flame initiation. The engine geometry and operating parameters were defined to represent a realistic low-thrust bipropellant configuration, suitable for experimental testing and CFD-based analysis. The engine is designed to operate under fuel-rich conditions in order to limit combustion temperature and thermal loads on the chamber walls [26]. Numerical simulations are carried out at representative nominal operating conditions derived from the engine design. These conditions are intended to match the target experimental operating point and are used as boundary conditions for the CFD calculations. Fig. 2 and Table 1 summarize the main geometric and operating parameters of the bipropellant platform adopted in the simulations.

Figure 1: 200N-class Oxygen–Methane rocket engine assembly [26].

Figure 2: Computationaldomain dimensions.

The nominal combustion chamber pressure is set to 15 bar, with total propellant mass flow rate of 83.70 g/s corresponding to a thrust level of approximately 200 N. The oxidizer-to-fuel mixture ratio (O/F) is fixed at 2.66, ensuring fuel-rich operation. Five injector configurations are investigated to assess the influence of the injector geometry on flow development and combustion behavior. All configurations adopt a ring-type fuel injection layout, with methane injected through 14 equally spaced 1.00 mm orifices arranged along a single circular ring on the injector plate with a radius of 13.2 mm, while the oxidizer injection strategy varies between five selected configurations: two different showerhead injectors and three single axial oxidizer orifices with different diameters. Table 2 summarizes the geometric parameters of the main oxidizer injectors.

Fig. 3 and Fig. 4 present the schematics of a representative oxidizer injector configuration for each of the two categories.

Figure 3: Showerhead #2 injector.

Figure 4: Single orifice #3 injector.

Considering these configurations, fully three-dimensional simulations are required to accurately represent the injector geometry and the resulting flow field within the combustion chamber. In particular, the Showerhead #2 configuration, as depicted in Fig. 3, features five axial injectors, while the remaining eight injectors are inclined at 45°, generating strong three-dimensional momentum components, asymmetric jet interactions, and enhanced transverse mixing. This injector arrangement induces complex azimuthal flow non-uniformities, jet-to-jet interactions, and recirculation structures that may cannot be captured by two-dimensional axisymmetric models. As a result, the three-dimensional approach is essential to correctly describe turbulence–combustion interaction, flame stabilization, and mixing processes. So, on the basis of the two-dimensional model developed and validated [26], three dimensional simulations were conducted and adopted as the reference framework for validation and detailed analysis.

3 Numerical Modeling and Methodology of Bipropellant Rocket Engine Combustion

3.1 Computational Domain

As seen in Fig. 1 on Section 2, the combustion chamber features a circular cross-section and an overall axisymmetric geometry. The computational domain extends from the injector faceplate to the nozzle throat and is fully three-dimensional to accurately capture injector geometry, azimuthal non-uniformities, and three-dimensional mixing and recirculation patterns. The analysis focuses exclusively on the combustion and mixing processes within the chamber, hence the divergent section of the nozzle is omitted from the computational domain. This choice is further supported by the adoption of a chemical equilibrium model, which provides a reliable description of the flow physics up to the throat; conversely, such an assumption would be physically inconsistent in the divergent section, where rapid expansion phenomena typically lead to non-equilibrium or frozen-flow conditions. The three-dimensional domain was discretized using a mesh primarily composed of tetrahedral elements, with prism layers applied along the walls of both the combustion chamber and the injector faceplates. Therefore, mesh refinement strategies were adapted accordingly to meet the specific resolution needs, to satisfy the near-wall resolution requirements of the turbulence models employed in this study. Particular attention was devoted to ensuring that the mesh resolution was sufficient to capture strong gradients of velocity, temperature and species concentration in the near-injector region and within recirculation zones, while maintaining numerical stability and solution convergence. The first cell height at solid walls was defined to achieve non-dimensional wall distance (y+) values consistently below 1. In Fig. 5, the distribution y+ is presented along the domain.

To ensure the generality of the results and facilitate comparison across different scales, the axial coordinate x has been normalized by the maximum length of the domain, L. This dimensionless representation provides a clearer overview of the mesh resolution quality relative to the entire device geometry, confirming that the viscous sublayer is appropriately resolved throughout the wall-bounded regions. The average value of y+ is 0.17, the maximum value is 0.97 (convergent section) along the chamber wall. This level of refinement is required to fully resolve the laminar sublayer, ensuring consistency with the low-Reynolds turbulence formulations adopted in this work and allowing for an accurate prediction of wall heat fluxes and skin friction. Mesh quality metrics, including skewness and aspect ratio, were also monitored to ensure adequate discretization accuracy and solver robustness. The discretization of the domain is depicted in Fig. 6, showing the overall mesh structure along with the specific refinement of the inflation layers near the convergent wall to accurately resolve the boundary layer.

Figure 5: Wall y+ values along the computational domain.

Figure 6: Full section of 3D mesh with a detail on the inflation near the convergent wall.

The prism layers used provide near-wall refinement necessary to accurately resolve boundary layer behavior and capture viscous effects near the walls. The mesh comprises approximately 3 million cells. The fixed domain explicitly models the full injectors geometry, including the spatial distribution and orientation of both oxidizer and fuel injection orifices, thereby allowing a more realistic representation of the injection and mixing processes. This approach enables the capture of three-dimensional effects such as jet–jet interactions and azimuthal flow non-uniformities, which are inherently neglected in the two-dimensional axisymmetric formulation.

Grid Sensitivity Analysis

Starting from the reference mesh, described in Section 3.1, two additional grids were generated for the grid sensitivity analysis: a coarse mesh of about 2 million cells and a finer mesh of about 6 million cells. The refinement strategy was applied consistently over the entire computational domain, preserving the same mesh topology and maintaining similar refinement patterns in critical regions such as the injector outlets, near-wall zones, and recirculation regions. The meshes were generated using a uniform global sizing approach, with a global element scale factor equal to 1.0 and a maximum element size set to 1.0. In order to better capture geometric details and strong flow gradients, curvature- and proximity-based refinement were enabled, with a minimum element size limited to 0.5. This allowed for automatic local refinement in regions characterized by high curvature or small geometric features, particularly in the injector region and in the vicinity of solid boundaries. Near-wall resolution was ensured through the introduction of prism layers on all solid surfaces. The prism mesh was constructed using an exponential growth law, with an initial layer height of 4 × 10−4 m and a growth ratio of 1.2. A total of 30 prism layers were employed, corresponding to an overall boundary layer thickness of approximately 0.5 mm. The same near-wall treatment was preserved across all grid levels to ensure consistency in the comparison. In particular, the prism layer structure and first cell height were chosen to maintain y+ values around 1 for all meshes, ensuring adequate resolution of the viscous sublayer. To assess the influence of spatial discretization, key integral flow quantities were monitored and compared across the three meshes. In particular, the chamber-averaged temperature and chamber pressure were analyzed. Fig. 7 shows a log-log plot of the numerical error versus the number of cells for the average values of the two quantities along a line positioned in the center of the post combustion chamber, taking into consideration an exemplary case. The discretization error was estimated by evaluating the relative variation of these quantities between successive mesh refinements comparing simulation outputs with asymptotic ones calculated using Richardson’s extrapolation. Pas=Pfine+Pfine−Prefng−1(1) where, in this case n = 1.5 and g is evaluated as:

Figure 7: Numerical error vs. number of grid cells.

Numerical results in terms of chamber temperature and pressure are summarized in Table 3. Based on these considerations, the reference mesh of approximately 3 million cells was selected as an optimal compromise between computational cost and solution accuracy, and it was therefore adopted for all the simulations presented in this work.

3.2 Model Set-Up and Boundary Conditions

The numerical simulations were carried out using ANSYS Fluent, adopting a steady-state RANS approach. As in the two-dimensional axisymmetric reference configuration [26] in the three-dimensional simulations, the injector entries were treated as mass-flow inlets. The propellant temperature at the inlets was fixed at 300 K. Turbulence at the inlets was prescribed using the intensity and hydraulic diameter method. The combustion chamber and injector walls were modelled as adiabatic, no-slip surfaces. A gauge pressure of 1.01 bar was imposed at the outlet. Pressure operating condition was set to 0 bar. Species mass fractions and temperature fields were initialized uniformly according to the ambient injection conditions. Fig. 8 shows the boundary condition on the simulation domain.

Figure 8: Schematic Representation of Boundary Conditions on the Simulation Domain.

Combustion was modeled using either the chemical equilibrium or the flamelet approach for non-premixed reacting flows. These methods were evaluated in combination with the standard k-ε or the SST k-ω turbulence models. Comparisons were performed on the three-dimensional reference case, and the most accurate formulation (SST k-ω + Flamelet) was subsequently applied to all simulations. A summary of the numerical setup is reported in Table 4.

3.3 Governing Equations and Physical Assumptions

This section presents the methodological framework adopted for the numerical simulations. Accurate prediction requires robust turbulence and combustion models. Turbulence and combustion models used in this study are introduced in Section 3.3.1.

3.3.1 RANS Turbulence Modelling

The turbulent flows in the rocket combustion chamber are modeled using the RANS equations, which result from time- or ensemble-averaging the Navier–Stokes system. The standard k-ε model solves transport equations for the turbulent kinetic energy k and its dissipation rate ε. The eddy viscosity is defined as: ν T = C μ k 2 ε . The transport equations for k and ε are [28,29]:

The standard values for the constants are Cμ = 0.09, Cε1 = 1.44, Cε2 = 1.92, σε = 1.3 [29]. The k-ω model defines the eddy viscosity as: ν T = C μ k ω , and solves the following transport equations:

The constant values are C ω 1 = 5 9 , C ω 2 = 5 6 , σω = σk = 2, Cμ = 0.09 [29].

The Shear Stress Transport (SST) model blends k-ω near the wall and k-ε in the free stream, limiting turbulent shear stress in strong vorticity regions. The turbulent viscosity is defined as follows: νT=minCμkω,Cμk2|Ω|F2(7) where F2 is the blending function that depends on k, ω, and Cμ. |Ω| is the vorticity magnitude [29].

In the present work, the SST k-ω model was preferred over standard k-ε formulations due to its better handling of flow separation and strong pressure gradients near the injector outlets. By blending the k-ω model in the near-wall region with the k-ε in the far-field, the SST version avoids the typical over-prediction of turbulence in regions with high acceleration or separation. Combined with a mesh resolution of y+ ≈ 1, this choice ensures a more accurate description of the shear layer development, as widely reported in literature for high-speed internal reacting flows where experimental benchmarks are often replaced by well-established turbulence models [30]. While the steady-state RANS approach provides a computationally efficient framework for the comparative analysis of different injector configurations, its inherent limitations in capturing transient phenomena must be acknowledged. In cases involving strong three-dimensional features, such as the inclined-hole injection in Showerhead #2, the time-averaging process of RANS might lead to the smoothing of unsteady structures, including vortex shedding, jet shear layer instabilities, and flame oscillations. Despite these limitations, the current approach is deemed suitable for identifying macroscopic trends in chamber pressure and temperature distributions, which are the primary focus of this comparative study. However, the present results should be interpreted as time-averaged representations of the flow field, and further high-fidelity LES studies could be targeted in future works to resolve the fine-scale temporal dynamics of the mixing process.

3.3.2 Non-Premixed Combustion Modelling

In a non-premixed combustion, fuel and oxidizer are injected separately and mix only within the combustion chamber. Two widely used approaches in this context are Chemical Equilibrium and Flamelet models. Turbulent fluctuations in species concentrations, temperature, and reaction rates are represented statistically using the Probability Density Function (PDF) approach, which forms the basis for coupling turbulence and chemistry. A key variable is the mixture fraction f, which quantifies the local mass fraction originating from the fuel stream, serving as the independent variable for both Chemical Equilibrium and Flamelet models [16,17]: f=Zi−Zi,OXZi,Fuel−Zi,OX(8) where Zi is the elemental mass fraction of element i. The scalar dissipation rate, χ, quantifies local mixing intensity: χ=2D|∇f|2(9) with D being a representative diffusion coefficient.

Under the assumption of equal diffusivities, the Favre-averaged transport equation for the mixture fraction and for the mixture fraction variance are [16]: ∂∂t(ρf¯)+∇·(ρV˜f¯)=∇·klCP+μtPrt∇f¯+Sm+Suser(10) ∂∂tρf′2¯+∇·ρV˜f′2¯=∇·klCP+μtPrt∇f′2¯+Cgμt(∇f¯)2−Cdρϵkf′2¯(11) where p is the density, kl is the laminar thermal conductivity, Cp is the specific heat, μt is the turbulent viscosity, Sm represents source terms from droplet vaporization, and Cg, Cd are model constants [16].

The Chemical Equilibrium model assumes instantaneous reactions relative to mixing timescales. The local thermochemical state (temperature, species mass fractions, density) is uniquely determined by the mixture fraction, pressure, and enthalpy. Turbulence–chemistry interactions are incorporated via a presumed PDF. This model is computationally efficient and suitable for steady-state conditions where chemical timescales are much shorter than mixing timescales [31]. The Flamelet model represents turbulent non-premixed flames as an ensemble of laminar flamelets embedded in a turbulent flow, each parameterized by the mixture fraction, f, and scalar dissipation rate, χ. Precomputed flamelet tables are used to reconstruct the local thermochemical state during CFD simulation [32]. Turbulence–chemistry interactions are accounted for by integrating flamelet states over a presumed joint PDF. Instead of solving the species transport in the physical coordinate system, the flamelet equations for mass fractions (Yi) and temperature (T) are mapped onto the mixture fraction (f) space, as follows: ρ∂Yi∂t=12ρχ∂2Yi∂f2+Si(12) ρ∂T∂t=12ρχ∂2T∂f2−1CP∑ihiSi+ρχ2CP∂CP∂f+∑iCP,i∂Yi∂f∂T∂f(13) where ℎi and Cp,i are the enthalpy and specific heat of species i, respectively, and Si is the species source term [32]. The PDF framework models turbulence–chemistry interactions statistically by considering the joint PDF of mixture fraction, temperature, and species mass fractions: p(f)Δf=limT→∞1T∑iτi(14) where p(f) is the probability density function, T is the total sampling time, and ????i represents the i-th time interval during which the mixture fraction f remains within the range [f, f + Δf].

In this study, the GRI-Mech 3.0 chemical mechanism is employed to generate pre-tabulated flamelet libraries for Oxygen–Methane combustion, enabling accurate reconstruction of species profiles and heat release without solving stiff chemical ODEs during CFD simulations [33].

4 Results and Discussion

4.1 Model Validation

The numerical framework, described in Section 3, was assessed through comparison with the experimental and numerical results published by Keller and Gerlinger [25]. To this end, the literature case study was reproduced, which examines a single-injector gaseous oxygen–methane combustion chamber operating at a nominal pressure of 20 bar. The benchmark propulsion system employs a shear-coaxial injector, in which oxygen is injected through a central duct, while methane flows through an annular outer passage. This configuration was selected as a validation benchmark due to its close resemblance to the thruster analyzed in the present study. The experimental combustor in the literature benchmark features a rectangular cross-section. However, the corresponding numerical study, and consequently the present work in order to replicate the literature results, utilizes a 3D cylindrical domain, maintaining the same fundamental area ratios and injection conditions. The computational domain includes the entire combustion chamber and nozzle length, with symmetry conditions applied along two orthogonal planes, allowing only one-quarter of the domain to be modeled while accurately capturing the three-dimensional flow structures. The proposed numerical framework relies on the SST k-ω turbulence model coupled with a non-adiabatic flamelet combustion formulation. This specific combination was selected for its superior reliability in capturing the complex turbulence–chemistry interactions and near-wall gradients typical of confined rocket combustors, as previously demonstrated in [26]. The validation focuses on a qualitative evaluation of the flame morphology and thermal field, alongside a quantitative assessment of axial pressure evolution and wall heat flux distributions. Fig. 9 compares the temperature fields obtained from the reference case reported in the literature with those predicted by the proposed SST k − ω turbulence model coupled with the flamelet combustion approach.

Figure 9: Temperature contour comparison between the literature case [25] and the validation case.

The slight discrepancies observed in the radial temperature distribution compared to the benchmark study are primarily due to the different combustion modeling strategies. While the reference work employs a Finite Rate Chemistry (FRC) approach, the present study utilizes a Flamelet-based model to maintain computational sustainability for the subsequent 3D multi-injector simulations. Although the Flamelet approach assumes a faster chemical response to mixing scales compared to FRC, it successfully captures the peak temperatures and the overall axial flame development by incorporating detailed kinetics through the GRI-Mech 3.0 mechanism. This confirms that the selected framework represents a robust trade-off between chemical accuracy and the computational efficiency required for comprehensive engine design studies. Quantitative validation was conducted by comparing axial pressure profiles along the combustion chamber and wall heat flux. Fig. 10 illustrates the axial pressure distribution along the chamber, highlighting the difference between the proposed model and the literature simulations. Considering that non-adiabatic wall conditions are imposed with a fixed wall temperature of 300 K, the resulting pressure profile shows a good agreement with the trend of the experimental data and closely matches the pressure distribution obtained in the 3D reference simulations. Despite its lower computational cost, the adopted flamelet approach captures the essential features of the combustion process with high fidelity. Interestingly, this model shows a closer agreement with experimental data than the detailed finite-rate mechanism, with 21 species and 98 reactions, used in the reference study. This suggests that, for the current domain configuration, the flamelet model effectively represents the turbulence–chemistry interactions that might be over-parameterized or miscalculated by more complex kinetic schemes. Furthermore, this highlights that an accurate description of the scalar dissipation rate and mixing scales, as provided by the flamelet approach, is more critical for this specific hardware than the inclusion of exhaustive elementary reactions.

Figure 10: Pressure trend comparison between the literature case [25] and the validation case.

The wall heat flux distribution along the combustion chamber is shown in Fig. 11, where numerical results obtained with the validated SST k-ω flamelet model under non-adiabatic wall conditions are compared with experimental measurements and literature simulation results.

Figure 11: Wall Heat Flux comparison between literature case [25] and validation case.

The good agreement observed in both magnitude and spatial trend provides additional support for the validity of the selected modeling framework. As seen in Fig. 11, the wall heat flux behavior computed with the proposed model shows a very good agreement with both experimental and literature simulations in the region near the faceplate and the first part of the chamber. Then, moving along the chamber and going through the nozzle, the output of the heat flux of the proposed model and the literature simulations are diverging, showing a slightly underestimation of the proposed model heat flux with respect to the experimental points. It is worth noting, however, that while the current approach tends toward a slight underestimation, the error remains significantly lower than that of the reference study, which exhibits a substantial overestimation of the heat flux. In conclusion, the validated model successfully reproduces the main features of the flame structure, flow organization, and pressure and heat flux trends reported in the reference study with a reduced computational cost, thereby confirming the robustness and suitability of the numerical framework for the analysis of the target Oxygen–Methane thruster.

4.2 Effects of the Injection Pattern on the Combustion

This chapter presents the results of the three-dimensional numerical simulations, using the validated SST k-ω turbulence model coupled with the flamelet combustion approach, on the five different injector configurations, described in Section 2, of the 200N-class Oxygen–Methane bipropellant rocket engine. Results in terms of combustion, temperature and pressure behavior, and overall performances are discussed. To facilitate a systematic comparison between the different injector designs, the results are organized into two primary categories based on their geometric configuration: ‘showerhead’ patterns, characterized by multiple distributed injection holes, and ‘single-hole’ configurations, featuring a single central hole. This categorization allows for a clearer analysis of how each fundamental injection strategy affects the combustion chamber’s performance. In Fig. 12 are shown the output of the simulations considering the two shower-head oxidixer injectors.

Figure 12: 200N engine three-dimensional CFD simulations: (a) Showerhead #1, (b) Showerhead #2. Temperature contours on the z-plane (top) and different x-planes (bottom).

As can be seen, Showerhead #1 configuration exhibits a relatively compact and nearly axisymmetric flame. The temperature field on the xy-plane displays a well-centered high-temperature core confined around the chamber axis, with a gradual radial spreading as the flow progresses downstream. The yz-plane cuts reveal that, at x = 10 mm, the hot region is still strongly localized, while further downstream the temperature distribution becomes smoother and more uniform, with limited azimuthal distortion. The associated recirculation pattern is predominantly toroidal and relatively weak, concentrated near the injector face and mainly driven by the shear between the oxidizer jets and the surrounding methane flow. In contrast, the Showerhead #2 configuration, featuring eight of the thirteen oxidizer orifices inclined at 45°, develops a more energetic and three-dimensionally structured flame. The xy-plane temperature contour in Fig. 12b displays a hotter and more elongated core, extending further downstream and occupying a larger portion of the chamber cross-section. The yz-plane sectors clearly highlight stronger azimuthal non-uniformities: already at x = 10 mm the hot regions are more asymmetric, and as x increases, the flame envelope becomes broader and more irregular, with pronounced high-temperature lobes. This behavior reflects the presence of a stronger and more complex recirculation system, where the inclined jets generate swirl-like motion and multiple interacting vortical structures that enhance entrainment and sustain combustion over a longer axial distance. The results of the three-dimensional simulations of the single-orifice injector cases are depicted in Fig. 13.

The contours show the xy-plane temperature contours and yz-plane cuts for the 2 mm, 4 mm, and 6 mm oxidizer orifices. In all three configurations, the oxidizer is injected through a single axial orifice, and the three-dimensional set-up explicitly resolves the discrete fuel injection pattern, consisting of fourteen methane jets arranged in a circular ring. This geometric detail introduces azimuthal modulation in the flow field and significantly affects the development of recirculation zones and flame symmetry. In Single Orifice #1, the 2 mm oxidizer jet produces a high-velocity core that dominates the axial flow and confines the combustion zone near the injector. In this configuration, the oxidizer injector operates under choked conditions, resulting in a strong throttling effect that fixes the mass flow rate and significantly increases the jet momentum at injection. The temperature field on the xy-plane shows a narrow and intense flame, while the yz-plane cuts reveal a compact high-temperature region with limited radial and azimuthal spreading. Recirculation is weak and localized, consistent with the strong jet momentum and minimal entrainment. In Single Orifice #2, the larger 4 mm oxidizer orifice reduces jet velocity and allows for more gradual mixing. The flame becomes broader and more symmetric downstream, with the contours of the yz-plane showing an increase in radial dispersion and the emergence of coherent recirculation structures. In Single Orifice #3 case, with a 6 mm oxidizer port, the xy-plane shows a wide and less intense high-temperature core, while the yz-plane cuts display azimuthal non-uniformities and well-developed recirculation zones. To quantitatively assess the impact of the injector geometry on the overall chamber performance, the main thermodynamic and performance indicators, namely chamber pressure, temperature, characteristic velocity, c*, and characteristic velocity efficiency, ηc*, were extracted at the inlet of the nozzle convergent for all five three-dimensional configurations. The characteristic velocity is defined as: c*=pcAtm˙(15) where pc is the pressure extracted at the nozzle convergent, At is the throat area, and m ˙ is the total mass flow rate. The efficiency is then calculated as: ηc*=c*cth*(16) where c t h * is the theoretical value obtained from chemical equilibrium analysis through Rocket Propulsion Analysis (RPA) software [34] at the same mixture ratio. These values are summarized in Table 5, which reports the results of the fully three-dimensional RANS simulations with discrete fuel and oxidizer injection.

Figure 13: 200N engine three-dimensional CFD simulations: (a) Single Orifice #1, (b) Single Orifice #2, (c) Single Orifice #3. Temperature contours on the xy-plane (top) and different yz-planes (bottom).

Among the configurations, Single Orifice #1 exhibits the highest overall performance, indicating that its injector design promotes strong mixing and effective heat release. However, this configuration operates with the oxidizer injector under choked flow conditions, implying that nominal operation requires high upstream pressure levels to sustain the imposed mass flow rate. This requirement may introduce practical limitations in terms of feed system complexity and overall engine operability, despite the favorable combustion performance. Showerhead #2 shows elevated temperatures and enhanced energetic flame characteristics, suggesting that the combination of axial and inclined oxidizer injectors generates pronounced three-dimensional recirculation and localized heat release. Showerhead #1 and Single Orifice #2 demonstrate intermediate performance, supporting stable combustion while producing moderate thermodynamic levels. Single Orifice #3 has the lowest overall performance, likely due to a slower jet and a more confined flame, which reduces heat release and characteristic velocity. To further quantify the performance shifts, Table 6 summarizes the percentage variations of the key variables relative to a baseline configuration, in this case Showerhead #1.

The quantitative results, shown in Table 6, highlight the fundamental trade-offs between the different injection strategies. The Single Orifice #1 configuration demonstrates the highest performance gain, with a 40.30% increase in combustion efficiency (ηc*) relative to the baseline. This performance boost is directly reflected in a 40.46% higher operating pressure compared to the baseline, highlighting that the superior mixing is driven by the extreme kinetic energy and shear intensity of the single-jet configuration. This gain in efficiency, however, comes at a higher pressure cost, underscoring a critical engineering trade-off. While the Single Orifice #1 configuration maximizes turbulent mixing, it also significantly increases the demands on the feed-system and structural requirements. In contrast, the Showerhead #2 presents a more balanced performance profile. It achieves a substantial 25.37% efficiency boost with a more moderate pressure increase (+26.27%), while simultaneously yielding the highest thermal gain, with an increasing of 16.20% in peak temperature. This suggests that the geometric vectoring of the 45° canted orifices optimizes the heat release distribution through macroscopic recirculation, rather than relying only on kinetic intensity. Finally, the negative variations observed for Single Orifice #2 and #3 confirm that reducing the injection velocity below a critical threshold leads to a decay in mixing quality, effectively making it less performant than the multi-orifice baseline.

4.3 Analysis of Turbulent Mixing Dynamics and Semi-Empirical Correlation

The interplay between turbulence production and mixing effectiveness is investigated by comparing the axial profiles of Turbulent Kinetic Energy (TKE) and the Mixture Fraction Standard Deviation (σf), as demonstrated in Fig. 14 and Fig. 15. In high-speed shear flows, the homogenization process is fundamentally governed by the turbulent energy transfer, where the kinetic energy from large-scale fluctuations is dissipated into smaller scales, significantly increasing the interfacial area between propellants. This mechanism accelerates the scalar dissipation rate, which is the primary driver for reducing the concentration gradients within the chamber. By correlating the spatial evolution of the turbulent field with the decay of σf, this analysis identifies how different injection architectures influence the flow field, either through high-velocity shear layers or vortex-induced entrainment, to achieve homogenization.

Figure 14: Comparison of Turbulent Kinetic Energy variation along the domain for different injection patterns.

Figure 15: Mixture Fraction Standard Deviation variation along the domain for different injection patterns.

The axial evolution of the Turbulent Kinetic Energy (TKE), illustrated in Fig. 14, highlights a fundamental difference in turbulence production mechanisms between the investigated designs. The Single Orifice #1 (green line) exhibits a distinctive behavior characterized by a sustained and nearly constant growth rate of TKE along the entire chamber length. Unlike the Showerhead configurations, which show a peak in turbulence intensity around x/L ≈ 0.3 followed by a decaying trend, the single-jet configuration maintains high shear-induced turbulence. This sustained production suggests that the high-momentum core of the Single Orifice #1 provides a continuous energy source for turbulent fluctuations, preventing the early dissipation observed in the multi-orifice patterns. Based on these considerations, the data reveals that the Single Orifice #1 configuration provides the most effective mixing environment. This configuration generates a continuous and superior growth of TKE along the chamber, reaching a peak of approximately 3424 m2/s2 at the exit plane (x/L = 1). This high level of turbulent intensity, driven by the significant momentum of the single axial jet, is the primary driver for the rapid decay of species non-homogeneity. Regarding the impact of the turbulent field on the mixing quality, it is quantified through the Mixture Fraction Standard Deviation (σf) in Fig. 15. All configurations initially exhibit a rise in σf up to x/L ≈ 0.3, representing the initial instability and breakup phase of the propellant streams. However, beyond this critical point, the Single Orifice #1 shows a superior mixing recovery, being the only configuration to exhibit a decisive and steep descent toward the minimum value at the chamber exit. While the Showerhead #2 manages to improve upon the baseline, the Single Orifice #1 clearly outperforms the Showerhead #1 (blue line), in which the standard deviation value remains significantly higher. This sharp divergence confirms that the sustained TKE production identified in Fig. 14 is effectively utilized to drive the rapid homogenization of the mixture. Consequently, Single Orifice #1 achieves the lowest σf values throughout the mixing process, concluding with a final value of 4.49 × 10−2 at x/L = 1. The enhanced performance of Showerhead #2 relative to Showerhead #1, despite the distributed flow, is attributed to the 45° canted orifices. This geometric configuration introduces a radial momentum component that promotes the formation of macroscopic vortex structures, facilitating species entrainment even at lower peak TKE levels compared to the Single Orifice #1. While the Showerhead #2 shows improved performance over the parallel Showerhead #1, it remains secondary to the strong mixing capability of the high-velocity single jet. This suggests that, for the investigated geometries, the shear-induced turbulence from a concentrated, high-momentum stream is more efficient at reducing σf than the distributed secondary flows provided by multi-orifice designs.

4.3.1 Physical Interpretation of Temperature Contours through Mixing Dynamics

The qualitative evidence provided by the temperature distributions in Fig. 12 and Fig. 13 is further supported and explained by the quantitative analysis of Turbulent Kinetic Energy (TKE) and Mixture Fraction Standard Deviation (σf). The morphological differences observed in the flame fronts are, in fact, the direct manifestation of the underlying turbulent mixing regimes. The concentrated and elongated flame core of the combustion simulation using Single Orifice #1 injector, depicted in Fig. 13, finds its physical justification in the sustained growth of TKE along the chamber. The high-speed shear layer prevents the early decay of turbulence, maintaining a high scalar dissipation rate that stretches the flame and drives the steep descent of σf toward the exit plane. The visual transition from cold core to reaction zone is thus the result of a strong kinetic mechanism that is less effective when using Single Orifice #2 and Single Orifice #3 injectors, where the lower injection velocities lead to wider and less intense thermal plumes. Similarly, the rapid radial expansion of the flame front in simulations with Showerhead #2 injector, represented in Fig. 12, is quantitatively explained by the introduction of the geometric enhancement factor Φ = 1.435. The non-axisymmetric thermal patterns observed in the yz-planes provide the visual counterpart to the macroscopic vortex structures generated by the 45° canted orifices. This structural flow manipulation promotes large-scale entrainment, allowing Showerhead #2 configuration to achieve homogenization levels comparable to the high-velocity Single Orifice #1. These considerations confirm that the observed temperature fields are strictly governed by the competition between shear-induced turbulence and vortex-driven transport. In summary, the detailed analysis of turbulent dynamics reveals two distinct flow mechanisms: a shear-driven sustained turbulence in the Single Orifice #1, which maximizes homogenization through high kinetic energy dissipation, and a vortex-induced convective mixing in the Showerhead #2, where geometric vectoring compensates for lower TKE levels. These mechanisms explain the performance trends observed in Table 6, clarifying how the injection strategy directly dictates the balance between pressure cost and combustion efficiency. To synthesize the physical link between the observed thermal fields and the turbulent quantities, it can be concluded that the detailed analysis of turbulent dynamics reveals two distinct flow mechanisms: a shear-driven sustained turbulence in the Single Orifice #1, which maximizes homogenization through high kinetic energy dissipation, and a vortex-induced convective mixing in the Showerhead #2, where geometric vectoring compensates for lower TKE levels. These mechanisms explain the performance trends observed in Table 6, clarifying how the injection strategy directly dictates the balance between pressure cost and combustion efficiency.

4.3.2 Semi-Empirical Correlation

To translate the numerical findings into a functional tool for rapid injector dynamics investigation, a semi-empirical correlation has been developed. This model scales the mixing quality at the chamber exit, quantified by the mixture fraction standard deviation (σf,out) at the final exit plane (x/L = 1), against the fundamental geometric and kinematic scales determined at the injection plane (x/L = 0). The choice of σf,out as the target metric, rather than the mass fraction of a single species, is dictated by its ability to represent the integral state of propellant homogenization regardless of local stoichiometric ratios. The structure of the proposed correlation follows the classical mixing framework established by Rupe [35], relating mixing effectiveness to the injection fluid dynamics and orifice geometry: σf,out=K·DeqL·Reeq−n·Φ(θeff)−1(17) where K is the global mixing constant, D e q = d h o l e n h o l e s is the equivalent hydraulic diameter, L is the combustion chamber length, n is the turbulent decay exponent, and Φ(θeff) is the mixing enhancement factor. This factor is defined as: Φ(θeff)=1+sinθeff(18) where θeff is the effective impingement angle, that for the single hole configurations and for the showerhead #1 injector is equal to 1. For the Showerhead #2 configuration, the directional influence is captured by the area-weighted effective angle that is computed as: θeff=arcsin(∑Aisinθi/∑Ai)(19) where Ai is the cross-sectional area of the i-th orifice, representing the weighting factor for the momentum distribution, and θi is the inclination angle of the i-th orifice axis relative to the longitudinal chamber axis (0° for purely axial holes and 45° for the canted peripheral holes in the Showerhead #2 configuration).

In this case, considering the configuration in which 8 of the 13 holes are inclined of 45°, the effective impingement angle is ≈25.8°, yielding Φ = 1.435.

The Reynolds number in Eq. (17) is computed as: Reeq=(ρO2uO2Deq)/μO2(20) where pO2, uO2 and μO2 are the oxygen density, velocity and laminar viscosity at the injector and Deq is the equivalent injector diameter.

The model parameters n and K are derived simultaneously through a logarithmic regression analysis of the purely axial configurations (Single Orifice #1, #2, #3), where Φ = 1. By linearizing the equation into the form:

The slope of the best-fit line identifies the exponent n ≈ 0.52, while the intercept determines the constant K. This value of n is physically consistent with the coaxial gas-jet mixing theory developed by Forstall and Shapiro [36], which demonstrates that scalar dissipation in fully turbulent regimes scales with the square root of the Reynolds number. Table 7 summarizes the validation of the correlation against the CFD datasets, showing a high degree of predictive accuracy.

The physical robustness of the model is highlighted by the comparison between the axial family and the canted designs. The Single Orifice #1 represents the ideal axial baseline: it achieves the highest mixing performance (σf,out = 0.0448) among axial cases solely through very high injection energy, leveraging an exceptionally high Reynolds number to compensate for the lack of radial transport. However, as the injection velocity decreases in Single Orifice #2 and Single Orifice #3, the mixing quality degrades significantly, confirming that axial injectors are strictly dependent on high pressure drops. In contrast, the Showerhead #2 architecture demonstrates a superior optimization of the mixing process. Despite operating at the lowest Reynolds number and possessing the largest equivalent diameter, Showerhead #2 achieves a homogenization level close to that of the high-velocity Single Orifice #1 and significantly better than the axial Showerhead #1. This proves that the introduction of radial momentum via the factor Φ effectively decouples mixing efficiency from extreme velocity requirements, making geometric vectoring a more robust strategy for high-performance combustion than pure turbulent kinetic energy.

5 Conclusions

In this work, a comprehensive numerical investigation of a 200 N-class Oxygen–Methane bipropellant rocket engine has been carried out with the objective of developing and validating a reliable CFD framework, with a reduced computational cost, for the analysis of injection and combustion processes in small-scale bipropellant rocket engines. Steady-state three-dimensional RANS simulations were performed using ANSYS Fluent, focusing on non-premixed gas–gas combustion under representative operating conditions. A detailed validation activity was conducted to assess the reliability of the proposed numerical framework. The model was used to replicate the single-injector GOx–GCH4 combustion chamber literature case [25]. By comparing the numerical results with the experimental and numerical data available in that study, the accuracy of the setup was verified. A combination of the SST k-ω turbulence model and a non-adiabatic flamelet combustion model is proposed. This approach accurately captures the reference flame structures while providing a high degree of correlation with experimental axial pressure and wall heat flux data. Once validated against the literature benchmark, the numerical framework was applied to the analysis of five injector configurations of the 200 N-class rocket engine presented in this work, including showerhead and single-orifice oxidizer layouts combined with discrete ring-type methane injection. Configurations featuring inclined oxidizer jets generate more complex three-dimensional flow structures and enhanced mixing, while single-orifice designs exhibit different flame stabilization and heat release patterns depending on jet momentum and orifice size. To maximize the effectiveness of the planned test bench activities, the performance parameters, including characteristic velocity and combustion efficiency, were used to discriminate among the different injector patterns. The results identify the configurations best suited to guarantee stable and efficient combustion for the target small-thrust operations. Overall, the present study establishes a validated and computationally efficient CFD methodology for the analysis of gas-gas bipropellant rocket engines, suitable for supporting ongoing and future experimental activities.