Subgrid Turbulent Flux Models for Large Eddy Simulations of Diffusion Flames in Space Propulsion

1. Introduction

Methane is an attractive option for the new generation of thrust chambers [1,2]. This fuel has promising capabilities due to its high specific impulse, easy production and handling [3]. However, several challenges related to ignition, mixing, combustion stability, and safe operation demand further research [4,5]. The design process of modern rocket engines relies on Computational Fluid Dynamics (CFD) to address these issues. Numerical simulations are a fast and inexpensive tool to investigate design modifications and optimize geometrical features. In addition, they can facilitate the understanding of test results since every relevant physical field is accessible. Consequently, CFD is an integral part of the design process of combustion chambers [6,7,8,9]. Traditionally, Reynolds Averaged Navier-Stokes (RANS) solvers have been extensively used in industry due to their low computational cost. In methane thrust chambers, multiple researchers have used this methodology for both the combustion chamber [10,11,12] and the cooling system [13,14,15,16,17]. The most relevant handicap of this approach is the suppression of transient behavior, which excessively restricts the degrees of freedom of the turbulent flow dynamics. Important features and physical processes are subsequently canceled, reducing the simulation fidelity. This oversimplification is particularly relevant in combustion chambers of rocket engines. In such a scenario, classical assumptions of turbulence are challenged due to the extreme conditions that take place. Turbulent diffusion flames for space propulsion applications may exhibit exotic phenomena such as reversed energy transfer [18], high anisotropy far from walls [19], and flame-generated turbulence [20,21,22,23]. Conventional RANS models are unfit to capture these processes due to their inherent steady nature. Scale-resolving methods are a promising alternative to approach the complexity of turbulent combustion in liquid rocket engines. With growing computational resources, the applicability window for Large Eddy Simulations (LES) has experienced a significant widening over the last years [24,25,26]. Despite their higher computational cost, this method circumvents most of the shortcomings associated with RANS by resolving the largest vortical structures. Under the current technological constraints, LES is the most accurate method in which the simulation of full-scale systems remains feasible [27] (although at a high computational cost). Due to its potential, LES for combustion applications has awakened considerable interest in the last few years.

The application of LES requires modeling unresolved turbulent quantities at a subgrid level. The extended approaches have been validated using empirical results or higher-fidelity Direct Numerical Simulations (DNS). However, the conditions for which most models were devised differ from the extreme environment found in combustion chambers of liquid rocket engines. In such applications, chemical scales are very thin and fast due to the high pressures and the usage of pure oxygen as oxidizer. These conditions have dramatic repercussions on the computational cost, which is driven by the necessary flame front resolution [28,29]. Artificially Thickened Flame (ATF) models [30,31] are a powerful alternative to address this challenge. In this approach, diffusivity and reaction rates are artificially varied to thicken the flame front while preserving flame speed. This method can provide invaluable computational cost savings enabling coarser grids. ATF has been applied to a wide variety of configurations, including non-premixed combustion [31], and particularly space propulsion engines [32]. A thickened flame alters the overlapping between vortical and chemical scales, altering the original interactions between turbulence and chemistry.When a flame is artificially thickened, the Damköhler and other dimensionless numbers are altered with respect to the original configuration. Such modifications have the potential of falsifying fundamental interaction such as the sensitivity of the flame displacement speed to stretching. These changes have a profound impact in the overall feedback between turbulence and combustion and can imparity the validity of the simulation if they are not carefully addressed. Efficiency functions are common solution used to mitigate the effect of neglected physical features [30,33,34]. Nonetheless, the performance and reliability of these models remains an unsolved challenge, which has motivated intensive research activity during the last years [35,36,37]. Other remarkable approaches for LES modeling in turbulent combustion applications include the G-equation model [38,39] and the Flame Surface Density (FSD) models [40].

The present paper focuses on the application of gradient models for closure of subgrid turbulent fluxes. This approach originates from the works of Clark et Al. [41] and Vreman et al. [42], which was later applied to combustion by Pierce and Moin [43]. Gradient models can determine subgrid turbulent fluxes assuming their similarity with the projection of the involved variables in the filtered space. This approach is particularly attractive if the flame dynamics are the dominating physical process at the filter level. Such a scenario is particularly interesting in high-pressure combustion applications, where the flame is very thin, and the unresolved eddies have low energy. In such situations, the influence of chemical processes to fluid dynamics is the main unresolved effect at a subgrid scale. Modeling of these processes without resorting to ATF allows for a partial resolution of the flame front, preserving the original physical flow properties and relevant dimensionless numbers (e.g. $Le$, and $Da$). The present work evaluates the applicability of this approach for turbulent flames in conditions representative of rocket combustors. The performance and the behavior of the model coefficients are assessed with databases obtained from high-fidelity numerical simulations. In addition, an adaptive strategy based on the local burning conditions is proposed to enhance the model's performance.

The rest of the paper is structured as follows. A brief theoretical introduction to LES and turbulent diffusion flames is provided in Section II. The methodology used to retrieve the subgrid turbulent dynamics is discussed in Section III, where the performed numerical simulations are detailed. The gradient model enhancement for subgrid turbulent fluxes is introduced and validated in Section 4. Finally, the results are summarized in Section 5.

2. Theoretical Background

The present section discusses basic concepts involving Large Eddy Simulations in diffusion flames. First, the mathematical foundations for the notion of LES filters are introduced. Secondly, the conservation equations of reacting flow subject to filtering are briefly commented on, emphasizing the role of unresolved subgrid fluxes and existing models. Finally, the most relevant characterization parameters in turbulent non-premixed combustion are presented to facilitate the discussion of the results.

2.1. Filters

In LES, the grid size $\mathsf{\Delta }$ is sufficiently small to resolve large eddies appropriately, but the smallest ones (i.e., the Kolmogorov eddies $\mathsf{\eta }$)are insufficiently resolved, as $\mathsf{\Delta }>(6/\mathsf{\pi })\mathsf{\eta }$ [44]. Therefore, unlike in DNS, it is not possible to assume a laminar behavior at subgrid levels, and statistical mechanics are relevant at these scales. The process of resolving a turbulent flow using a coarser resolution than necessary can be seen as applying a low-pass filter. In general, the application of a filter $\mathrm{F}$ on a quantity $\mathrm{q}$ can be defined as:

$$\overline{\mathrm{q}}=\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\mathrm{q}\left({\mathrm{x}}^{\prime },{\mathrm{y}}^{\prime },{\mathrm{z}}^{\prime }\right)\mathrm{F}\left(\mathrm{x}-{\mathrm{x}}^{\prime },\mathrm{y}-{\mathrm{y}}^{\prime },\mathrm{z}-{\mathrm{z}}^{\prime }\right)\mathrm{d}\mathrm{x}\prime \mathrm{d}\mathrm{y}\prime \mathrm{d}\mathrm{z}\prime$$

(1)

In variable-density flows, it is more convenient to perform the spatial integration weighing for density, as suggested by Favre [45]. This operation can be expressed as:

$$\tilde{\mathrm{q}}=\frac{1}{\overline{\mathsf{\rho }}}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\mathsf{\rho }\left({\mathrm{x}}^{\prime },{\mathrm{y}}^{\prime },{\mathrm{z}}^{\prime }\right)\mathrm{q}\left({\mathrm{x}}^{\prime },{\mathrm{y}}^{\prime },{\mathrm{z}}^{\prime }\right)\mathrm{F}\left(\mathrm{x}-{\mathrm{x}}^{\prime },\mathrm{y}-{\mathrm{y}}^{\prime },\mathrm{z}-{\mathrm{z}}^{\prime }\right)\mathrm{d}\mathrm{x}\prime \mathrm{d}\mathrm{y}\prime \mathrm{d}\mathrm{z}\prime$$

(2)

Which is denoted as the Favre-filtered value of $\mathrm{q}$. There is a wide variety of filters that can be applied [28,46,47]. In the frame of the present work, top-hat filters have been used because of their simplicity. These mathematical operators are explicitly defined as:

$$\mathrm{F}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)=\left\{\begin{array}{c}1/{\mathsf{\Delta }}^3, \mathrm{if}\max \left(\left|\mathrm{x}-{\mathrm{x}}_c\right|,\left|\mathrm{y}-{\mathrm{y}}_c\right|,\left|\mathrm{z}-{\mathrm{z}}_c\right|\right)\le \mathsf{\Delta }/2 \\ 0, \mathrm{otherwise}\end{array}\right.,$$

(3)

where $({\mathrm{x}}_{\mathrm{c}},{\mathrm{y}}_{\mathrm{c}},{\mathrm{z}}_{\mathrm{c}})$ are the spatial coordinates of the considered location. This filter corresponds to the average over a cubic box centered at $({\mathrm{x}}_{\mathrm{c}},{\mathrm{y}}_{\mathrm{c}},{\mathrm{z}}_{\mathrm{c}})$. In the frame of LES solvers, cut-off and Gaussian filters are other common alternative filtering options, but they are not detailed in this text for the sake of brevity.

2.2. Filtered Conservation Equations and Modeling

The presented filters can be applied to the transient conservation equations of a turbulent reacting flow to obtain a filtered LES version, which is formally very similar to their RANS counterparts [28]:

Conservation of mass:

$$\frac{\partial \overline{\mathsf{\rho }}}{\partial \mathrm{t}}+\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{u}}}_{\mathrm{i}}\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}=0$$

(4)

Conservation of momentum:

$$\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{u}}}_{\mathrm{i}}\right)}{\partial \mathrm{t}}+\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}+\frac{\partial \overline{\mathrm{p}}}{\partial {\mathrm{x}}_{\mathrm{j}}}=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\overline{\mathsf{\tau }}}_{\mathrm{i}\mathrm{j}}-\overline{\mathsf{\rho }}\left(\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}\right)\right)$$

(5)

Conservation of chemical species:

$$\frac{\partial \overline{\mathsf{\rho }}\tilde{Y} }{\partial \mathrm{t}}+\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{Y}}}_{\mathrm{k}}\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left(\overline{{\mathrm{V}}_{\mathrm{k},\mathrm{i}}{\mathrm{Y}}_{\mathrm{k}}}-\overline{\mathsf{\rho }}\left(\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{Y}}_{\mathrm{k}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{Y}}}_{\mathrm{k}}\right)\right)+{\overline{\dot{\mathsf{\omega }}}}_{\mathrm{k}},$$

(6)

where ${\mathrm{Y}}_{\mathrm{k}}$ and ${\dot{\omega }}_k$ denote the mass fraction and the chemical production rate of the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ chemical species. In addition, ${\mathrm{V}}_{\mathrm{k},\mathrm{i}}$ represents the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ component of the diffusion velocity ${\mathrm{V}}_{\mathrm{k}}$ for the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ chemical species.

Conservation of energy:

$$\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{h}}}_{\mathrm{s}}\right)}{\partial \mathrm{t}}+\frac{\partial \left(\overline{\mathsf{\rho }}{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{h}}}_{\mathrm{s}}\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}=\frac{\partial \overline{p}}{\partial \mathrm{t}}+\overline{{\mathrm{u}}_{\mathrm{i}}\frac{\partial \mathrm{p}}{\partial {\mathrm{x}}_{\mathrm{i}}}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left(\overline{\mathsf{\lambda }\frac{\partial \mathrm{T}}{\partial {\mathrm{x}}_{\mathrm{i}}}}-\overline{\mathsf{\rho }}\left(\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{s}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{h}}}_{\mathrm{s}}\right)\right)+\overline{{\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}\frac{\partial {\mathrm{u}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}}-\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left(\overline{\mathsf{\rho }\sum _{\mathrm{k}=1}^{\mathrm{N}}{\mathrm{V}}_{\mathrm{k},\mathrm{i}}{\mathrm{Y}}_{\mathrm{k}}{\mathrm{h}}_{\mathrm{s},\mathrm{k}}}\right)+{\overline{\dot{\mathsf{\omega }}}}_{\mathrm{T}}$$

(7)

In this equation, $\mathsf{\lambda }$ denotes the thermal conductivity, and ${\dot{\mathsf{\omega }}}_{\mathrm{T}}$ the enthalpy change rate due to chemical reactions. Several terms in the presented equations are unresolved and require modeling. The present work focuses on the subgrid turbulent fluxes, which have the general form $\widetilde{{\mathrm{u}}_{\mathrm{i}}\mathrm{q}}$. The most relevant ones are the subgrid stresses, which constitute the main unresolved term in the conservation of momentum equation:

$$\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}$$

(8)

This element is the only term requiring modeling in constant-density non-reactive flows. Therefore, it plays a key role in any application, and it has been the subject of intense research during the last decades. Smagorinsky [48] pioneered a closure strategy based on the Boussinesq approximation:

$$\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}=-2{\mathsf{\nu }}_{\mathrm{T}}\left({\tilde{\mathrm{S}}}_{\mathrm{i}\mathrm{j}}-\frac{{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}}{3}{\tilde{\mathrm{S}}}_{\mathrm{k}\mathrm{k}}\right)=-2{\mathrm{C}}_{\mathrm{s}}^2{\mathsf{\Delta }}^2‖{\tilde{\mathrm{S}}}_{\mathrm{i}\mathrm{j}}‖\left({\tilde{\mathrm{S}}}_{\mathrm{i}\mathrm{j}}-\frac{{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}}{3}{\tilde{\mathrm{S}}}_{\mathrm{k}\mathrm{k}}\right),$$

(9)

where ${\tilde{\mathrm{S}}}_{\mathrm{i}\mathrm{j}}=\frac{1}{2}\left(\frac{\partial {\tilde{\mathrm{u}}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{i}}}+\frac{\partial {\tilde{\mathrm{u}}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)$ denotes the filtered strain rate tensor, and ${\mathrm{C}}_{\mathrm{s}}$ is a model constant, originally assumed to be universal. The application of this model requires the a priori knowledge of ${\mathrm{C}}_{\mathrm{s}}$, which can be extracted from experimental data or DNS. This motivated the research of this value under different setups [49,50,51]. These studies revealed a good overall model performance, but the optimized constant ${\mathrm{C}}_{\mathrm{s}}$ varied by a factor of 2 between different experiments. In addition, this constant has been observed to be sensitive to the used filter width $\mathsf{\Delta }$ [52,53]. To circumvent these challenges, Germano and Lilly proposed a dynamic approach in which ${\mathrm{C}}_{\mathrm{s}}$ was no longer a constant but a dynamically adaptable coefficient. The actual value is then determined from the problem resolution at a test level with larger filter widths ${\mathsf{\Delta }}_{test}>\mathsf{\Delta }$. This strategy has proved excellent performance in various applications [54]. In compressible turbulence, modeling the diagonal part of the subgrid stresses i.e. $\widetilde{{\mathrm{u}}_{\mathrm{k}}{\mathrm{u}}_{\mathrm{k}}}-\widetilde{{\mathrm{u}}_{\mathrm{k}}}\widetilde{{\mathrm{u}}_{\mathrm{k}}}$ is required as well. For low Mach number flows, the diagonal component is usually absorbed into the filtered pressure field $\overline{\mathrm{p}} \left[28\right].$ Yoshizawa [55] proposed a model for the trace of the subgrid Reynolds stresses formally similar to Smagorinky's approach. However, the influence of this term is very limited, and its relevance is primarily restricted to highly compressible flows [56].

In reacting flows, the subgrid fluxes of chemical species (i.e. $\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{Y}}_{\mathrm{k}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{Y}}}_{\mathrm{k}}$) are unresolved as well. These scalars are usually modeled resorting to the gradient assumption [28], which requires determining the subgrid turbulent viscosity and Schmidt number. These values are commonly retrieved using strategies analogous to the Smagorinski formulation, and assuming constant turbulent Schmidt number. The subgrid heat flux (i.e. $\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{s}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{s}}$) is a relevant unresolved term within the energy conservation equation (7). A popular approach suggested by Eidson [57] models this term assuming proportionality between the energy transfer from resolved to subgrid levels and the filtered temperature gradient $\nabla \tilde{\mathrm{T}}$.

The introduced models were mainly conceived for incompressible non-reactive flows with low Damköhler numbers. Nonetheless, in high-pressure combustion applications, where the flame dynamics are very fast, strong interactions between turbulence and combustion may appear at a subgrid level. This circumstance can limit the validity of extended subgrid models due to the onset of unforeseen interactions. For example, the Boussinesq approximation prevents capturing reversed turbulent kinetic energy transfer, which is known to be relevant in turbulent combustion [18,58]. In addition, turbulent flames can display spectacular anisotropic features, even at the smallest vortical scales. Therefore, the combustion community has dedicated a significant effort in developing models able to cope with the challenges inherent to reacting flows.

Clark et. al [41] and Vreman et. al [42] derived gradient-based models which allow to express the subgrid stresses using the following similarity law:

$$\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}=\mathrm{C}{\mathsf{\Delta }}^2\left(\frac{\partial \widetilde{{\mathrm{u}}_{\mathrm{i}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\frac{\partial \widetilde{{\mathrm{u}}_{\mathrm{j}}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right)$$

(10)

This idea was further exploited by Pierce and Moin [43], who combined it with the notion of scale-similarity [59] to deduce the following general expression for subgrid covariances:

$$\widetilde{\mathsf{\varphi }\mathsf{\psi }}-\tilde{\mathsf{\varphi }}\tilde{\mathsf{\psi }}={\mathrm{C}}_{\mathsf{\varphi }\mathsf{\psi }}{\mathsf{\Delta }}^2\left(\frac{\partial \tilde{\mathsf{\varphi }}}{\partial {\mathrm{x}}_{\mathrm{k}}}\frac{\partial \tilde{\mathsf{\psi }}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right)$$

(11)

Models following this expression have been successfully used in turbulent combustion LES to determine relevant unresolved statistics, such as the subgrid variance of the turbulent mixture fraction $\widetilde{{\mathrm{Z}}^{″2}}$ [43,60,61] or the subgrid scalar dissipation rate $\tilde{\mathsf{\chi }}$ [61,62].

2.3. Turbulent Diffusion Flames

In turbulent diffusion flames, the propellants are injected separately and rely on turbulent and molecular diffusion to mix and burn. A comprehensive analysis of these physical systems poses a formidable challenge exceeding the scope of the present work. The reader is advised to review the available literature [63,64] for a detailed overview of turbulent non-premixed combustion. In this section, only the most relevant indicators to describe non-premixed flames are introduced due to their relevance in the model discussion.

2.3.1. Flame Index FI

Real combustion applications are neither purely premixed nor non-premixed. Instead, there is usually a coexistence between these two conditions. The flame index definition proposed by Yamashita et al. [65] aims to evaluate the local regime at which combustion takes place:

$$\mathrm{F}\mathrm{I}=\frac{\nabla {\mathrm{Y}}_{\mathrm{f}}\cdot \nabla {\mathrm{Y}}_{\mathrm{o}\mathrm{x}}}{\left|\nabla {\mathrm{Y}}_{\mathrm{f}}\right|\left|\nabla {\mathrm{Y}}_{\mathrm{o}\mathrm{x}}\right|},$$

(12)

where $\mathrm{Y}$ represents the mass fraction, and the subscripts $\mathrm{f}$ and $\mathrm{o}\mathrm{x}$ stand for fuel and oxidizer, respectively. If combustion occurs in premixed conditions, both reactants are depleted simultaneously, and their gradients are aligned, yielding $\mathrm{F}\mathrm{I}\approx 1$. The opposite happens if combustion is locally non-premixed, with the gradients of fuel and oxidizer being opposed. In such a case, the flame index tends to $\mathrm{F}\mathrm{I}\approx -1$. An example of the instantaneous flame index is provided in Figure 1. As can be seen, this parameter presents an almost binary distribution. This observation is in agreement with the trends found in one recent study [66]. The pattern near injection displayed in Figure 1 corresponds to the typical triple-flame structure [67], with a non-premixed branch surrounded by one lean premixed branch and one fuel-rich premixed branch.

2.3.2. Mixture Fraction $Z$

In a turbulent diffusion flame, the mixture composition continuously fluctuates in space and time due to the competition between the injection boundary conditions and turbulent diffusion processes. The mixture fraction $Z$ is a conserved scalar characterizing the mixing state. This parameter essentially represents the local proportion of fuel coming from the fuel stream with respect to the total flow. In the frame of hydrocarbon combustion, the mixture fraction of the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ chemical species, with ${\mathrm{n}}_{\mathrm{C}}$, ${\mathrm{n}}_{\mathrm{H}}$, and ${\mathrm{n}}_{\mathrm{O}}$ carbon, hydrogen, and oxygen atoms can be expressed as:

$${\mathrm{Z}}_{\mathrm{k}}=\frac{{\mathrm{A}}_{\mathrm{C}}{\mathrm{n}}_{\mathrm{C}}+{\mathrm{A}}_{\mathrm{H}}{\mathrm{n}}_{\mathrm{H}}}{{\mathrm{A}}_{\mathrm{C}}{\mathrm{n}}_{\mathrm{C}}+{\mathrm{A}}_{\mathrm{H}}{\mathrm{n}}_{\mathrm{H}}+{\mathrm{A}}_{\mathrm{O}}{\mathrm{n}}_{\mathrm{O}}}\approx \frac{12{\mathrm{n}}_{\mathrm{C}}+{\mathrm{n}}_{\mathrm{H}}}{12{\mathrm{n}}_{\mathrm{C}}+{\mathrm{n}}_{\mathrm{H}}+16{\mathrm{n}}_{\mathrm{O}}},$$

(13)

where ${\mathrm{A}}_{\mathrm{C}},{\mathrm{A}}_{\mathrm{H}}$ and ${\mathrm{A}}_{\mathrm{O}}$ denote the atomic weights of carbon, hydrogen, and oxygen, respectively. The individual mixture fractions can be combined to obtain the global value:

$$\mathrm{Z}={\sum }_{\mathrm{k}=1}^{\mathrm{n}}{\mathrm{Y}}_{\mathrm{k}}{\mathrm{Z}}_{\mathrm{k}},$$

(14)

where ${\mathrm{Y}}_{\mathrm{k}}$ is the mass fraction of the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ chemical species. Obviously, this parameter tends to unity towards the fuel and to zero close to the oxidizer injection. For methane–oxygen combustion, the stoichiometric oxidizer-to-fuel ratio is ${\left(\mathrm{O}/\mathrm{F}\right)}_{\mathrm{s}\mathrm{t}}\approx 4$, with a corresponding stoichiometric mixture fraction ${\mathrm{Z}}_{\mathrm{s}\mathrm{t}}\approx 0.2$. An example of the instantaneous mixture fraction can be observed in Figure 2. In this graphic, oxygen is injected at the center-left and methane at the top and bottom. The overall configuration is fuel-rich ($O/F=2$), this can be verified by observing the composition at downstream conditions, where the flow is well mixed, and the mixture fraction ranges from stoichiometric to mildly fuel-rich.

In the case of infinitely fast chemistry, combustion progress is strictly limited by the mixing configuration, which can be characterized by the mixture fraction and its dissipation rate $\mathsf{\chi }=\mathrm{D}\frac{\partial \mathrm{Z}}{\partial {\mathrm{x}}_{\mathrm{k}}}\frac{\partial \mathrm{Z}}{\partial {\mathrm{x}}_{\mathrm{k}}}$, representing the mixing intensity. If chemistry is very fast compared to turbulence, the flame dynamics can be decoupled from vortical interactions. In such cases, it is possible to solve the chemical problem independently and use tabulated results. This is the main assumption in the flamelet approach [68], which allows for representing all chemistry with the single scalar $\mathrm{Z}$. The capability of the mixture fraction to solely describe the mixing state illustrates the relevance of this parameter.

2.3.3. Progress Variable

The progress variable $\mathrm{c}$ represents the degree of combustion process completion. It ranges from zero to unity, corresponding to reactants, and products, respectively. In diffusion flames, it is often described using an expression of the following sort, as proposed by Bray et al. [69]:

$$\mathrm{C}=\frac{{\mathrm{Y}}_{\mathrm{c}}}{{\mathrm{Y}}_{\mathrm{c}}^{\mathrm{E}\mathrm{q}}\left(\mathrm{Z}\right)},$$

(15)

where ${\mathrm{Y}}_{\mathrm{c}}$ is the reaction progress variable [63,70], and ${\mathrm{Y}}_{\mathrm{c}}^{\mathrm{E}\mathrm{q}}\left(\mathrm{Z}\right)$ is the equilibrium value at a specific given mixture fraction. The definition of a suitable function for the reaction progress variable is strongly dependent on the used propellant combination and the operating conditions. For oxygen–methane combustion at high pressures, the following expression has been used in previous research [71,72]:

$${\mathrm{Y}}_{\mathrm{c}}=\frac{{\mathrm{Y}}_{\mathrm{C}{\mathrm{O}}_2}}{{\mathrm{M}}_{\mathrm{C}{\mathrm{O}}_2}}+\frac{{\mathrm{Y}}_{\mathrm{C}\mathrm{O}}}{{\mathrm{M}}_{\mathrm{C}\mathrm{O}}}+\frac{{\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}}{{\mathrm{M}}_{{\mathrm{H}}_2\mathrm{O}}}+\frac{{\mathrm{Y}}_{\mathrm{O}}}{{\mathrm{M}}_{\mathrm{O}}}+\frac{{\mathrm{Y}}_{\mathrm{H}}}{{\mathrm{M}}_{\mathrm{H}}}+\frac{{\mathrm{Y}}_{\mathrm{O}\mathrm{H}}}{{\mathrm{M}}_{\mathrm{O}\mathrm{H}}},$$

(16)

where ${\mathrm{M}}_{\mathrm{k}}$ denotes the molecular weight of the ${\mathrm{k}}^{\mathrm{t}\mathrm{h}}$ species. This definition is based on the formulation proposed by Pierce and Moin [73], which was applied in recent works [66] with similar setups as the current paper. Several variations have been incorporated to derive the expression in (16). The motivation behind these modifications was to account for specific chemical phenomena that take place in high-pressure combustion for space propulsion applications. An example of the instantaneous progress variable is displayed in Figure 3. This instantaneous 2D cross-section cut corresponds to the same instant as those in Figure 1 and Figure 2. Near injection, the progress variable reaches unity at the stoichiometric line, where combustion occurs in a non-premixed regime. At downstream positions, the flow is more homogeneous, and combustion primarily develops in fuel-rich premixed conditions. These figures illustrate the complexity and spatial variability of a real turbulent diffusion flame.

3. Simulation Setup

Direct numerical simulations of a methane-oxygen turbulent non-premixed flame at 20 bar were conducted to generate a database for subgrid turbulence dynamics investigation. The present section describes the simulation framework and the strategy followed to filter the resolved fields. The results of the described simulations are used to investigate the subgrid scale dynamics in the following sections.

3.1. Strategy

Two simulations are necessary to obtain the aimed turbulent diffusion flame. First, a precursor simulation is conducted to obtain realistic turbulent inlet boundary conditions. These fields are patched as the inlet boundary condition for the main simulation, where turbulent mixing and combustion take place. This second simulation is the one relevant within the scope of the present work. The overall scheme can be visualized in Figure 4. The described database has been used in previous research [74], where the reader can find a more detailed explanation concerning the simulation setup. This text is limited to a brief description for the sake of conciseness.

In the precursor simulation, periodic Synthetic Turbulence Generation (STG) takes place at the inlet, with the methodology discussed in [75]. This approach is a modification of the periodic STG proposed by Morschbach et al. [76], which is based on the work of Shur et al. [77]. These methods derive from the original formulation developed by Kraichnan [78], based on the superposition of several harmonics derived from a reference turbulent kinetic energy spectrum to generate synthetic turbulence. The STG development into realistic turbulence can be validated by examining classical statistics of turbulent flows such as higher-order derivatives and structure functions as proved in a previous work [75]. The turbulent characteristics were selected considering the results of previous RANS simulations of a scale methane rocket combustor [15]. It should be mentioned that the combustion pressure and the large eddy scales are roughly an order of magnitude smaller than what would be expected in real applications. These deviations are primarily motivated by the computational cost. Nonetheless, the chosen conditions are closer to realistic rocket combustion than what is usually found in the scientific literature. A summary of the turbulent characteristics at injection is provided in Table 1.

The main simulation aims at recreating mixing and combustion nearby injection and far away from the combustion chamber walls. A stationary diffusion flame is simulated in a cuboid domain with dimensions 0.2 x 0.2 x 1.2 mm and resolved with 192 x 192 x 1152 cubic cells with constant geometry. This configuration corresponds to a grid spacing ${\mathsf{\Delta }}_0=1.04167\mathsf{\mu }\mathrm{m}$. This cell size ensures enough resolution for the smallest Kolmogorov eddies [44]. The used resolution was proved to ensure convergence of the turbulent statistics in a mesh sensitivity analysis virtually identical setups [72,79]. This grid size provides a resolution of the stoichiometric premixed flame front with roughly 5.3 cells. Nevertheless, in turbulent methane-oxygen diffusion flames the thinnest realized flamelets at least double this value, as evaluated in previous works [72,80]. Therefore, the effective flame front resolution is in the order of 10-20 cells. A resolution with a finer grid size ($0.52083\mathsf{\mu }\mathrm{m}$ and 384 x 384 x 2304 cells) was conducted to perform a mesh sensitivity analysis on the results, with marginal differences found as presented in section 4.2.1.

The transient reactive turbulent flow In the main simulation is computed using the reactive solver EBI-DNS [81,82,83,84]. This code package is embedded in the open-source software OpenFOAM [85,86] and was developed by the Engler-Bunte Institute (EBI) for combustion technology. EBI-DNS solves the unsteady compressible conservation equations using the Finite Volume Method (FVM) [87,88]. The capabilities of this software have been established in a large number of research works conducted during the last years [89,90,91,92,93]. In a recent work [94], Zirwes et al. validated the EBI-DNS code comparing its accuracy with state-of-the art explicit reactive solvers in a benchmark case of hydrogen flames interacting with Taylor-Green vortices. Cantera [95] is used to compute detailed chemistry and transport properties. This code implements the mixture-averaged transport model as described by Kee et al. [96]. The diffusive mass flux is determined with the mixture-average model derived from the Maxwell-Stefan equations [96]. Pressure and temperature diffusion (i.e. Soret effect) are neglected. The methane oxidation process is calculated with the finite rate method using a complex chemical mechanism developed by Slavinskaya et al. [97]. This reaction scheme was conceived for space propulsion applications at high pressures, and it has been extensively used in several numerical research in this field [10,98]. The relevant combustion parameters are displayed in Table 2. Numerical diffusion is a major concern in unsteady numerical simulations, especially if implicit schemes are used. Zirwes et. al [93] demonstrated that the used solver exhibits negligible numerical diffusion with a second order spatial discretization scheme and low CFL numbers. A mesh sensitivity analysis in an almost identical setup [72] further supports the low risk of numerical diffusion polluting the solve’'s fidelity. In this previous work, marginal changes in averaged first and second order statistics were found for finer and coarser meshes. Since numerical diffusion is strongly dependent on the grid size, it is highly improbable that it plays a relevant role in the present context. The virtually identical results obtained with the finer grid size further support the low risk of numerical diffusion.

3.2. Filtering

The fully resolved results can be filtered using different filter widths $\mathsf{\Delta }$ to obtain validation data subgrid models. This enables the measurement of the actual subgrid turbulent quantities. An output example for the filtering process can be visualized in Figure 5, where the filtered temperature field with different widths $\mathsf{\Delta }$ can be compared against the entirely resolved direct numerical simulation. A summary of the used filters is provided in Table 2. As can be seen, large filter widths over DNS fields produces an excessive pixelation. Such results are physically unrealistic and cannot be the outcome of a numerically sound LES. However, if a top-hat filter with an excessive width is applied over a well resolved field, such unnatural results are feasible. The inaccuracy of such filter widths will be further commented on in the results section.

Figure 5. Examples of different filtering widths: Direct numerical simulation $\mathsf{\Delta }={\mathsf{\Delta }}_0$ (a), Filter width $\mathsf{\Delta }=4{\mathsf{\Delta }}_0$ (b), Filter width $\mathsf{\Delta }=12{\mathsf{\Delta }}_0$ (c), Filter width $\mathsf{\Delta }=24{\mathsf{\Delta }}_0$ (d).

Figure 5. Examples of different filtering widths: Direct numerical simulation $\mathsf{\Delta }={\mathsf{\Delta }}_0$ (a), Filter width $\mathsf{\Delta }=4{\mathsf{\Delta }}_0$ (b), Filter width $\mathsf{\Delta }=12{\mathsf{\Delta }}_0$ (c), Filter width $\mathsf{\Delta }=24{\mathsf{\Delta }}_0$ (d).

Table 3. Summary of the used filters.

Simulation	Number of cells [x y z]	$\mathbf{Grid} \mathbf{size} [\mathsf{\Delta }/{\mathsf{\delta }}_{\mathbf{L}0}$]	$\mathbf{Grid} \mathbf{size} [\mathsf{\Delta }/{\mathsf{\Delta }}_{\mathbf{D}\mathbf{N}\mathbf{S}}$]
$\mathrm{D}\mathrm{N}\mathrm{S}$	192 x 192 x 1152	0.189	1
$\mathrm{F}1$	96 x 96 x 756	0.378	2
$\mathrm{F}2$	64 x 64 x 384	0.568	3
$\mathrm{F}3$	48 x 48 x 288	0.759	4
$\mathrm{F}4$	32 x 32 x 192	1.135	6
$\mathrm{F}5$	24 x 24 x 144	1.514	8
$\mathrm{F}6$	16 x 16 x 96	2.271	12
$\mathrm{F}7$	12 x 12 x 72	3.028	16
$\mathrm{F}8$	8 x 8 x 48	4.514	24
$\mathrm{F}9$	4 x 4 x 24	9.083	48

Table 3. Summary of the used filters.

Simulation	Number of cells [x y z]	$\mathbf{Grid} \mathbf{size} [\mathsf{\Delta }/{\mathsf{\delta }}_{\mathbf{L}0}$]	$\mathbf{Grid} \mathbf{size} [\mathsf{\Delta }/{\mathsf{\Delta }}_{\mathbf{D}\mathbf{N}\mathbf{S}}$]
$\mathrm{D}\mathrm{N}\mathrm{S}$	192 x 192 x 1152	0.189	1
$\mathrm{F}1$	96 x 96 x 756	0.378	2
$\mathrm{F}2$	64 x 64 x 384	0.568	3
$\mathrm{F}3$	48 x 48 x 288	0.759	4
$\mathrm{F}4$	32 x 32 x 192	1.135	6
$\mathrm{F}5$	24 x 24 x 144	1.514	8
$\mathrm{F}6$	16 x 16 x 96	2.271	12
$\mathrm{F}7$	12 x 12 x 72	3.028	16
$\mathrm{F}8$	8 x 8 x 48	4.514	24
$\mathrm{F}9$	4 x 4 x 24	9.083	48

The filtering of physical fields over the fully resolved DNS domain provides new fields with a different number of cells, as illustrated in Figure 5. Each of these coarser cells provides one data point to investigate the subgrid turbulent dynamics. The number of obtained data points ranges from 384 using $\mathrm{F}9$ to 5308416 using $\mathrm{F}1$. A total of 40 timesteps, with a temporal spacing of $\mathsf{\mu }\mathrm{s}$ were filtered to obtain sufficient observations of subgrid turbulent quantities. This implies 15360 data points for the worst case ($\mathrm{F}9$), sufficient to perform a meaningful statistical analysis. Due to its simplicity, top-hat filters have been used. The filtering operations were performed using gaussian quadratic integration.

4. Gradient Model Enhancement

The present section introduces an extension of the gradient model for unresolved subgrid turbulent fluxes aiming at high-pressure combustion applications. The application scope corresponds to insufficiently resolved flame fronts, in which both turbulent and chemical processes interact at a subgrid scale. The discussion is structured into three parts. First, the model is presented along with its physical implications. In the second sub-section, the model's performance is studied using the available DNS as a validation resource. Finally, the dynamic behavior of the model coefficient is analyzed and scaling laws are proposed.

4.1. Model Description

In LES of variable-density flows, unresolved fluxes present the following general form:

$$\widetilde{{\mathrm{u}}_{\mathrm{i}}\mathrm{q}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}\tilde{\mathrm{q}}=\widetilde{{\mathrm{u}}_{\mathrm{i}}^{″}{\mathrm{q}}^{″}},$$

(17)

where $\mathrm{q}$ represents the advected quantity. Here $\mathrm{q}″$ denotes the instantaneous subgrid fluctuations over the filtered reference value $\tilde{\mathrm{q}}$ at a given position i.e. ${\mathrm{q}}^{″}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)=\mathrm{q}\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)-\tilde{\mathrm{q}}(\mathrm{x},\mathrm{y},\mathrm{z})$. The application of the gradient model yields the following closure equation:

$$\widetilde{{\mathrm{U}}_{\mathrm{i}}\mathrm{q}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}\tilde{\mathrm{q}}=\widetilde{{\mathrm{u}}_{\mathrm{i}}^{″}{\mathrm{q}}^{″}}={\mathrm{C}}^2{\mathsf{\Delta }}^2\nabla {\tilde{\mathrm{u}}}_{\mathrm{i}}\cdot \nabla \tilde{\mathrm{q}},$$

(18)

where $\mathrm{C}$ is a model constant that shall be determined using DNS or experimental data. This coefficient implicitly embeds the correlation between the projection of the involved quantities $\nabla {\tilde{\mathrm{u}}}_{\mathrm{i}}\cdot \nabla \tilde{\mathrm{q}}$ and the modeled subgrid covariance. A high correlation takes place if the subgrid dynamics are similar to the ones in larger scales. In diffusion flames combustion may take place in non-premixed or premixed conditions depending on the local conditions. This is particularly relevant in space propulsion applications, where the high reactivity of the propellants increases the viability of combustion in locally premixed conditions [99]. In LES, the actual local conditions may be evaluated using the filtered flame index:

$$\widetilde{\mathrm{F}\mathrm{I}}=\frac{\nabla {\tilde{\mathrm{Y}}}_{\mathrm{f}}\cdot \nabla {\tilde{\mathrm{Y}}}_{\mathrm{o}\mathrm{x}}}{\left|\nabla {\tilde{\mathrm{Y}}}_{\mathrm{f}}\right|\left|\nabla {\tilde{\mathrm{Y}}}_{\mathrm{o}\mathrm{x}}\right|}$$

(19)

Since the combustion dynamics are significantly different depending on the considered regime, this information can be incorporated to enhance the models' accuracy. Hence, it is possible to extend the previous expression to define the following subgrid scale model:

$$\widetilde{{\mathrm{u}}_{\mathrm{i}}\mathrm{q}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}\tilde{\mathrm{q}}=\widetilde{{\mathrm{u}}_{\mathrm{i}}^{″}{\mathrm{q}}^{″}}=\left\{\begin{array}{c}{\mathrm{C}}_{\mathrm{m}}^2{\mathsf{\Delta }}^2\nabla {\tilde{\mathrm{u}}}_{\mathrm{i}}\cdot \nabla \tilde{\mathrm{q}}, \mathrm{if} \widetilde{\mathrm{F}\mathrm{I}}>0\\ {\mathrm{C}}_{\mathrm{d}}^2{\mathsf{\Delta }}^2\nabla {\tilde{\mathrm{u}}}_{\mathrm{i}}\cdot \nabla \tilde{\mathrm{q}}, \mathrm{if} \widetilde{\mathrm{F}\mathrm{I}}<0\end{array}\right.,$$

(20)

where ${\mathrm{C}}_{\mathrm{m}}$ and ${\mathrm{C}}_{\mathrm{d}}$ are coefficients requiring determination. The physical motivation for the derived model is grounded on the equilibrium assumption and dominance of global flame dynamics at the filtered level. Therefore, it is expected that the coefficients $\mathrm{C}$ are not constant but depend on the overlapping between the grid size and the chemical scales. This dependency can be modeled with a general power law of the sort:

$$\mathrm{C}\left(\mathsf{\Delta }\right)=\mathsf{\alpha }{\left(\frac{{\mathsf{\delta }}_{\mathrm{L}0}}{\mathsf{\Delta }}\right)}^{\mathsf{\beta }},$$

(21)

where the exponent $\beta$ is assumed to be universal. The coefficient $\alpha$ can be determined dynamically using a coarser test scale, following a similar method as suggested by Germano [100] and other researchers. The function $\mathrm{C}\left(\mathsf{\Delta }\right)$ is detailed in the following text section.

4.2. Model Validation

The presented model was calibrated using the available results from the fully resolved direct numerical simulations. The obtained coefficients for premixed and diffusion-driven combustion i.e. ${\mathrm{C}}_{\mathrm{m}}$, and ${\mathrm{C}}_{\mathrm{d}}$ are presented along with an evaluation of the model performance. Finally, the dynamic determination of the coefficients $\mathrm{C}\left(\mathsf{\Delta }\right)$ is discussed along with its physical implications.

4.2.1. Reynolds Stresses $\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{u}}}_{\mathrm{j}}$ and $\widetilde{{\mathrm{u}}_{\mathrm{k}}{\mathrm{u}}_{\mathrm{k}}}-{\tilde{\mathrm{u}}}_{\mathrm{k}}{\tilde{\mathrm{u}}}_{\mathrm{k}}$

The subgrid Reynolds stresses are probably the most important unresolved quantities in LES. This modeling is detailed here to illustrate the parameter adjustment process and the evaluation of the model performance.

For each recorded timestep, subgrid and filtered quantities are available. The role of an LES model is to predict the former with the latter, using expressions of the sort of (18). For every cell, it is possible to determine the filtered flame index along with the observed subgrid quantities and the elements proposed in this model. Storing this information, it is possible to generate a database like the one presented in Figure 8. Here, data from a single timestep is displayed to prevent the loss of perspective. However, the actual statistical analyses were performed using data from several timesteps to minimize the statistical error. In the graph of Figure 8, the vertical axis corresponds to the quantity that has to be modeled (i.e. the subgrid Reynolds stress tensor), whereas the horizontal one contains the model-dependent parameters. The color of each measurement is related to the filtered flame index $\widetilde{\mathrm{F}\mathrm{I}}$. As can be seen, there is a clear correlation between both variables. For each case, the slope of the scattered data cloud corresponds to the model coefficient $\mathrm{C}$, which is used to denote ${\mathrm{C}}_{\mathrm{m}}$ or ${\mathrm{C}}_{\mathrm{d}}$ indistinctively. The cells in a premixed regime ( $\widetilde{\mathrm{F}\mathrm{I}}>0$) present a slightly smaller slope compared with the ones that burn in non-premixed conditions ( $\widetilde{\mathrm{F}\mathrm{I}}<0$). This is the motivation for using different coefficients, as mentioned in the previous section.

Figure 6. Example of instantaneous results for modeling the subgrid stress ${\mathsf{\tau }}_{\mathrm{x}\mathrm{y}}$ $\mathsf{\Delta }=6{\mathsf{\Delta }}_0$.

Figure 6. Example of instantaneous results for modeling the subgrid stress ${\mathsf{\tau }}_{\mathrm{x}\mathrm{y}}$ $\mathsf{\Delta }=6{\mathsf{\Delta }}_0$.

Besides the different sensitivities, it is evident that dispersion is greater for the cells in the non-premixed regime. This trend is consistently observed for most of the modeled subgrid fluxes, indicating a worse model performance in the regions where combustion is driven by diffusion. This behavior is probably caused by the fact that microscale mixing plays a significant role in the subgrid dynamics when the flow is locally non-premixed. Since these processes are partially uncoupled with the filtered (i.e. larger scale) quantities, they incorporate additional noise, worsening the model performance. For the different used grid sizes, the model constants at both premixed and non-premixed conditions were fitted using the least squares method. The obtained coefficients and model performance are presented in Figure 7.

The coefficient of determination ${\mathrm{R}}^2$ is a valuable indicator of the model quality. However, this parameter is very sensitive to the outlier observations. Hence, it is more convenient to visualize the complete predicted and observed fields as displayed in Figure 8, to have a more holistic perspective concerning the model performance. In this graph, the measured and modeled subgrid stresses ${\mathsf{\tau }}_{\mathrm{x}\mathrm{y},\mathrm{s}\mathrm{g}\mathrm{s}}$ are presented with different filter widths. For a fine resolution $\mathsf{\Delta }=3{\mathsf{\Delta }}_0$, the model perfectly predicts the subgrid stress, as one could expect from ${\mathrm{R}}^2>0.9$. For the filter width $\mathsf{\Delta }=16{\mathsf{\Delta }}_0,$ the model has an average coefficient of determination in the order ${\mathrm{R}}^2\approx 0.55$. Although noticeable disparities appear at certain locations, the model accurately captures the subgrid quantities both quantitatively and qualitatively. In general, for ${\mathrm{R}}^2>0.7$, the model can be rendered as excellent, and for $0.7>{\mathrm{R}}^2>0.4$, the performance may be classified as good.

Figure 8. Observed and modeled subgrid Reynolds stress ${\tau }_{xy,sgs}$: Observed value with filter width $\mathsf{\Delta }={3\mathsf{\Delta }}_0$ (a), predicted value with filter width $\mathsf{\Delta }=3{\mathsf{\Delta }}_0$ (b), observed value with filter width $\mathsf{\Delta }={16\mathsf{\Delta }}_0$ (c), Predicted value with filter width $\mathsf{\Delta }=12{\mathsf{\Delta }}_0$ (d).

Figure 8. Observed and modeled subgrid Reynolds stress ${\tau }_{xy,sgs}$: Observed value with filter width $\mathsf{\Delta }={3\mathsf{\Delta }}_0$ (a), predicted value with filter width $\mathsf{\Delta }=3{\mathsf{\Delta }}_0$ (b), observed value with filter width $\mathsf{\Delta }={16\mathsf{\Delta }}_0$ (c), Predicted value with filter width $\mathsf{\Delta }=12{\mathsf{\Delta }}_0$ (d).

In the inlet regions, there is a significant discrepancy between the model prediction and the observed subgrid stresses. This is caused by the fact that the flow is non-reactive in these locations. This violates the central premise of the developed model. For this reason, points with a filtered progress variable $\tilde{c}<0.02$ have been excluded from the statistical analysis and model performance evaluations.

A finer grid has been used to perform a mesh sensitivity analysis over the entire procedure. The grid size of the refined mesh is half the one of the reference DNS. The obtained value for the coefficient $C_m$ in both the reference simulation and the refined one can be compared in Figure 9. Only the coefficients for the smaller filter sizes are displayed since the results for larger filters are virtually indistinguishable. As can be seen, there is a minor error in the coefficient for the filter with $\mathsf{\Delta }=0.378{\delta }_{L0}$. More specifically, the finer grid provides a value a 6.1 % larger compared to the reference case, for the smallest filter width. This increase originates from the small amount of unresolved turbulent processes in the reference case. The discrepancy between the reference case and the one with a refined grid tends asymptotically to zero as larger filter widths are considered, illustrating that the amount of subgrid information contained in both simulations is virtually identical.

Due to the large computational cost and manipulation complexity of the finer database, the mesh sensitivity analysis was only performed over the subgrid Reynolds stress tensor. Nonetheless, these results prove that the reference grid size is small enough to ensure accurate results.

4.2.2. Enthalpy Flux $\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{h}}_{\mathrm{s}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{h}}}_{\mathrm{s}}$

The subgrid enthalpy flux is one of the many terms that require modeling in the energy conservation equation. The devised model was fitted following the previously described procedure. The model constants and behavior present similar values to those obtained for the Reynolds stresses for small filter widths. The main difference between these results and the ones reported for the subgrid Reynolds stresses is that the model constants and performance are more alike for premixed and non-premixed conditions. This similarity is probably caused by the greater sensitivity of Reynolds stresses to composition changes. Microscale mixing processes have an immediate effect on the subgrid Reynolds stresses due to the viscosity differences. The conversion of mixing variability into enthalpy changes requires chemical reactions, which are limited by the local mixture fraction. This additionally required step can dampen the disruptive effect of microscale mixing into the model performance.

Figure 10. Coefficients and model performance for the subgrid enthalpy flux.

Figure 10. Coefficients and model performance for the subgrid enthalpy flux.

4.2.3. Species Fluxes $\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{Y}}_{\mathrm{k}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}{\tilde{\mathrm{Y}}}_{\mathrm{k}}$

The investigation of the error performance for unresolved species fluxes is significantly challenging. A comprehensive approach would require considering every chemical species in the reaction mechanism. To prevent the loss of sight, only a few selected species are analyzed. The analysis is divided into major species and intermediate i.e., minor species. The summary of the performed adjustment is displayed in Figure 11. The most relevant difference with the previously reported results is the model performance at small filter widths. This result might appear counterintuitive at first impression. If small grid sizes are considered, the behavior of some quantities is almost laminar, and some unresolved terms tend to be zero. Hence, the fitting problem $y=ax$ has to be solved in a scenario where $y$ and $x$ converge to zero. The low model performance is caused by the ill-posed nature of this mathematical problem in such a context. To verify this claim, the negligibility of the unresolved fluxes has to be demonstrated. With this goal in mind, the following indicator of the subgrid scale fluxes intensity is defined:

$${\mathsf{\varphi }}_{\mathrm{s}\mathrm{g}\mathrm{s}}\left(\mathrm{q}\right)=\left|\frac{\widetilde{{\mathrm{u}}_{\mathrm{i}}{\mathrm{q}}_{\mathrm{k}}}-{\tilde{\mathrm{u}}}_{\mathrm{i}}\tilde{\mathrm{q}}}{{\tilde{\mathrm{u}}}_{\mathrm{i}}\tilde{\mathrm{q}}}\right|,$$

(22)

which tends to zero when no subgrid modeling for the flux of the quantity $q$ is required. This parameter was evaluated for the subgrid water flux at different filter widths. For $\mathsf{\Delta }=2{\mathsf{\Delta }}_0$ $,{\mathsf{\varphi }}_{\mathrm{s}\mathrm{g}\mathrm{s}}\left({\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}\right)$ ranges from ${10}^{-5}$ to ${10}^{-2}$ through all the domain, with a median value ${\mathsf{\varphi }}_{\mathrm{s}\mathrm{g}\mathrm{s}}\left({\mathrm{Y}}_{{\mathrm{H}}_{2\mathrm{O}}}\right)=1.3\cdot {10}^{-4}$. These results indicate the negligibility of scalar flux at the specified filter width. Neglecting any modeling would have a minimal impact on the species conservation transport equation. For $\mathsf{\Delta }=6{\mathsf{\Delta }}_0$, the subgrid scale intensity becomes relevant, with ${\mathsf{\varphi }}_{\mathrm{s}\mathrm{g}\mathrm{s}}$ presenting a median value ${\mathsf{\varphi }}_{\mathrm{s}\mathrm{g}\mathrm{s}}\left({\mathrm{Y}}_{{\mathrm{H}}_2}\mathrm{O}\right)\approx 0.07$, which justifies the need for modeling. With this filter width, the model works already very well with $R^2\in \left[\mathrm{0.7,0.8}\right]$ depending on the considered species.

In the analysis summarized in Figure 11, there appears to be no significant difference between the major and minor species. However, for premixed combustion conditions, the growth of the model coefficient ${\mathrm{C}}_{\mathrm{m}}$ for species $\mathrm{O},\mathrm{H},$ and $\mathrm{O}\mathrm{H}$ appears to stop after $\mathsf{\Delta }/{\mathsf{\delta }}_{\mathrm{L}0}>3$. The profile of these species strongly depends on the local mixture fraction, and it has a non-monotonic profile through the flame front. These circumstances may decrease the correlation between the resolved processes and the subgrid dynamics. One possible solution would be to consider different coefficients for lean and fuel-rich conditions, enabling adaptation to the varying chemistry behavior in these contexts.

4.3. Dynamic Model Formulation

The previous analyses revealed an important sensitivity of the model coefficients $C_{m/d}$ to the grid size $\mathsf{\Delta }$. This dependency can be addressed with a dynamic approach, in which the problem is resolved and filtered with a coarser filter width to derive the model coefficients. However, scaling laws are needed to determine the coefficients at the targeted resolution. The scaling laws for the model coefficients is studied in the present section using the available results. The behavior of coefficient $\mathrm{C}$ for premixed combustion (i.e. ${\mathrm{C}}_{\mathrm{m}}$) as a function of the grid size $\mathsf{\Delta }$ is provided in Figure 12 for different subgrid fluxes. Due to its similarity, the function of the coefficient ${\mathrm{C}}_{\mathrm{d}}$ is not reproduced to prevent overclutter. As can be seen, the coefficient follows identical trends for all the unresolved fluxes. This result hints the existence of a universal power law dependency with a value $\mathsf{\beta }\approx 0.36$. If this parameter is known a priori, results at test levels can be used to determine $\mathsf{\alpha }$ in (21), and consequently obtain the functions ${\mathrm{C}}_{\mathrm{d}}\left(\mathsf{\Delta }\right)$ and ${\mathrm{C}}_{\mathrm{m}}\left(\mathsf{\Delta }\right)$.

For small filter widths, a noticeable deviation with respect to the main trend can be observed. More specifically, the function appears to become linear for $\mathsf{\Delta }<{\mathsf{\delta }}_{\mathrm{L}0}$. Hence, a linear modeling of the sort $\mathrm{C}={\mathsf{\alpha }}_{\mathsf{\lambda }}+{\mathsf{\beta }}_{\mathsf{\lambda }}({\mathsf{\delta }}_{\mathrm{L}0}/\mathsf{\Delta })$ appears to be better suited in these cases. The results of assuming this dependency for the constant modeling can be visualized in Figure 13. As it can be seen, linearization provides an almost perfect fitting for the region where the grid size is smaller than the laminar flame front thickness at stoichiometric conditions ${\mathsf{\delta }}_{\mathrm{L}0}$.

Based on these results, if the filter width is known to be smaller than the laminar flame thickness, a linear dependency should be assumed. One drawback of this approach is that it requires performing two different tests at already small filter widths, which can be computationally expensive. This problem can be circumvented by linking the behavior in the linear region ($\mathsf{\Delta }<{\mathsf{\delta }}_{\mathrm{L}0}$) with the universal dependency law. If ${\mathrm{C}}^1$ continuity is assumed for $\mathrm{C}\left(\mathsf{\Delta }\right),$ then the derivatives of both models should meet at the transition point ${\mathsf{\Delta }}_t$:

$$\frac{\partial \mathrm{C}}{\partial \mathsf{\Delta }}=\frac{-{\mathsf{\beta }}_{\mathsf{\lambda }}{\mathsf{\delta }}_{\mathrm{L}0}}{{\mathsf{\Delta }}_{\mathrm{t}}^2}=\frac{-\mathsf{\beta }\mathsf{\alpha }{\mathsf{\delta }}_{\mathrm{L}0}^{\beta }}{{\mathsf{\Delta }}_{\mathrm{t}}^{\mathsf{\beta }+1}},$$

(23)

where ${\mathsf{\Delta }}_t$ indicates the filter width at which the transition from potential to linear law occurs. The continuity condition can be used to obtain ${\mathsf{\beta }}_{\mathsf{\lambda }}$ as a function of this point$:$

$${\mathsf{\beta }}_{\mathsf{\lambda }}=\mathsf{\beta }\mathsf{\alpha }{\left(\frac{{\mathsf{\delta }}_{\mathrm{L}0}}{{\mathsf{\Delta }}_{\mathrm{t}}}\right)}^{\mathsf{\beta }-1}\approx \mathsf{\beta }\mathsf{\alpha }$$

(24)

In this approximation, it was assumed that the transition filter width corresponds to the laminar flame thickness at stoichiometric conditions. The validity of this premise can be tested with the available results. With the obtained values for $\mathsf{\alpha }$ in Figure 12, if we set ${\mathsf{\Delta }}_{\mathrm{t}}={\mathsf{\delta }}_{\mathrm{L}0}$, ${\mathsf{\beta }}_{\mathsf{\lambda }}\approx 0.14$ and ${\mathsf{\beta }}_{\mathsf{\lambda }}\approx 0.086$ are obtained for the steam mass fraction flux and the other unresolved fluxes. These values provide a rather small discrepancy when compared with the actual sensitivity coefficients displayed in Figure 13. The expression in (24) provides an almost perfect prediction for the coefficients ${\mathsf{\beta }}_{\mathsf{\lambda }}$ if the transition filter width is taken as ${\mathsf{\Delta }}_{\mathrm{t}}\approx {\mathsf{\delta }}_{\mathrm{L}0}/1.4$, which indicates that the change to a linear regime occurs for scales slightly narrower than the laminar flame thickness at stoichiometric conditions. This statement can be qualitatively verified attending to Figure 13.

The discussed variations in the function $\mathrm{C}\left(\mathsf{\Delta }\right)$ embed profound physical motivations. It shall be considered that the size of the realized flamelets ${\delta }_c$ is roughly 3 times larger than laminar flame thickness at stoichiometric conditions. Hence, the found threshold maybe rewritten as roughly ${\mathsf{\Delta }}_{\mathrm{t}}\approx {\mathsf{\delta }}_{\mathrm{c}}/4$. If The filter width is significantly smaller than the flame thickness, it is safe to assume that the filtering operator $F$ is entirely applied either within the flame or outside it. In such cases, there is consistency between the dynamics at a subgrid level and at large scales. For larger values, the possibility of applying the filter in regions which are only partially reactive becomes relevant. In this scenario, the existence of different subgrid dynamics of reactive and non-reactive zones degrades the similarity between subgrid fluxes and filtered quantities. This effect is displayed in the schematic of Figure 14. If the filter width is significantly smaller than the reaction thickness, then it is likely to be applied within a volume which falls entirely either within or outside the flame.

5. Validation with

The development and fitting of the presented model is based on the individual evaluation of instantaneous subgrid fluxes. However, it is relevant to investigate the integral performance of the proposed model in the frame of a complete numerical simulation. This section is devoted to the assessment of the model’s global validity using a reference DNS database.

5.1. Simulation Description

The case chosen for the global validation is a recessed methane-oxygen injector which has been studied in recent research works [80,101]. These simulations recreate the turbulent combustion in methane – oxygen injectors, using realistic geometrical features, and have been validated with experimental data in a comparable setup. The proposed models in this paper are limited to subgrid fluxes. Despite their remarkable relevance, there are multiple additional unresolved subgrid statistics defining the transient development of the conservation equations. The strong intertwining between the various unclosed terms complicates the global assessment of the enhancements in numerical accuracy provided by an improvement in the subgrid turbulent flux modelling. The following three simulations were used in the present study to investigate the sensitivity of the numerical accuracy to the usage of subgrid adaptive gradient models:

• Direct numerical simulation (DNS).
• Large eddy simulation with no modelling of subgrid turbulence ($\mathrm{L}\mathrm{E}{\mathrm{S}}_0$).
• Large eddy simulation with modelling of subgrid turbulence at downstream positions ($\mathrm{L}\mathrm{E}{\mathrm{S}}_{\mathrm{T}\mathrm{F}}$).

The simulation $\mathrm{L}\mathrm{E}{\mathrm{S}}_0$ uses no modelling at all for subgrid quantities, and it is essentially DNS with insufficient resolution. The third simulation i.e. $\mathrm{L}\mathrm{E}{\mathrm{S}}_{\mathrm{T}\mathrm{F}}$ models the subgrid fluxes with the models presented in the previous section at downstream positions. Both simulations $\mathrm{L}\mathrm{E}{\mathrm{S}}_0$, and $\mathrm{L}\mathrm{E}{\mathrm{S}}_{\mathrm{T}\mathrm{F}}$ use a grid size $\mathsf{\Delta }=8{\mathsf{\Delta }}_{\mathrm{D}\mathrm{N}\mathrm{S}}$.

The model coefficients take the following value:

$$\mathrm{C}={\mathrm{f}}_{\mathrm{c}}{\mathrm{f}}_{\mathrm{g}}{\mathrm{C}}_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}$$

(25)

where $f_c$, and $f_g$ are coefficients defining the local application of the subgrid models, and $C_{model}$ is the reference coefficient for the modelled subgrid flux. The coefficient $f_c$ controls the restriction of the model application to reactive regions using the following formulation:

$${\mathrm{f}}_c=\left\{\begin{array}{cc}3{\left(\frac{\tilde{c}}{0.05}\right)}^2-2{\left(\frac{\tilde{c}}{0.05}\right)}^3,\hfill & if \tilde{c}<0.05\hfill \\ 1,\hfill & if \tilde{c}\ge 0.05\hfill \end{array}\right.$$

(26)

Which enables a smooth transition from the non-reactive to the reactive case. The factor $f_g$ is devised in a similar fashion, defining a transition from $z=d_0$ to $z=4d_0$, as illustrated in Figure 15:

$${\mathrm{f}}_{\mathrm{g}}=\left\{\begin{array}{cc} 0,\hfill & if z\le d_0\hfill \\ 3{\left(\left(\frac{z}{d_0}-1\right)/3\right)}^2-2{\left(\left(\frac{z}{d_0}-1\right)/3\right)}^3,& if{d}_0<z<4d_0\hfill \\ 1,\hfill & if z\le d_0\hfill \end{array}\right.$$

(27)

The main motivation of this transition is to avoid the application of the presented model in regions close to walls, where dedicated treatment of wall turbulence is better approached with other alternatives.

5.2. Results Comparison

The average temperature over time and radial direction is used to compare the output of the different simulations. This scalar can be taken as a robust marker of the overall combustion development. The axial evolution of this quantity for the three considered simulations is displayed in Figure 16. As can be seen, the complete neglection of subgrid effects is translated into a substantial underprediction of the combustion progress ($LE{S}_0$). This is caused by the insufficiently small grid size which effectively reduces the turbulent Reynolds number, thereby suppressing microscale diffusion and mixing. The incorporation of the models for subgrid turbulent fluxes increases the axial temperature gradient, bringing the DNS and $\mathrm{L}\mathrm{E}{\mathrm{S}}_{\mathrm{T}\mathrm{F}}$ outputs closer. Nonetheless, a steady state offset between these two simulations remains. This discrepancy is attributed to absence of any modelling within $z\in [-d_0,d_0]$, and the unresolved subgrid quantities that are not turbulent fluxes, which remain non-modelled through all the domain of the simulation $LE{S}_{TF}$. These results demonstrate the capability of the presented gradient models to improve the accuracy of LES in reactive turbulent flows.

6. Conclusions and Outlook

The present research has addressed the application of the gradient model for subgrid turbulent fluxes in turbulent combustion conditions representative of space propulsion applications. An extension of the gradient model considering the local burning regime has been introduced. Hence, the determination of subgrid covariances depends on whether combustion takes place at premixed or non-premixed conditions. It is found that the gradient model performs significantly better if combustion is locally premixed, and that the magnitude of subgrid fluxes is slightly greater in zones where combustion occurs in non-premixed conditions. These trends are attributed to the lower complexity of premixed flamelets compared to the ones driven by diffusion. Since mixing is required for the development of non-premixed flamelets, the subgrid dynamics have a greater potential to decoupled from the processes at larger scales, reducing the model’s validity. The dependency of the model’s coefficients on the grid size was studied comparing the fitted values for different subgrid fluxes. For filters widths below the laminar flame thickness at stoichiometric conditions (${\mathsf{\delta }}_{\mathrm{L}0}$), a linear sensitivity of the model coefficients with respect to the filter widths is found. This linearity degrades if filter widths larger than ${\mathsf{\delta }}_{\mathrm{L}0}$ are considered, and a power law in the order $\mathrm{C}~{\mathsf{\Delta }}^{0.36}$ provides a good general fitting. This change is attributed to the simultaneous presence of various flame regions at a subgrid level, which decouples the small- and large-scale dynamics. Overall, the proposed model can predict subgrid turbulent fluxes with high precision for filter widths up to twice the laminar flame thickness at stoichiometric conditions. The validity of the presented models has been investigated in complete Large Eddy Simulations as well. It was found that the introduction of subgrid gradient models at positions sufficiently far away from the walls improves numerical accuracy. Future research works shall apply the proposed methodology at different combustion conditions with different turbulent and chemical dynamics to assess the universality of the presented models.