Evaluation of Large-Eddy Simulation Coupled with an Homogeneous Equilibrium Model for the Prediction of Coaxial Cryogenic Flames under Subcritical Conditions

Abstract

Large Eddy Simulations of liquid O${}_2$/gaseous H${}_2$ coaxial flames at subcritical pressure conditions are reported in this paper. These simulations reproduce the experimental Mascotte cases A1, A10 and A30, operating at 1, 10 and 30 bar, respectively, and for which temperature measurements and experimental visualisations are available. The main objective of this work is to assess the accuracy of the multi-fluid Homogeneous Equilibrium Model (HEM) described in Pelletier et al. (Computers & Fluids, 2020) for rocket engine applications. Of particular interest is the comparison with the experimental temperature measurements from Grisch et al. (Aerospace science and technology, 2003). To that purpose, numerical simulations are conducted with care, in order to ensure a proper statistical convergence and estimate the influence of the grid resolution for each case. Despite the crude assumptions—no surface tension and no atomisation model, for instance—that are made with the HEM used in this work, results are found to be in reasonable agreements with the measurements for case A10, even with the coarser grid. For case A30, a fine mesh resolution is required to capture the low intensity recirculation zone downstream of the inner jet necessary to reproduce the shape of the experimental profile. Finally, case A1 simulations, with the lowest Weber number, show large departures with the experimental measurements. This is expected to be due to a deficiency of the model to properly reproduce the two-phase dispersed flow.

1. Introduction

The jet formed during liquid injection is strongly dependent on the operating temperature and pressure conditions. When pressure and temperature are low enough with respect to the critical point, the fluid undergoes a classical break-up process. The interfaces correspond to discontinuities between the liquid and gas phases. As pressure or temperature are increased to supercritical conditions, surface tension vanishes and the phase discontinuity that is observed at lower pressure is no longer present. Instead, the jet evolves in the presence of a continuous interface between the high-density stream and the surrounding gases. Thermodynamic properties feature strong - but continuous - variations across a diffuse interface. Under such conditions, the fluid injection regime is often referred to as transcritical and the jet mixing is controlled by turbulence and is analogous to that of a variable density jet [1].

Supercritical and transcritical injection modeling has been deeply investigated during the last decades [2,3]. One important ingredient of such models is the appropriate description of the real-gas thermodynamics that drive the fluid behavior in such high pressure/low temperature situations. Thermodynamics generally rely on the use of a cubic equation of state (EoS), such as the Soave Redlich Kwong EoS [4]. In the context of Large-Eddy Simulations, good results have been obtained for the simulation of cryogenic coaxial flames, with satisfactory agreements in terms of flame topology between simulations and experiments [5,6,7].

In the subcritical domain, the strategies encountered in the literature for the modeling of compressible two-phase flows mostly focus on diffuse interface methods and more specifically on multi-fluid methods [8,9]. These methods seem well-adapted by construction to handle the interface appearance and disappearance that may occur in the targeted applications. They rely on an ensemble averaging of the phases properties. The resulting sets of equations are hyperbolic provided that convex thermodynamic closures are used [9]. Under the assumption of pressure, temperature and chemical potentials equilibria between phases, simplified multifluid models can be derived, that are convenient to treat numerically [10]. Such models have been recently coupled with cubic EoS and a multicomponent two-phase equilibrium solver in a classic finite-volume framework [11,12], showing very good results with available experimental data in the context of turbulent multi-component jet injection.

In the present contribution, coaxial cryogenic flames operating at subcritical pressures are simulated using the finite-element LES solver AVBP with the multi-fluid approach developed in [13]. The main objective is to assess the use of an homogeneous equilibrium model for the simulation of such flames. This is done comparing simulations with the temperature measurements performed by Grisch et al. [14]. The experimental cases of reference are first detailed in Section 2. The numerical setup is then given in Section 3. The sensitivity of the numerical results to the grid resolution is discussed in Section 4. The flame topologies are then analyzed in Section 5. Comparisons with the experimental measurements are eventually offered in Section 6. Finally, conclusions are presented in Section 7.

2. Experimental Reference Cases and Computational Domain

The experimental setup simulated in this paper corresponds to the single injector configuration of the Mascotte test-bench operated at Onera [15,16]. Low-velocity liquid oxygen surrounded by high-velocity gaseous hydrogen is injected through a coaxial injector in the chamber at 1, 10 and 30 bar. Oxygen is injected at 80 K, and hydrogen is injected at 289 K. Under such conditions, hydrogen is gaseous and oxygen is liquid, and the flow is expected to undergo a classical two-phase flow atomisation process. The injection conditions of interest have been experimentally studied in Grisch et al. [14] and Kendrick et al. [17]. They are detailed in

Table 1. Injection conditions used in the simulations. $\dot{m}$ is the mass flow rate at O${}_2$ and H${}_2$ inlets, $P_{ch}$ is the chamber pressure, and We = (${\rho }_{H2}$ u${}_{H_2}^2$ d)/$\sigma$, where ${\rho }_{H2}$, $\sigma$ and d are the hydrogen density at injection, the oxygen surface tension coefficient and the inner injector diameter, respectively.

Case	${\dot{\mathit{m}}}_{\mathit{H}2}$ [g/s]	${\mathit{u}}_{\mathit{H}2}$ [m/s]	P${}_{\mathbf{ch}}$ [bar]	We [-]	Experimental Data
A1	15	680	1	13,000	CARS [14], OH*-emission [18]
A10	23.7	300	10	28,000	CARS [14], OH*-emission [17,18]
A30	25.2	170	30	84,000	CARS [14]

. The mass flow rate of liquid oxygen is the same for all the cases and is set to 50 g/s, while the hydrogen mass flow rate is adjusted depending on the case, from 15 g/s for case A1 to 25 g/s for cases A10 and A30. The change in the operating conditions (injection velocity, surface tension trough pressure) leads to a large range of We numbers: from 13,000 for case A1 to 84,000 for case A30. Temperature measurements are available for cases A1, A10, and A30 [14], and OH*-emission images have been obtained in [17,18] for cases A1 and A10.

The computational domain consists of a 400 mm long rectangular chamber, with a 50 mm × 50 mm cross section. A longitudinal slice of the computational domain for case A10 is shown in Figure 1. The oxygen injector diameter is d = 5 mm, while the hydrogen annular injector diameter is d${}_{\mathrm{outer}}$ = 16 mm, 12 mm and 10 mm for cases A1, A10 and A30, respectively [14]. The outlet nozzle is not considered, and pressure is directly imposed at outlet through a non-reflecting boundary condition [19]. Mass flow rates, temperature and mass fractions are prescribed at inlet, also using characteristic boundary conditions [19]. Turbulent velocity perturbations are added at inlets following the method depicted in [7]. It relies on prescribing mean and fluctuating velocities profiles at inlet (given by Figure 9 in [7]). The most energetic turbulent length scale is set to half of the inner injector diameter or half of the annular injector width. The same profiles are used for all the cases and meshes considered in this work. Walls are considered adiabatic.

3. Numerical Setup

3.1. Governing Equations

Both Equation (3) and (4) models have been integrated in the AVBP solver (see [13] for details). The three-equation model, used in the simulations presented below, is detailed in this section. This model being similar to Euler equations used in gaseous and supercritical flows, similar closures are used here [7]. In the context of Large-Eddy Simulation, the corresponding Favre-filtered, fully compressible Navier–Stokes equations are given by:

$$\begin{array}{cc}\hfill \frac{\partial \overline{\rho }{\tilde{Y}}_k}{\partial t}+\frac{\partial \overline{\rho }{\tilde{Y}}_k{\tilde{u}}_j}{\partial x_j}& =-\frac{\partial \overline{J_{k,j}}}{\partial x_j}-\frac{\partial J_{k,j}^t}{\partial x_j}+\overline{{\dot{\omega }}_k}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}& =-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial \overline{{\tau }_{i,j}}}{\partial x_j}+\frac{\partial {\tau }_{i,j}^t}{\partial x_j}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial \overline{\rho }\tilde{E}}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_j\tilde{E}}{\partial x_j}& =-\frac{\partial \overline{p}{\tilde{u}}_j}{\partial x_j}+\frac{\partial {\tilde{u}}_i\overline{{\tau }_{i,j}}}{\partial x_j}-\frac{\partial \overline{q_j}}{\partial x_j}-\frac{\partial q_j^t}{\partial x_j}+\overline{{\dot{\omega }}_T}\hfill \end{array}$$

(3)

where $\overline{\varphi }$ and $\tilde{\varphi }$ denote spatial and mass-weighted (Favre) spatially filtered quantities. p is the pressure, T the temperature, $\rho$ the density, $Y_k$ is the mass fraction of the species k, $u_i$ represents the velocity vector components, $x_i$ the spatial coordinates, t is the time, E the total sensible energy, ${\tau }_{i,j}^t$ the sub-grid scale (SGS) stress tensor, $q_j^t$ the SGS energy fluxes, $J_{k,j}^t$ the SGS species fluxes, ${\dot{\omega }}_k$ the species reaction rate and ${\dot{\omega }}_T$ the heat release rate. The species ${\overrightarrow{J}}_k$ and heat $\overrightarrow{q}$ fluxes use classical gradient approaches. The fluid viscosity and the heat diffusion coefficient are calculated following Chung et al.’s method [20], and mass diffusion coefficients are deduced from heat diffusivity by assuming a unity Lewis number. Soret and Dufour effects are neglected. The sub-grid scale (SGS) energy and species fluxes are modeled using the gradient transport assumption with turbulent Prandtl and Schmidt numbers both set to 0.7. The SGS turbulent viscosity is modeled with the Wall-Adapting Large Eddy (WALE) model [21]. In the present work, combustion is modeled assuming infinitely fast reactions and pure diffusion regime operation [7].

3.2. Thermodynamic Closure

Mixture properties are computed using the Soave–Redlich–Kwong (SRK) EoS [4]:

$$p(\rho ,T,Y)=\frac{\rho rT}{1-b_m\rho }-\frac{a_m\left(T\right){\rho }^2}{1-b_m^2{\rho }^2}$$

(4)

with the mixture covolume $b_m$ and attractive coefficient $a_m$ computed following the van der Waals mixing laws [22]. In the subcritical domain, single-phase states can become unstable, leading to phase separation. The instability can be mechanical or chemical in the case of multicomponent mixtures, corresponding respectively to a loss of thermodynamic convexity along the pressure direction or along the chemical composition directions. Convexity can be restored by computing an equilibrium between a liquid and a vapor phase. In the present contribution, an approximation of the exact multi-component equilibrium is used assuming that, for each species, the liquid and vapor phase mass fractions are equal. Such an approximation is reasonable as long as phase separation is not dominated by a chemical instability. Practical methods for its calculation as well as more elaborated equilibrium closures are detailed in [13,23].

3.3. Meshing Strategy

An Adaptive Mesh Refinement (AMR) strategy is applied in this study. It relies on the use of mean solutions to locate poorly resolved regions of the flow. These regions are then locally refined using the MMG3D library [24]. For a given flow variable $\varphi$, the mesh scaling factor ${\Phi }^{*}$ is defined by [25]:

$$\begin{array}{c}\hfill {\Phi }^{*}\left(\varphi \right)=\left(1-\frac{\overbrace{{\Phi }^G\left(\varphi \right)}-{\Phi }_{min}^G\left(\varphi \right)}{{\Phi }_{max}^G\left(\varphi \right)-{\Phi }_{min}^G\left(\varphi \right)}\right)(1-ϵ)+ϵ\end{array}$$

(5)

where $ϵ$ determines the smallest scaling factor (sets to 0.5 here). $\overbrace{{\Phi }^G}$ indicates that the variable ${\Phi }^G$ is filtered using a Gaussian filter, applied 5 times. The criteria ${\Phi }^G$ are obtained using a Gaussian filtering of a mean flow variable $\varphi$. In this work, this filter is approximated by a second order derivative:

$$\begin{array}{c}\hfill {\Phi }^G\left(\varphi \right)=\parallel <\widehat{\varphi }>-<\varphi >\parallel \approx \frac{{\Delta }_x^2}{24}\parallel \nabla ·\nabla (<\varphi >)\parallel \end{array}$$

(6)

where $\varphi$ is a flow variable (i.e., such as velocity or sound speed for ex.), $\widehat{·}$ indicates Gaussian filtering, and ${\Delta }_x$ is the characteristic cell size. The mesh scaling factor $\Phi$ is computed using the average velocity, sound speed and heat release rate:

$$\begin{array}{c}\hfill \Phi =min\left[{\Phi }^G(u,v,w),{\Phi }^G\left(ss\right),{\Phi }^G\left({\dot{\omega }}_T\right)\right]\end{array}$$

(7)

The final mesh scaling factor is eventually obtained after a sequence of propagation of the smallest factor over five cells. This is to avoid any confinement of the flow by the grid. Different grid resolution could be obtained applying the method several times. The resulting grids are shown in Figure 2 and Figure 3, and their corresponding computational cost is summarized in

Table 2. Meshes. The number of refinement iterations corresponds to the number of times that the strategy depicted in Section 3.3 is used to refine the grid. Meshes and CPU costs given in this table correspond to case A10 as they are weakly sensitive to the case.

Mesh	Number of Refinement Iterations	Number of Nodes	CPU [kh] (for 10 ms)
M0	0	750,000	7.2
M1	1	3,300,000	92
M2	2	13,000,000	900

. Grids are essentially refined in the mixing layer near the injector exit, along the injector walls, at the position where the flow impacts the chamber walls and downstream of the inner jet where recirculations occur.

4. Influence of the Grid Resolution

The objective of this section is to investigate the sensitivity of the results to the grid resolution. The convergence of the temperature field is of particular interest to allow a proper comparison with the experimental temperature measurements. Ideally, mean results should be independent of the grid resolution so that proper conclusions may be drawn. Studies are carried out using longitudinal slices of mean temperature and radial profiles of temperature and axial velocity. Radial profiles are chosen at the experimental measurement axial positions.

The grids considered for each injection case are given in

Table 3. Grid sensitivity. Grids considered in the mesh convergence study for each injection point. ${\tau }_{av}$ is the averaging time.

Case	Mesh	${\mathit{\tau }}_{\mathbf{av}}$
	M0	36 ms
A1	M1	15 ms
	M2	9 ms
A10	M0	68 ms
	M1	33 ms
	M0	56 ms
A30	M1	39 ms
	M2	9 ms

with their averaging time ${\tau }_{av}$. Results for case A1, A10 and A30 are shown in Section 4.1–Section 4.3, respectively. Additional studies are provided in the Annex. First, the statistical convergence for each case on its final grid is given in Appendix A. It shows that the physical times used for transients and averages are sufficient to enable clear conclusions for each case. Finally, results obtained using a locally refined grid as depicted in Section 3.3 and a homogeneous grid refinement for case A30 are plotted in Appendix B. Virtually no differences are observed, which confirms that the local meshing strategy is well-adapted to the cases under study.

4.1. Case A1

Figure 4 shows temperature and axial velocity profiles for the three grid resolutions. While a reasonable grid convergence is obtained between meshes M1 and M2 at x = 80 mm, the grid resolution impact is high for the first profile at x = 40 mm. To understand this behavior, it is convenient to examine the flame topology. To that purpose, longitudinal slices of averaged temperature are shown in Figure 5. The flame length (qualitatively identified with the end of the blue (cold) region for example) strongly varies with the grid resolution. It is the shortest with mesh M0 and the longest with mesh M1. In addition, the flame initial expansion, near the injector, is shifted little downstream as the grid is refined explaining the strong sensitivity of the radial profiles at 40 mm with the grid resolution as a small error in this region leads to a strong modification of the profiles. This is not the case anymore further downstream where the flame radius is evolving more slowly with the axial position. A finer grid resolution is not affordable with the current resources available and thus the grid convergence cannot be rigorously demonstrated. However, it will be shown in Section 6.1 that the flow in this region is in good agreement with the experimental data, suggesting a sufficiently refined grid resolution. In addition, the profile at 80 mm for $r>15$ mm—the region where accurate measurements are available (see Section 6.1)—is weakly dependent on the grid resolution.

One explanation for the poor grid convergence is the strong coupling between the mean pressure distribution in the chamber and the annular injection velocity. The mean pressure is essentially driven by the pressure loss and the heat release distribution. Given the high compressibility of hydrogen under the thermodynamic conditions of the simulation, an error on the mean pressure at the injector exit directly translates in an error in hydrogen injection velocity, since the mass flow rate is conserved. This in turn modifies the entire flow fields as shear, and thus the initial flame spread are changed. The mean pressure in the first part of the chamber is plotted in Figure 6, where a departure of more than 20% is noticed between the meshes leading to a modification of the annular injection velocity as shown in Figure 7. This poor convergence could be due to the strong assumptions used to model the multiphase flow, especially the liquid vapor velocity equilibrium hypothesis that necessitates a very fine grid resolution to properly represent the huge shear between the two fluids. It is expected that out-of-equilibrium models [9] with proper relaxation models could improve the numerical predictions.

4.2. Case A10

Temperature profiles for meshes M0 and M1 are plotted in Figure 8. The two mesh refinements lead to close results, both in terms of mean and rms values. Similar results are obtained for axial velocity (not shown here). It is concluded that the results are nearly independent of the grid resolution for mesh M1.

4.3. Case A30

Temperature and velocity radial profiles are eventually plotted for case A30 (Figure 9) for mesh M0, M1 and M2. At x = 50 mm, meshes M0 and M1 show close results, but a sudden departure is observed with mesh M2, the latter featuring a larger initial expansion than M0 and M1. At x = 100 mm, a strong departure between the three grids is measured on the temperature, especially for r < 10 mm. Longitudinal slices of temperature in Figure 10 indicate that the flame length is reduced as the grid is refined. The initial opening of the flame is similar for the meshes M0 and M1 but increases as the mesh is further refined. Both the flame length reduction and the modification of the initial spreading rate are attributed to the formation of a large scale recirculation zone downstream of the inner jet with mesh M2 (Figure 11). One major consequence of this recirculation is a flattening of the temperature profile at 100 mm. As for case A1, a mesh convergence cannot be rigorously demonstrated here. However, the apparition of the recirculation seems necessary to explain the flat profile measured experimentally (see Section 6.3). This makes us think that mesh M2 may be necessary to capture properly the flow characteristics.

5. Flow Visualisations and Analysis

Before performing the comparison with the experimental results, it is first interesting to qualitatively analyse the flow topology for each simulated case. Results are presented in terms of instantaneous and mean flow fields in Section 5.1 and Section 5.2, respectively.

5.1. Instantaneous Fields

Longitudinal slices of temperature are shown in Figure 12. Regions of phase coexistence are also indicated with white iso-contours on the pictures. The flame topology strongly differs between the three flames. Case A1 features a very short dense core, violently shaken by the very high velocity annular stream. The two-phase flow region that follows is expected to actually correspond to an atomized spray that spreads down to one third of the chamber length. As pressure is increased, the liquid core length increases and the atomized region size diminishes. The flame length is also shorter. Finally, case A30 shows very small regions of phase coexistence, the oxygen being quickly mixed with the surrounding gases.

5.2. Average Fields

The average fields of temperature and oxygen mass fraction (Figure 13 top and middle) confirm the observation made with the instantaneous fields in Section 5.1: the flame length decreases as the pressure is increased. The oxygen penetration length is around 18d for case A30, 27d for case A10 and 35d for case A1. The mean axial velocity field indicates the presence of a large scale recirculation region for case A30 for 10d < x < 15d. This region is smaller for case A10 and no longer present for case A1. It is interesting to notice that, for each case, the first experimental profile (shown with a dashed line) is located just before the location where the flame interacts with the walls, at its maximum spreading rate. Thus, a small error on the flame expansion could strongly affect the results. The second profile is in the middle of the flame for cases A1 and A10, i.e., where the oxygen is still present and the flame is still opened. However, for case A30, the profile at 100 mm is already located at a position where the flame is closed. This distinct behavior is due to the large scale central recirculation region that develops for case A30. Finally, the last profile at 200 mm is located at the end of the flame for case A10, where all the oxygen has been burnt.

6. Comparison with Experiments

The flame topology extracted from LES is now compared with experimental results from [14,17]. Temperature measurements were obtained performing CARS thermometry on H${}_2$ and H${}_2$O molecules [14]. Each measurement is characterized by a validation rate, defined as the ratio between the number of spectra successfully processed and the total number of laser shots during a run. In addition, the resulting temperature may differ depending on the probed molecule. This situation is generally encountered when the validation rate is “low” (typically below 50 %). It is then decided in this work, for each measurement point, to keep the value with the highest validation rate. It is expected that the higher the validation rate, the more confidence one can have in the measurements. Finally, when comparing with the experimental visualisations, the objective is to assess if the mean flame position is properly retrieved. To this end, comparing OH* emission signal and OH field can be considered as an accurate approximation [18] and is used in this work.

6.1. Case A1

Temperature profiles are plotted in Figure 14. The mean profile at 40 mm indicates that the flame is qualitatively well-positioned. It seems that the maximum flame temperature is not captured by the experiments because of the limited discretisation. Given the low validation rate of the measurements at r = 0 and r = 5 mm, these values were discarded for the analysis. On the contrary, the other points correspond to a very high validation rate and are thus expected to be accurate. The mean temperature shows a reasonable agreement with the experimental measurements, even though a large departure of 500 K is present at r = 10 mm. As discussed in Section 4.1, this profile is highly sensitive to errors on the flame position. A small error on the flame spreading angle could lead to large errors on the radial temperature profile. The mean numerical results show large departures with experimental measurements for the profile at 80 mm. Given their low validation rates, it is difficult to conclude concerning the points at r = 0, r = 5 mm, and r = 10 mm. At r = 15 mm and r = 20 mm, errors are very large. One possible explanation is a wrong prediction of the atomisation and the droplet dynamics and evaporation, as no dedicated model is used here for the primary and secondary atomisations and the dispersed phase/flame interaction. Finally, the last profile at 400 mm shows better agreements between the simulation and the measurements, and the profile shape is reasonably recovered. However, there is an under-prediction of about 300 K at r = 0, 5 and 10 mm. The nearly flat numerical profile indicates an under-estimated flame length. This observation seems in agreement with the over-estimated temperature at 80 mm.

There is a large over-estimation of the rms values obtained by the simulations compared with the measurements at x = 40 mm and x = 80 mm. It is believed to be related to the spatial averaging induced by the experimental probe volume (1 mm-long and 50 $\mathsf{\mu }$m in diameter for H${}_2$ and 2 mm-long and 100 $\mathsf{\mu }$m in diameter for H${}_2$O [14]). It has a limited impact on mean profiles but could strongly decrease rms values if the flame thickness is smaller than the probe volume. This is expected to occur in the highly stratified and stretched regions of the flow, thus essentially near the injector. Indeed, further downstream at x = 400 mm, when the flow is more homogeneous, numerical rms values are in good agreement with measurements.

Figure 15 shows the radial position of the maximum OH* emission from experiments and OH mass fraction from the simulation. The flame shows a slightly under-predicted initial spreading rate. The axial position of the maximum flame radius is also positioned about 1d downstream of the measured location. Nevertheless, results are satisfactory, and the initial flame shape is reasonably retrieved. This is in agreement with the previous observation made on the temperature measurements at x = 40 mm.

In conclusion, the results for case A1 are mixed. The near-injector flame spreading angle and temperature are qualitatively captured, but large departures exist concerning the temperature at 80 mm between experiments and numerical simulations. These errors are attributed to a lack of models to properly represent the dispersed phase and droplet/flame interactions.

6.2. Cases A10

The temperature profiles for case A10 are plotted in Figure 16. Simulations are found to be in good agreement with the experiments for the profiles at x = 100 mm and x = 200 mm. These points feature high validation rates (greater than 70%). However, there is a large departure closer to the injector, at x = 50 mm. Unfortunately, given the low validation rates of these points, it is difficult to draw a clear conclusion on which data are erroneous. This near injector region can however be qualitatively validated thanks to the experimental visualisations. The comparison between simulation and experiment is shown in the left part of Figure 17. The flame shape is well-reproduced by the simulation. Its initial spreading rate is close to the experimental one. The axial position of its maximum opening is however slightly shifted in the simulation compared with the experiment. This observation is confirmed by comparing the radial position of the flame with the experimental measurement from [18] (Figure 17 right). The agreement is very good up to 9d. The maximum radius is larger in the simulation (3d) than in the experiment (2.5d). Contrary to the observations made with the temperature profiles, it seems from Figure 17 that the simulated flame is a little longer than the experimental one, a feature that was also observed in a previous work [26]. However, results are in general satisfactory for this case. These results are in line with those obtained by [27].

6.3. Case A30

The temperature profiles for case A30 are now compared with the experimental measurements in Figure 18. Because of its low validation rate, the point at r = 0 mm of the profiles at 50 mm is not accounted for in the discussion. The profile at 50 mm shows large departures with the reference, and the expansion of the flame is over-predicted. This is confirmed by adding in the comparison another profile taken two diameters upstream of the measurement location. A very good agreement is now obtained for this profile taken at 40 mm. The shape of the profile at 100 mm is well-retrieved. As discussed earlier, it is found that capturing the large scale recirculation was necessary to obtain this nearly flat profile. However, a large departure of 500 K is noticed. To investigate a possible impact of heat losses at the wall on the temperature, simulations are conducted assuming a constant wall temperature of 500 K. The results are compared with the adiabatic case in Figure 19 and show virtually no differences. Departure for r < 10 mm for the profile at 100 mm is attributed to a lack of statistical convergence as the isothermal simulation has been averaged with only 8 ms (after a transient time of 8 ms) for CPU cost reasons. At the moment, there is thus no clear explanation for the large temperature departure between simulations and experiments noticed at 100 mm.

7. Conclusions

Large-eddy simulations of three LOx/H${}_2$ Mascotte test cases operating at subcritical conditions are performed in this work. The numerical strategy relies on a homogeneous equilibrium model coupled with a cubic equation of state and a fast chemistry approach. The three simulated cases essentially differ by the chamber pressure (1, 10 and 30 bar) and the associated Weber number (13,000, 28,000 and 84,000, respectively).

Case A10 at 10 bar features a very good grid convergence, the coarser grid already giving the converged mean temperature profile. It also shows a good agreement with the experiments. On the contrary, cases A1 (1 bar) and A30 (30 bar) are strongly impacted by a refinement of the grid, and it was not possible to ensure a proper mesh convergence for these two cases. They also show limited agreements with the experimental temperature measurements. This suggests that the sub-grid models fail to properly reproduce the actual sub-grid contributions for these cases. However, errors do not seem to have the same origin for the two cases.

For case A30, results indicate the apparition of a large recirculation of the flow downstream of the inner jet for the finer grid. This behavior is in agreement with the experimental measurements featuring a nearly flat temperature profile at 100 mm, the footprint of an efficient mixing that is only observed numerically when the recirculation is present. While the shape of the profiles is qualitatively recovered, a large over-estimation of the mean temperature, for the profiles at 50 mm and 100 mm, is observed for the simulations. No explanations are found at this moment, since adding wall heat losses shows virtually no impact on the temperature in the region of interest.

The numerical results for case A1 show a very short intact core length for the liquid jet and a large region of atomized spray, up to one third of the chamber length. The poor convergence for this case is attributed to the strong assumptions used to model the multiphase flow, especially the liquid vapor velocity equilibrium hypothesis that requires a very fine grid resolution to properly represent the huge shear between the two fluids. The absence of dedicated models for the dispersed phase may also explain the limited agreements with the experimental measurements.

Current research focus on extending the present model to an out-of-equilibrium multi-fluid model [9]. In addition to the numerical challenges for solving such a system, closures are needed to properly represent the transfer between the phases as well as the two-phase turbulent combustion.