Next Article in Journal
Specific Relation Attention-Guided Graph Neural Networks for Joint Entity and Relation Extraction in Chinese EMR
Previous Article in Journal
Study on the Influence of Ultrasound Homogenisation on the Physical Properties of Vegan Ice Cream Mixes
Previous Article in Special Issue
A Transition of Ignition Kernel Delay Time at the Early Stages of Lean Premixed n-Butane/Air Turbulent Spherical Flame Propagation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Differential Subgrid Stress Model and Its Assessment in Large Eddy Simulations of Non-Premixed Turbulent Combustion

1
Central Aerohydrodynamic Institute (TsAGI), 1 Zhukovsky Str., 140180 Zhukovsky, Russia
2
Independent Researcher, 91120 Palaiseau, France
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2022, 12(17), 8491; https://doi.org/10.3390/app12178491
Submission received: 20 July 2022 / Revised: 19 August 2022 / Accepted: 20 August 2022 / Published: 25 August 2022
(This article belongs to the Special Issue Advances in Turbulent Combustion)

Abstract

:
We present a new subgrid stress model for the large eddy simulation of turbulent flows based on the solution of transport equations for stress tensor components. The model was a priori term-by-term calibrated against an open DNS database on forced isotropic turbulence (Johns Hopkins University database). After that, it was applied in a large eddy simulation of non-premixed turbulent combustion. To demonstrate the impact of the new subgrid stress model on scalar fields, we excluded the backward effect of heat release on the subgrid stresses, considering an isothermal reaction (i.e., diluted mixture; the density variations associated with chemical heat release can be neglected) and a Burke–Schumann reaction sheet approximation. A periodic box filled with a homogeneous turbulent velocity field and a three-layer top-hat mixture fraction field was studied. Four simulations were performed in which a fixed model for mixture fraction and its variance was combined with either the proposed subgrid stress model or one of the standard models, including Smagorinsky, dynamic Smagorinsky and WALE. Qualitatively correct backscatter was observed in a simulation with the new model. The differences in the statistics of the mixture fraction and reactive component fields caused by the new subgrid stress model were analyzed and discussed. The importance of using an advanced subgrid stress model was highlighted.

1. Introduction

Presently, almost all large eddy simulations (LESs) are carried out with the use of either algebraic or differential subgrid viscosity models [1]. Non-negative subgrid viscosity, ν s g s , prevents these models from capturing the regions of “reverse” energy flux from subgrid scales to resolved ones, or backscatter, observed experimentally [2] and numerically [3]. A dynamic procedure alleviates this problem by allowing negative ν s g s values [4]. However, this approach suffers from numerical instability, and a common way to stabilize a dynamic model is to introduce ν s g s clipping [5] or spatial averaging [6]. A more theoretically justified way to address the issue of backscatter may be to employ a differential subgrid stress model (DSM) in which a differential equation is solved for each component of a subgrid stress tensor. In DSM, the principal axes of stress and strain tensors as well as their eigenvalues are decoupled. This allows the simulation of inefficient and even locally reversed energy cascades, as recently studied by analyzing DNS data [3].
Although there is evidence that the correct choice of subgrid model may be crucial for a successful LES [7,8], there are still few papers on this topic [9,10,11,12,13]. Pioneering research was published by J.W. Deardorff [9], who formulated the first DSM in 1973 and conducted atmospheric turbulence simulations on a 40 × 40 × 20 gird, demonstrating improved results compared with eddy viscosity models. After more than two decades, C. Fureby et al. published an updated version of a DSM [10] and discussed its potential benefits. Since the 2010s, several papers have been published [11,12,13] on hybrid RANS/LES models which combine a differential Reynolds stress model near a wall with a DSM built on the same basis. Thus, research on DSMs is presently entering its active phase.
Subgrid models also determine the scalar dissipation field and the smallest resolved eddy structure. Both factors have an essential influence on the characteristics of non-premixed turbulent flames, controlled by mixing fuel with an oxidizer. Thus, an advanced subgrid stress model can also improve the predictions of turbulent diffusion combustion.
The current paper contributes to the development of a DSM. A systematic approach is presented to choose appropriate closures and calibrate their coefficients using an open DNS database. The resulting model is compared with the classical Smagorinsky model [14], its dynamic version [4] and a popular WALE model [15] in large eddy simulations of turbulent diffusion combustion. To simplify the task and exclude the backward influence of combustion on the turbulence, a case of so-called isothermal combustion [16] was used. Additional simplifications include consideration of a one-stage unidirectional reaction, proceeding at an infinite speed (Burke–Schumann reaction sheet approximation [17]). The average rates of chemical reactions were determined using the classical non-premixed flamelet approach [18]. Estimates of the differences in the statistics of the mixture fraction fields were obtained, analyzed and discussed.
The structure of the paper is as follows. In Section 2, closures for the terms in the subgrid stress tensor equation are proposed, and coefficient calibration is discussed and performed. In Section 3, an LES of homogeneous turbulence in a periodic box is presented. A three-layer initial field of mixture fraction, corresponding to a fuel layer placed between two oxidizer layers, is considered. Mixture composition fields, developing due to isothermal combustion, are obtained at each time moment by processing the instantaneous fields of mixture fraction and dispersion. The new DSM is compared to the standard subgrid viscosity models. In Section 4, the results are discussed, followed by our conclusions.

2. Subgrid Stress Tensor Equation Closure and Calibration

2.1. Subgrid Stress Tensor Equation

We base this research on the incompressible Navier–Stokes equations, as follows:
u i x i = 0 ,
u i t + u k u i x k = 1 ρ p x i + ν 2 u i x k x k .
Summation over repeated indices is implied.
After spatial filtering denoted by an overbar, an additional term called subgrid stress tensor appears in the following momentum equation:
u ¯ i t + u ¯ k u ¯ i x k + u i u k ¯ x k = 1 ρ p ¯ x i + ν 2 u ¯ i x k x k .
Hereinafter, cumulants [19], or central moments [20], are used to abbreviate the following notation:
a b ¯ = a b ¯ a ¯ b ¯ ,           a b c ¯ = a b c ¯ a ¯ b c ¯ b ¯ c a ¯ c ¯ a b ¯ a ¯ b ¯ c ¯ .
By combining (1) and (2) and applying a filtering operator, one can obtain the following differential equation for the subgrid stress tensor:
u i ' u j ' ¯ t + x k [ u i ' u j ' ¯   u ¯ k + u i ' u j ' u k ' ¯ + 1 ρ ( p u i ' ¯ δ j k + p u j ' ¯ δ i k ) ν u i ' u j ' ¯ x k ]                                                                         = u i ' u k ' ¯ u ¯ j x k u j ' u k ' ¯ u ¯ i x k + 1 ρ p ( u i ' x j + u j ' x i ) ¯ 2 ν u i ' x k u j ' x k ¯
The terms in this equation can be interpreted as follows. On the LHS, under the divergence symbol, there are advection by a resolved flow field, turbulent transport due to subgrid-scale velocity fluctuations, turbulent transport due to subgrid-scale pressure fluctuations and molecular transport, respectively. On the RHS, the first two terms represent the subgrid stress tensor production. This is followed by a pressure–strain correlation which acts as redistribution between energy components in the principal frame of u i u j ¯ . The last term in (3) represents the subgrid stress viscous dissipation. Turbulent transport terms, the pressure–strain correlation and the dissipation tensor require closure models.

2.2. Closure

We are not aiming to construct the best possible closure model for unclosed terms in (3). Our goal is rather to optimize the coefficients of some typical form of closures based on DNS data. According to this, we use locally isotropic assumption, leading to
2 ν u i x k u j x k ¯ = 2 3 ε δ i j ,           ε = ν u i x k u j x k ¯ ,
where ε is the subgrid kinetic energy dissipation rate. It is natural to express ε in terms of subgrid kinetic energy, k = u i u i ¯ / 2 , and characteristic filter size, Δ, as follows:
ε = C E k 3 / 2 Δ .
The simplest model assumes C E = c o n s t ; however, due to the importance of accurate representation of the dissipation process, we will consider a dependence of C E on invariants of subgrid stress anisotropy tensor a = ( a i j ) :
C E = C E ( 0 ) + C E ( 1 ) A 2 + C E ( 2 ) A 3 .
Here, a i j , A 2 and A 3 are defined as
a i j = u i u j ¯ k 2 3 δ i j ,           A 2 = a i j a j i ,           A 3 = a i j a j k a k i .
Equation (5) can be considered as Taylor expansion of C E in the corresponding variables.
The pressure–strain correlation tensor, denoted hereinafter as Φ = ( Φ i j ) , is usually expanded in powers of a . The most general form of linear model for Φ is the following [21]:
Φ = C Φ 1 a ε + k [ C Φ 2 S + C Φ 3 ( a S + S a 2 3 tr ( a S ) I ) + C Φ 4 ( Ω a a Ω ) ] ,
with tensor multiplication obeying the matrix rule, ( A B ) i j = A i k B k j . Here, I = ( δ i j ) is a unity tensor; S and Ω are resolved strain rate and rotation tensors, respectively,
S i j = 1 2 ( u ¯ i x j + u ¯ j x i ) ,               Ω i j = 1 2 ( u ¯ i x j u ¯ j x i ) .
Again, the simplest model assumes C Φ 1 , , C Φ 4 to be constants, but an inclusion of dependence on A 2 and A 3 may be considered as means to improve the model.
For the triple-velocity correlation, we use a gradient model modified by means of resolved velocity gradients, the presence of which follows from exact the equation for u i u j u k ¯ [22], as follows:
u i u j u k ¯ = { i j k } [ k 1 / 2 Δ ( C T 1 δ k l + C T 2 a k l ) + C T 3 Δ 2 u ¯ k x l ] u i u j ¯ x l .  
Summation is performed over cyclic permutations of the i, j and k indices:
{ i j k } A i j k = A i j k + A j k i + A k i j .
The coefficients C T 1 , C T 2 , C T 3 can be either constants or functions of A 2 and A 3 .
It remains to consider the pressure–velocity correlation. According to [23], one can approximate it by contracting the triple-velocity correlation tensor:
1 ρ p u i ¯ = C P T u i u m u m ¯ .
There is one empirical coefficient: C P T .

2.3. DNS Database and Filtering

To calibrate the coefficients in the formulated models, a DNS dataset that was representative of typical subgrid-scale turbulence was required. We used a computation of forced isotropic turbulence [24] from the Johns Hopkins Turbulence Database [25] for this purpose. This dataset was originally generated with a pseudo-spectral algorithm in a periodic domain of size L × L × L with uniform grid of 1024 3 points, L = 2 π . In this flow, the Reynolds number, with respect to the Taylor microscale, is equal to 433; the Kolmogorov microscale is equal to 0.46 h, where h is grid spacing, the integral length scale is equal to 0.217 L; the integral time scale is equal to 0.2 T, where T is simulation period after attaining a statistically steady state. All the numbers are nondimensional.
We extracted 11 snapshots separated by 0.1 T and filtered them using a box filter of different sizes, Δ. Three filter sizes were considered: Δ = 64   h , Δ = 128   h and Δ = 256   h . In Figure 1, their positions on a radial kinetic energy spectrum are shown. While Δ = 64   h is close to the bottom part of the inertial interval, two larger filters are significantly inside it, which corresponds to a correct LES.
To process the filtered data effectively, we extracted an “LES grid” from each filtered field with the spatial step h L E S = Δ / 4 . The derivatives of filtered fields on the RHS of (7) and (8) were calculated by a central difference formula on LES grids, for example:
( u ¯ x ) α , β , γ u ¯ α + 1 , β , γ u ¯ α 1 , β , γ 2 h L E S .
A comparison of spatially averaged, scale-to-scale energy flux, P = u i u j ¯   S i j , calculated first using finite differences (10) and after that employing “exact” values of u ¯ i / x j , showed that the error induced by the finite spatial step was within 3%.

2.4. Calibration Method

In each of the closures in (4), (5), (7), (8) and (9), the empirical coefficients enter linearly. The naïve method of their calibration is to obtain coefficient values from an overdetermined linear equation system by the least squares method, with each equation corresponding to a single data point at a single time instance and a single tensor component. Suppose we have a one-dimensional data vector, b , formed in this way, and that its model is b ˜ = C 1 a 1 + + C n a n . Here, a i and b are exact filtered DNS data. Then, the least squares solution of overdetermined linear system is as follows:
C 1 a 1 + + C n a n = b
which, with respect to C 1 , , C n , gives the coefficients which minimize the Euclidean distance b ˜ b 2 , where
x p = ( i | x i | p ) 1 / p .
The corresponding model b ˜ can be interpreted as orthogonal projection of b onto the linear span of { a 1 , , a n } . However, a drawback follows: if the “basis vectors” a 1 , , a n are inadequate to represent the data vector b , the model vector b ˜ will have significantly smaller norm than that of b . In our preliminary tests, for example, we were systematically obtaining the models for Φ in the form of (7), such that Φ ˜ 2 / Φ 2 0.4 0.5 .
Having faced this problem, we developed a calibration strategy, which is—to our knowledge—a new approach. Let us introduce a correlation coefficient between b ˜ and b :
Corr ( b ˜ ,   b ) = ( b ˜ α b ˜ α ) ( b α b α ) b ˜ b ˜ 2   b b 2
where angle brackets denote averaging over whole filtered dataset. Then, we can perform an iteration process. Denoting the coefficient approximation at the mth iteration by C 1 m , , C n m , at the first step, we update their values by moving along the correlation coefficient gradient, as follows:
C ˜ i m + 1 = C i m + K Corr ( b ˜ ,   b ) C i ,
choosing the constant K to empirically maximize the convergence rate by maintaining the stability of the process (typically, K ranges from 0.5 to 10). In (11), we calculate the correlation coefficient derivative on the RHS by the following central difference formula:
Corr ( b ˜ ,   b ) C i Corr ( + ( C i + Δ C ) a i + ,   b ) Corr ( + ( C i Δ C ) a i + ,   b ) 2 Δ C
with coefficient increment, Δ C , equal to 10 4 . At the second step, we multiply the new coefficients by the norm ratio, C i m + 1 = N C ˜ i m + 1 , where N = b 1 / C ˜ 1 m + 1 a 1 + + C ˜ n m + 1 a n 1 , to ensure the correct norm of the model. Here, we prefer the L 1 norm to avoid giving large weights to large data components. This procedure is repeated until the coefficients converge to a stationary state.
The idea of the method is visualized in Figure 2. A correlation coefficient contour is depicted as a function of two model coefficients, C 1 and C 2 . The white solid line corresponds to pairs C = { C 1 , C 2 } , which yield a correct model norm, b ˜ 1 = b 1 . Let the coefficients at the mth iteration be at the C m point on the plot. The point, after moving along the correlation coefficient gradient, is denoted as C ˜ m + 1 . Re-normalization corresponds to returning to the white solid line along the white dashed line (passing through the coordinate frame origin), so the new coefficients are C m + 1 . The iterations converge to point C on the plot, where the coefficient correlation is maximized, while the model norm is correct. To avoid falling to a local maximum, different initial conditions should be considered; however, all our tests revealed unique solutions independently of widely varied initial coefficient sets.
Finally, it is worth noting that, to avoid programming errors, the presented method was implemented independently by two authors. The identity of the results was verified.

2.5. Calibration Results

At the first step, we calibrated the dissipation rate model of the form (4) and (5). In Table 1, the calibrated coefficients are presented for different filter sizes together with their correlation coefficients.
The model correlates well with an exact ε field, yielding Corr of about 0.71–0.86. After averaging the coefficient values for different filter sizes and keeping two decimal digits, we propose the following dissipation model for use:
ε = max { C E ( 0 ) + C E ( 1 ) A 2 + C E ( 2 ) A 3 ,   C E l i m } k 3 / 2 Δ , C E ( 0 ) = 0.8 ,       C E ( 1 ) = 1.4 ,       C E ( 2 ) = 2.1 ,       C E l i m = 0.21 .
Note that the inclusion of A 2 and A 3 is justified because the correlation coefficient of a constant coefficient model (with only C E ( 0 ) non-zero) is approximately 0.11 lower. The maximum function is introduced to avoid possible undershoots of the dissipation rate in strongly anisotropic flow regions. We selected the value of C E l i m by referring to the probability density function of “exact” C E = ε Δ / k 3 / 2 distribution, as calculated from DNS data (see Figure 3). C E values lower than C E l i m = 0.21 are observed, with less than 1% probability.
At the second step, we investigated the constant coefficient model (7). The results are shown in Table 2.
Averaged coefficient values are C Φ 1 = 1.8 , C Φ 2 = 1.3 , C Φ 3 = 0.39 and C Φ 4 = 0.69 . Unfortunately, the preliminary simulations using this model showed that it is unstable due to the value of C Φ 2 being too high. Our attempt to introduce the dependence of C Φ 1 , , C Φ 4 on A 2 and A 3 did not help. The model remained unstable, while its correlation coefficient increased only by 0.2%. After a series of attempts, we found a workable solution, which was to replace the C Φ 2 k S term in the model with the following:
C Φ 2 ( 0 ) k S + C Φ 2 ( 1 ) k Δ ( S Ω Ω S ) + C Φ 2 ( 2 ) k Δ ( Ω 2 1 3 tr ( Ω 2 ) I )                                                                                         + C Φ 2 ( 3 ) Δ 2 ( S Ω Ω + Ω Ω S 2 3 tr ( S Ω Ω ) I tr ( Ω 2 ) S ) .
The introduced nonlinear combinations of S and Ω are widely used in the framework of EARSM [26]. The values of model coefficients obtained during the calibration are collected in Table 3.
Compared with the constant coefficient model, the correlation coefficient increased, on average, by 1%; however, the extended model turned out to produce stable solutions.
After averaging the data in Table 3, we formulated the following model:
Φ = C Φ 1 a ε + C Φ 2 ( 0 ) k S + C Φ 2 ( 1 ) k Δ ( S Ω Ω S ) + C Φ 2 ( 2 ) k Δ ( Ω 2 1 3 tr ( Ω 2 ) I ) + C Φ 2 ( 3 ) Δ 2 ( S Ω Ω + Ω Ω S 2 3 tr ( S Ω Ω ) I tr ( Ω 2 ) S ) + C Φ 3 k ( a S + S a 2 3 tr ( a S ) I ) + C Φ 4 k ( Ω a a Ω ) , C Φ 1 = 3.0 ,       C Φ 2 ( 0 ) = 1.04 ,       C Φ 2 ( 1 ) = 0.136 ,       C Φ 2 ( 2 ) = 0.058 ,         C Φ 2 ( 3 ) = 0.23 ,       C Φ 3 = 0.34 ,       C Φ 4 = 1.2 .
Consider now a triple-velocity correlation model (8). Again, starting with constant coefficients, C T 1 , , C T 3 , we obtained the results collected in Table 4. Interestingly, the corresponding term norms relate as follows:
C T 1 k 1 / 2 Δ { i j k } u i u j ¯ x k 2 : C T 2 k 1 / 2 Δ { i j k } a k l u i u j ¯ x l 2 : C T 3 Δ 2 { i j k } u ¯ k x l u i u j ¯ x l 2 = 6 : 1 : 53 .
The second term in the model is clearly negligible, so we omitted it. The recalibration of the remaining coefficients is presented in Table 5. Fortunately, the high correlation coefficient was not affected by the simplification.
An attempt to make transport coefficients dependent on A 2 and A 3 was unsuccessful, giving no significant improvements. The final model reads as follows:
u i u j u k ¯ = { i j k } [ C T 1 k 1 / 2 Δ δ k l + C T 3 Δ 2 u ¯ k x l ] u i u j ¯ x l ,   C T 1 = 0.019 ,       C T 3 = 0.064 .  
Making the model coefficients dependent on A 2 and A 3 allowed us to increase Corr only by 0.2–0.3%, which is impractical. Note that, according to (15), the leading term in the model (16) is the one containing the resolved velocity gradients. More traditional models, which do not contain explicit dependence on velocity gradients, perform significantly worse. Our tests show that a correlation coefficient for a model with only C T 1 non-zero has a Corr in the range 0.27–0.33, depending on Δ .
Finally, the optimal value of C P T in the pressure velocity correlation model given by (9) was found to be 0.482 for Δ = 64   h , 0.422 for Δ = 128   h and 0.344 for Δ = 256   h . The corresponding model,
1 ρ p u i ¯ = C P T u i u m u m ¯ ,             C P T = 0.42 ,
has correlation coefficient in the range 0.33–0.39, which is not promising. However, making C P T a function of local dimensionless parameters improved Corr by as little as 0.02, so we do not introduce such extensions to the final model.
Equations (12), (14), (16) and (17) form a closure of (3), which we propose in this paper.

2.6. Standard Subgrid Viscosity Models

In the current paper, we compare the performance of the proposed DSM with that of Smagorinsky model [14], dynamic Smagorinsky model [4] and WALE model [15], each extended with diagonal stress terms. These three models express the subgrid stress tensor as follows:
u i u j ¯ = 2 3 k δ i j 2 ν s g s S i j .
In the Smagorinsky model, subgrid viscosity is defined as
ν s g s = ( C s H ) 2 S ,
where H is the computational grid spacing (equal to characteristic filter size Δ ), S = 2 S i j S i j and C s is a constant coefficient. In dynamic Smagorinsky model, ν s g s = C d y n H 2 S , where
C d y n = min { max { 1 α m l S m l H 2 S S ˜ k n S k n α 2 H 2 S ˜ S ˜ p q S p q ,   C d y n m i n } ,   C d y n m a x } ,
i j = u ¯ i u ¯ j ˜ u ¯ ˜ i u ¯ ˜ j .
Tilda denotes averaging over a larger star-shaped filter of seven cells. Its size exceeds that of the original filter by a factor of α = 7 3 . The coefficient is clipped on both sides to stabilize the model.
In the WALE model,
ν s g s = ( C W H ) 2 ( 𝓈 i j d 𝓈 i j d ) 3 / 2 ( S i j S i j ) 5 / 2 + ( 𝓈 i j d 𝓈 i j d ) 5 / 4 ,
where
𝓈 i j d = ( g 2 ) i j + ( g 2 ) j i 2 1 3 ( g 2 ) k k δ i j ,           g i j = u ¯ i x j ,           ( g 2 ) i j = g i k g k j
and C W is a constant coefficient. There are different approaches to determine the coefficient values [27,28]. In the next section, we find C s , C d y n m i n , C d y n m a x and C W from the requirement of obtaining the correct behavior of the energy spectrum in the inertial range.
In subgrid viscosity models, the subgrid kinetic energy can be evaluated by referring to the second invariant, A 2 , of the anisotropy tensor, see Equation (6). According to the filtered DNS data, the volume-averaged A 2 value is 0.31 ± 0.04 for the filter widths chosen in Section 2.3. From (18), A 2 is expressed as 2 ( ν s g s S / k ) 2 . Equating this to 0.31, we obtain
k = ( 2 0.31 ) 1 / 2 ν s g s S 2.54 ν s g s S .
We use this formula in the Smagorinsky and WALE models. In dynamic Smagorinsky model, we restrict the kinetic energy to non-negative values:
k = max { 2.54 ν s g s S ,   0 } .
Finally, note that, when the DSM is employed, we approximate the subgrid viscosity which appears in the mixture fraction equations (see below) by the Smagorinsky formula (19) and determine k according to its definition, k = u i u i ¯ / 2 .

3. Large Eddy Simulations of Isothermal Combustion

3.1. Problem Statement and Numerical Setup

To study the subgrid model’s influence on the mixing and combustion process, we consider a three-layer top-hat mixture fraction field, transported by a homogeneous turbulent velocity field. Four simulations were performed using the new DSM, and the Smagorinsky, dynamic Smagorinsky and WALE models. In all cases, the flow corresponds to a non-premixed turbulent combustion regime without the backward effect of heat release. Since we use a compressible solver, we solve a compressible counterpart of (1) and (2), supplemented by an energy equation (molecular transport omitted due to formal R e limit):
ρ ¯ t + ρ ¯ u ¯ i x i = 0 ,
ρ ¯ u ¯ i t + x k [ ρ ¯ u ¯ i u ¯ k + p ¯ δ i k + ρ ¯ u i u k ¯ ] = 0 ,
ρ ¯ E ¯ t + x k [ ρ ¯ E ¯ u ¯ k + p ¯ u ¯ k + ρ ¯ u i u k ¯ u ¯ i + ρ ¯ u i u i u k ¯ 2 ] = 0 .
This system is closed either by six equations obtained by introducing density into (3), or by algebraic relation (18). The expressions for total energy and state equation are as follows:
E ¯ = u ¯ i u ¯ i 2 + R T ¯ κ 1 ,             p ¯ = ρ ¯ R T ¯ ,
where R 287   J / ( kg K ) , κ = 1.4 . To minimize the compressibility effects, we make sure that the resolved field Mach number, M = ( u ¯ i u ¯ i ) 1 / 2   / ( κ R T ) 1 / 2 , stays less than 0.1 during the simulations (angle brackets denote averaging over the whole computational domain). Accordingly, we neglect the moments involving density fluctuations in (22)–(25).
The equation system is solved numerically using multiblock finite-volume in-house code zFlare [8] with a purely central second-order kinetic energy-preserving scheme [29] and an explicit three-step Heun time-integration method [30]. The time step, Δ t , satisfies the CFL condition, ( u ¯ i + ( κ R T ) 1 / 2 ) Δ t / H 1 , i = 1 ,   2 ,   3 , at each cell.
The computational domain is a triply periodic box, [ L / 2 ,   L / 2 ] 3 , L = 2 π   m . Initially, we fill the domain with a realization of the velocity field, according to [31]. The velocity field is represented as a sum of transverse harmonic waves with von Karman energy spectrum. The initial integral length scale, as determined by the first zero of the transverse correlation function, is chosen to be L 0 = L / 4 , and the initial resolved energy is E k 0 = u i u i ¯ / 2 = 1000   m 2 / s 2 , which leads to initial large eddy turnover time scale of T 0 = L 0 / E k 0 1 / 2 = 0.157   s . For DSM, the initial uniform isotropic subgrid stress field is specified, u i u i ¯ = 2 k 0 δ i j / 3 , where k 0 = 150   m 2 / s 2 . At the first stage of the simulation, the initial field is freely decayed, transforming into a developed turbulent field until the resolved energy drops to E k 1 = 100   m 2 / s 2 . At the second stage, we impose a three-layer mixture fraction field, as follows:
z ¯ = { 1 , L / 6 x L / 6 , 0 , otherwise ,                   z 2 ¯ = 10 9 ,
where z ¯ is a spatially filtered mixture fraction and z 2 ¯ its subgrid-scale variance. We continue the simulation solving two additional equations, as follows:
ρ z ¯ t + x k [ ρ z   ¯ u ¯ k + ρ z u k ¯ ] = 0 ,
ρ z 2 ¯ t + x k [ ρ z 2 ¯   u ¯ k + ρ z 2 u k ¯ ] = ρ P ¯ z ρ χ ¯ .
To ensure that the subgrid stress model comparison is straightforward, we use the same closures of these equations for all considered models. Subgrid transport terms are modeled using a simple gradient hypothesis, as follows:
z u k ¯ = σ z 1 ν s g s z ¯ x k ,                 z 2 u k ¯ = σ z 2 ν s g s z 2 ¯ x k ,
where ν s g s is determined as described in Section 2.6, and σ z 1 = σ z 2 = 1 for simplicity. Accordingly, the production term P ¯ z = 2 z u k ¯ z ¯ / x k is approximated as
P ¯ z = 2 σ z 1 ν s g s z ¯ x k z ¯ x k .
Finally, the scalar dissipation term χ ¯ = 2 ν   z / x k z / x k ¯ is represented as
χ ¯ = C χ ε k z 2 ¯ = C χ C E 0 k Δ z 2 ¯ ,     C χ = 1.5 ,                 C E 0 = 0.47 ,
where the C χ value lies in the range proposed in the literature [32] and C E 0 is determined from filtered DNS by calibrating a constant coefficient dissipation rate model (4).
We intentionally use simple models of z ¯ and z 2 ¯ which can be directly coupled to subgrid viscosity models and DSM. A drawback of this approach is missing backscatter in mixture fraction field (i.e., regions where P ¯ z < 0 ). An advanced model with differential equations for z u k ¯ , which allows the turbulent transport of a mixture fraction to be misaligned with z ¯ / x k , will be investigated in the future.
Spatially filtered values of mass fractions can be found using a subgrid-scale probability density function (PDF) of the mixture fraction z, as follows:
Y j ¯ = 0 1 Y j ( z ) P D F ( z ) d z ,
where Y j ( z ) corresponds to the exact solution of the stationary laminar flamelet equations [18] in the Burke–Schumann approximation [17] (see Table 6).
In Table 6, Y F , Y O and Y P are the mass fractions of the fuel, the oxidizer and the product of a single unidirectional reaction F + O P , respectively, and z s t is the value of the mixture fraction at the stoichiometric surface, where the flame front is located. Let us denote the mass of the oxidizer required for the combustion of unit of mass of the fuel as L 0 . Since the model problem with arbitrary substances is considered, we take L 0 = 1 , which leads to flame front coordinate z s t = 1 / ( 1 + L 0 ) = 0.5 .
A presumed beta PDF is used:
P D F ( z ) = z β 1 ( 1 z ) γ 1 0 1 z ^ β 1 ( 1 z ^ ) γ 1 d z ^ .
The values of β and γ , entering the exponents, can be calculated using the following equations [33,34]:
β = z ¯ ( z ¯ ( 1 z ¯ ) z 2 ¯ 1 ) , ° γ = ( 1 z ¯ ) ( z ¯ ( 1 z ¯ ) z 2 ¯ 1 ) .
To ensure the conditions of β > 0 , γ > 0 , a restriction, z 2 ¯ < z ¯ ( 1 z ¯ ) , is imposed at each time step. If z ¯ = 0 or z ¯ = 1 , the mass fraction values from the corresponding flamelet boundaries ( Y j ( 0 ) and Y j ( 1 ) ,   respectively ) are taken as Y ¯ j . In addition, at very small values of z 2 ¯ , when the beta PDF (29) approaches a delta function δ ( z z ¯ ) , the equation Y ¯ j = Y j ( z ¯ ) is employed instead of (28).
The use of PDF of the form (29) and the Burke–Schumann hypothesis, as well as the absence of the reverse effect of chemical reactions on the turbulent flow, allowed us to analytically calculate the filtered mass fraction values [35] as a postprocessing step after obtaining a solution for z ¯ and z 2 ¯ , as follows:
Y ¯ F = A ( 1 I 1 ) 0 1 z β ( 1 z ) γ 1 d z 0 1 z β 1 ( 1 z ) γ 1 d z + B ( 1 I 2 ) ,               Y ¯ O = 1 + L 0 Y ¯ F ( 1 + L 0 ) z ¯ ,
I 1 = 0 z s t z β ( 1 z ) γ 1 d z 0 1 z β ( 1 z ) γ 1 d z ,     I 2 = 0 z s t z β 1 ( 1 z ) γ 1 d z 0 1 z β 1 ( 1 z ) γ 1 d z ,
where A = ( 1 + L 0 ) / L 0 , B = 1 / L 0 . To calculate the beta-function integrals, library Boost [36] was used in the computational code.

3.2. Subgrid Model Performance

To determine the constants in subgrid viscosity models, we conducted a series of simulations (only the first stages) in a computational domain of 64 × 64 × 64 cells, varying the constants in the range 0.1 C s 0.2 , 0.1 C d y n m i n 0 , 0.05 C d y n m a x 0.2 and 0.3 C W 0.6 . We analyzed the longitudinal energy spectrum E u u ( k x ) at the end of each simulation, as follows:
E u u ( k x ) = | L / 2 L / 2 u ( x ,   y ,   z ) e ı k x d x | 2 y z ,
where y z denotes the averaging over all lines, y ,   z = c o n s t . The values corresponding to the Kolmogorov spectrum, E u u ( k x ) k x 5 / 3 , in the range of 7 k x L / 2 π 24 were selected, which appeared to be equal to C s = 0.13 , C d y n m i n = 0.03 , C d y n m a x = 0.1 and C W = 0.4 . The same procedure was repeated for DSM, in which the ratio Δ / H was varied ( H is the computational grid spacing first appeared in Equation (19)). The optimum choice was found to be Δ = 0.4   H . In Figure 4, the optimal spectra are illustrated for all models along with the resolved kinetic energy decay. The oscillating character of energy is due to small dilatational motions present in the flow field. After a relaxation period of approximately T 0 / 2 , the energy decays according to a power law of E k ~ t 1.1 . All the models give similar statistics. However, this does not mean that the resolved fields are statistically equivalent. The energy decay rate and the energy spectrum correspond to the first two statistical moments of the velocity field. This is insufficient to fully describe the velocity field, e.g., its intermittency is left uncharacterized [3].
Interestingly, despite the discussed similarity in the resolved flow statistics, the subgrid turbulence behaves very differently with these models. At the end of the first stage, the Smagorinsky model predicts the volume-averaged subgrid kinetic energy to be 0.027 E k 1 , the WALE model predicts it to be 0.030 E k 1 , the dynamic Smagorinsky model gives 0.10 E k 1 , while the DSM model gives 0.055 E k 1 . The high value obtained with the dynamic Smagorinsky model was surprising to us. This may indicate that the choice of ν s g s clipping in (20) is somewhat arbitrary and should be performed in some other way; another possible explanation is the use of k clipping in (21), which tends to increase the volume-averaged k value due to the presence of flow regions, where ν s g s < 0 . The DSM predicts the subgrid kinetic energy to be approximately two times higher than the Smagorinsky and WALE models. Most interestingly, DSM simulates backscatter regions where energy flux is locally directed from subgrid scales toward resolved ones which is consistent with experimental [2] and DNS [3] data. The same holds for the dynamic Smagorinsky model, although it overestimates the frequency of backscatter events. In Figure 5, the probability density functions of scale-to-scale energy flux P = u i u j ¯   S i j obtained with all models are plotted and compared with filtered DNS data. DSM is qualitatively correct in this distribution, which is encouraging while negative P values are absent with Smagorinsky and WALE models. The dynamic Smagorinsky model predicts too wide distribution of the energy flux. Note that all models produce approximately the same volume-averaged energy flux which follows from the similarity of lines on the right of Figure 4.
The fact that the energy fluxes are close to each other while subgrid energies differ by a factor of 2 (not taking into account the dynamic Smagorinsky model) may be interpreted from the perspective of [3]. In [3], the authors analyzed the same DNS dataset as in the current paper and concluded that the turbulent cascade is relatively “ineffective” in transporting the energy across the scales. The principal axes of u i u j ¯ and S i j are statistically misaligned, which leads to a cascade efficiency of only 25% in the inertial range. In other words, a major part of the subgrid stresses drives spatial transport instead of participating in the spectral energy transfer. DSM can capture this effect due to the lack of algebraic connection between u i u j ¯ and S i j . On the contrary, subgrid viscosity models align these tensors perfectly, which leads to errors in either energy transfer rate or in the subgrid stress levels.

3.3. Mixture Fraction Statistics

In Figure 6, the fuel-consumption rate, W = d M F / d t , is shown where M F is the total fuel mass obtained by integrating the ρ Y F field over the volume of the computational domain. Consumption rate is normalized by W 0 , an initial value of W in DSM simulation. It is seen that the profiles for all subgrid viscosity models are qualitatively similar while the DSM simulation departs from them near the start of the process and after t / T 0 = 0.05 . The maximum ratio of fuel-consumption rates obtained with different models is as high as 1.8 at t / T 0 0.09 . This indicates the importance of accurate description of subgrid processes in the problem considered.
We selected four time instances, t = { 0.02 T 0 ,   0.05 T 0 ,   0.1 T 0 ,   0.15 T 0 } , for comparative analysis. Below (Figure 7, Figure 8, Figure 9 and Figure 10), we consider the instantaneous z ¯ field (black isoline corresponds to z ¯ = 0.5 ), z ¯ Y Z ( x ) , profile which is obtained by averaging z ¯ over YZ-planes, resolved mixture fractionmixture fraction variance z 2 ¯ r e s ( x ) = z ¯ 2 Y Z z ¯ Y Z 2 , subgrid-scale mixture fraction variance z 2 ¯ s g s ( x ) = z 2 ¯ Y Z and its production rate P ¯ z Y Z ( x ) , and scalar dissipation χ ¯ Y Z ( x ) . The last two statistics are calculated according to (26), (27).
Again, a comparison of the z ¯ Y Z and z 2 ¯ r e s profiles shows that the resolved scalar statistics are similar with all the considered models. However, the subgrid-scale characteristics are significantly different. A typical level of mixture-fraction subgrid-scale variance is two times lower for the DSM than for the Smagorinsky model. A possible explanation of this difference is as follows. Consider a balance between production and dissipation of z 2 ¯ . In local equilibrium, P ¯ z χ ¯ , so
2 σ z 1 ν s g s z ¯ x k z ¯ x k C χ 2.54 ν s g s S Δ z 2 ¯               z 2 ¯ ν s g s S Δ z ¯ x k z ¯ x k Δ H .
In DSM, Δ / H = 0.4 , while in three other models, Δ / H = 1 . This leads to ( z 2 ¯ ) D S M / ( z 2 ¯ ) o t h e r 0.4 , which is even lower than the observed ratio of 0.5. To improve the estimate, one has to take into account the diffusion terms in z 2 ¯ equation and the possible differences in the statistics of the z ¯ / x k field, as predicted by different models.
These differences in the subgrid-scale characteristics may indicate that more advanced turbulent combustion models may behave very differently when employed with different classes of subgrid stress models. The DSM looks promising here because its subgrid stress field is not algebraically connected to the resolved strain field, which is more realistic compared with the subgrid viscosity models. This may be important in large eddy simulations of turbulent combustion with finite-reaction-rate chemistry.

4. Discussion and Conclusions

A new differential subgrid stress model (DSM) for large eddy simulations of turbulent flows was developed, a priori term-by-term calibrated and validated. For this purpose, an open DNS database [24,25] was used. The least squares method was found to be inappropriate due to its tendency to underestimate the modeled field norm. To avoid this problem, a new calibration method was proposed which maximizes the correlation coefficient of a modeled field while maintaining its correct norm. Even with this method, it was found that the calibrated model may appear to be unstable, so an appropriate choice of closure functional form is important in addition to determining the coefficient values. In the present study, only the inclusion of terms which were nonlinear in strain and rotation rates allowed us to stabilize the pressure–strain model. To our knowledge, all well-known pressure–strain models such as SSG [37] are linear in their mean (resolved) velocity gradients, which follows from the exact form of p S i j ¯ in homogeneous turbulence. Thus, the new model can be considered as an extension of classical models to the inhomogeneous turbulence conditions found in the inertial range. Another result which we considered to be unexpected was the importance of the resolved velocity gradients in the triple-velocity correlation model (8). Although it is not a common practice to include such terms, this allowed us to increase the correlation coefficient from 0.3 to 0.85. The idea was inspired by [22], in which the importance of the velocity gradients in u i u j u k ¯ closure was highlighted. The observations made in the current paper may be important for the development of future advanced subgrid stress models.
The model assessment was carried out in an LES of non-premixed turbulent combustion. The impact of the new subgrid stress model on the scalar fields was shown by a comparison of the DSM with the Smagorinsky, dynamic Smagorinsky and WALE models. There appeared to be no significant differences in the resolved velocity and mixture fraction fields, which may indicate that the backward influence of the subgrid scales on the resolved ones was relatively small in the considered problem. However, the new model was able to simulate backscatter in contrast to the Smagorinsky and WALE models, an effect observed in DNS. In the DSM, the backscatter relates to the misalignment between the principal axes of u i u j ¯ and S i j . This misalignment also leads to the relative inefficiency of the turbulent energy cascade [3], which manifests itself as an increased subgrid stress level of DSM compared with the Smagorinsky and WALE models, while the energy transfer rate is approximately the same. The ability to capture such flow features may play a major role in more complex flows where backscatter is a flow-determining effect. Note that the backscatter is also captured by the dynamic Smagorinsky model, but in this case, it is caused by the presence of negative ν s g s values. It was shown that the probability density function of energy flux obtained with the dynamic Smagorinsky model is too wide compared with the filtered DNS data.
Concerning the mixture fraction statistics, while their resolved parts behave similarly independent of a subgrid model, their subgrid parts significantly depend on the stress field. In the conducted simulations, the subgrid scalar variance was found to be two times lower with the DSM than with the Smagorinsky and WALE models. This was due to a smaller production term and a significantly larger scalar dissipation term. We believe that these differences may be important when advanced turbulence–combustion interaction models are employed.

Author Contributions

Conceptualization, V.S. and V.V.; methodology, A.T. and R.B.; programming, R.B., L.U. and A.T.; data analysis and interpretation, A.T., R.B. and V.S. All authors have read and agreed to the published version of the manuscript.

Funding

The research was funded by a grant from the Russian Science Foundation No. 21-71-10105, https://rscf.ru/en/project/21-71-10105/ (accessed on 19 July 2022).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The reference DNS data for model calibration is publicly available at the Johns Hopkins Turbulence Database, http://turbulence.pha.jhu.edu/ (accessed on 19 July 2022).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Sagaut, P. Large Eddy Simulations of Incompressible Flows, 3rd ed.; Springer: Berlin, Germany, 2008; 558p. [Google Scholar]
  2. Higgins, C.W.; Parlange, M.B.; Meneveau, C. Alignment trends of velocity gradients and subgrid-scale fluxes in the turbulent atmospheric boundary layer. Bound.-Layer Meteorol. 2003, 109, 59–83. [Google Scholar] [CrossRef]
  3. Ballouz, J.G.; Ouellette, N.T. Tensor geometry in the turbulent cascade. J. Fluid Mech. 2018, 835, 1048–1064. [Google Scholar] [CrossRef]
  4. Germano, M.; Piomelli, U.; Moin, P.; Cabot, W.H. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 1991, 3, 1760–1765. [Google Scholar] [CrossRef]
  5. Zang, Y.; Street, R.L.; Koseff, J.R. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows. Phys. Fluids A 1993, 5, 3186–3196. [Google Scholar] [CrossRef]
  6. Zhao, H.; Voke, P.R. A dynamic subgrid-scale model for low-Reynolds-number channel flow. Int. J. Numer. Methods Fluids 1996, 23, 19–27. [Google Scholar] [CrossRef]
  7. Montecchia, M.; Brethouwer, G.; Johansson, A.V.; Wallin, S. Taking large-eddy simulation of wall-bounded flows to higher Reynolds numbers by use of anisotropy-resolving subgrid models. Phys. Rev. Fluids 2017, 2, 034601. [Google Scholar] [CrossRef]
  8. Troshin, A.; Bakhne, S.; Sabelnikov, V. Numerical and physical aspects of large-eddy simulation of turbulent mixing in a helium-air supersonic co-flowing jet. In Progress in Turbulence IX. iTi 2021; Örlü, R., Talamelli, A., Peinke, J., Oberlack, M., Eds.; Springer: Cham, Switzerland, 2021; pp. 297–302. [Google Scholar] [CrossRef]
  9. Deardorff, J.W. The use of subgrid transport equations in a three-dimensional model of atmospheric turbulence. J. Fluids Eng. 1973, 95, 429–438. [Google Scholar] [CrossRef]
  10. Fureby, C.; Tabor, G.; Weller, H.G.; Gosman, A.D. Differential subgrid stress models in large eddy simulations. Phys. Fluids 1997, 9, 3578–3580. [Google Scholar] [CrossRef]
  11. Zhuchkov, R.N.; Utkina, A.A. Combining the SSG/LRR-ω differential Reynolds stress model with the detached eddy and laminar-turbulent transition models. Fluid Dyn. 2016, 51, 733–744. [Google Scholar] [CrossRef]
  12. Liu, Y.; Zhou, Z.; Zhu, L.; Wang, S. Numerical investigation of flows around an axisymmetric body of revolution by using Reynolds-stress model based hybrid Reynolds-averaged Navier–Stokes/large eddy simulation. Phys. Fluids 2021, 33, 085115. [Google Scholar] [CrossRef]
  13. Wang, G.; Wang, S.; Li, H.; Fu, X.; Liu, W. Comparative assessment of SAS, IDDES and hybrid filtering RANS/LES models based on second-moment closure. Adv. Mech. Eng. 2021, 13, 1–14. [Google Scholar] [CrossRef]
  14. Smagorinsky, J. General circulation experiments with the primitive equations: I. The basic experiment. Mon. Weather Rev. 1963, 91, 99–164. [Google Scholar] [CrossRef]
  15. Nicoud, F.; Ducros, F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow Turbul. Combust. 1999, 62, 183–200. [Google Scholar] [CrossRef]
  16. Cook, A.W.; Riley, J.J.; Kosály, G. A laminar flamelet approach to subgrid-scale chemistry in turbulent flows. Combust. Flame 1997, 109, 332–341. [Google Scholar] [CrossRef]
  17. Williams, F.A. Combustion Theory, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2018; 708p. [Google Scholar]
  18. Peters, N. Laminar diffusion flamelet models in non-premixed turbulent combustion. Prog. Energy Combust. Sci. 1984, 10, 319–339. [Google Scholar] [CrossRef]
  19. Rota, G.-C. On the combinatorics of cumulants. J. Comb. Theory Ser. A 2000, 91, 283–304. [Google Scholar] [CrossRef]
  20. Germano, M. A proposal for a redefinition of the turbulent stresses in the filtered Navier–Stokes equations. Phys. Fluids 1986, 29, 2323–2324. [Google Scholar] [CrossRef]
  21. Hanjalić, K.; Launder, B. Modelling Turbulence in Engineering and the Environment; Cambridge University Press: Cambridge, UK, 2011; 402p. [Google Scholar]
  22. Younis, B.A.; Gatski, T.B.; Speziale, C.G. Towards a rational model for the triple velocity correlations of turbulence. Proc. R. Soc. Lond. A 2002, 456, 909–920. [Google Scholar] [CrossRef]
  23. Lumley, J.L. Computational modeling of turbulent flows. Adv. Appl. Mech. 1979, 18, 123–176. [Google Scholar] [CrossRef]
  24. Li, Y.; Perlman, E.; Wan, M.; Yang, Y.; Meneveau, C.; Burns, R.; Chen, S.; Szalay, A.; Eyink, G. A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence. J. Turbul. 2008, 9, N31. [Google Scholar] [CrossRef]
  25. Available online: http://turbulence.pha.jhu.edu/ (accessed on 4 July 2022).
  26. Menter, F.R.; Garbaruk, A.V.; Egorov, Y. Explicit algebraic Reynolds stress models for anisotropic wall-bounded flows. In Proceedings of the 3rd European Conference for Aero-Space Sciences, Versailles, France, 6–9 July 2009. [Google Scholar]
  27. Lily, D.K. The representation of small-scale turbulence in numerical simulation experiments. In Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences, Yorktown Heights, NY, USA, 14–16 November 1966; IBM Form No. 320-1951. pp. 195–210. [Google Scholar]
  28. Canuto, V.M.; Cheng, Y. Determination of the Smagorinsky–Lilly constant CS. Phys. Fluids 1997, 9, 1368–1378. [Google Scholar] [CrossRef]
  29. Jameson, A. Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes. J. Sci. Comput. 2008, 34, 188–208. [Google Scholar] [CrossRef]
  30. Lambert, J.D. Computational Methods in Ordinary Differential Equations (Introductory Mathematics for Scientists and Engineers); John Wiley & Sons: Hoboken, NJ, USA, 1973; 278p. [Google Scholar]
  31. Shur, M.L.; Spalart, P.R.; Strelets, M.K.; Travin, A.K. Synthetic turbulence generators for RANS-LES interfaces in zonal simulations of aerodynamic and aeroacoustic problems. Flow Turbul. Combust. 2014, 93, 63–92. [Google Scholar] [CrossRef]
  32. Peters, N. Turbulent Combustion; Cambridge University Press: Cambridge, UK, 2000; 324p. [Google Scholar]
  33. NIST/SEMATECH e-Handbook of Statistical Methods. Available online: https://itl.nist.gov/div898/handbook/eda/section3/eda366h.htm (accessed on 4 July 2022).
  34. Poinsot, T.; Veynante, D. Theoretical and Numerical Combustion, 2nd ed.; RT Edwards: Philadelphia, PA, USA, 2005; 538p. [Google Scholar]
  35. Lien, F.S.; Liu, H.; Chui, E.; McCartney, C.J. Development of an analytical β-function pdf integration algorithm for simulation of non-premixed turbulent combustion. Flow Turbul. Combust. 2009, 83, 205–226. [Google Scholar] [CrossRef]
  36. Available online: https://www.boost.org/ (accessed on 4 July 2022).
  37. Speziale, C.G.; Sarkar, S.; Gatski, T.B. Modelling the pressure–strain correlation of turbulence: An invariant dynamical systems approach. J. Fluid Mech. 1991, 227, 245–272. [Google Scholar] [CrossRef]
Figure 1. Radial kinetic energy spectrum in DNS and different filter cut-off scales.
Figure 1. Radial kinetic energy spectrum in DNS and different filter cut-off scales.
Applsci 12 08491 g001
Figure 2. Iteration process scheme: correlation coefficient contour and b ˜ 1 = b 1 isoline are shown.
Figure 2. Iteration process scheme: correlation coefficient contour and b ˜ 1 = b 1 isoline are shown.
Applsci 12 08491 g002
Figure 3. CE probability density function for Δ = 256 h The region on the left corresponding to 1% of total probability is shown in red.
Figure 3. CE probability density function for Δ = 256 h The region on the left corresponding to 1% of total probability is shown in red.
Applsci 12 08491 g003
Figure 4. Longitudinal spectra at the end of the first stage (left) and resolved kinetic energy evolution during the first stage (right).
Figure 4. Longitudinal spectra at the end of the first stage (left) and resolved kinetic energy evolution during the first stage (right).
Applsci 12 08491 g004
Figure 5. Normalized probability density function of scale-to-scale energy flux.
Figure 5. Normalized probability density function of scale-to-scale energy flux.
Applsci 12 08491 g005
Figure 6. Fuel-consumption rate.
Figure 6. Fuel-consumption rate.
Applsci 12 08491 g006
Figure 7. Time instance t = 0.02T0. Instantaneous z ¯ field obtained with the DSM (upper left), z ¯ Y Z ( x ) (upper right), z 2 ¯ r e s (middle left), z 2 ¯ s g s (middle right), P ¯ z Y Z ( x ) (lower left), χ ¯ Y Z ( x ) (lower right).
Figure 7. Time instance t = 0.02T0. Instantaneous z ¯ field obtained with the DSM (upper left), z ¯ Y Z ( x ) (upper right), z 2 ¯ r e s (middle left), z 2 ¯ s g s (middle right), P ¯ z Y Z ( x ) (lower left), χ ¯ Y Z ( x ) (lower right).
Applsci 12 08491 g007
Figure 8. Time instance t = 0.05T0. Designations as in Figure 7.
Figure 8. Time instance t = 0.05T0. Designations as in Figure 7.
Applsci 12 08491 g008
Figure 9. Time instance t = 0.10T0. Designations as in Figure 7.
Figure 9. Time instance t = 0.10T0. Designations as in Figure 7.
Applsci 12 08491 g009
Figure 10. Time instance t = 0.15T0. Designations as in Figure 7.
Figure 10. Time instance t = 0.15T0. Designations as in Figure 7.
Applsci 12 08491 g010
Table 1. Coefficients in the dissipation rate model (4) and (5).
Table 1. Coefficients in the dissipation rate model (4) and (5).
Filter   Size ,   Δ / h C E ( 0 ) C E ( 1 ) C E ( 2 ) Correlation Coefficient
640.823–1.3081.8480.861
1280.804–1.3802.1030.808
2560.783–1.4212.3050.717
Table 2. Coefficients in the constant coefficient model (7).
Table 2. Coefficients in the constant coefficient model (7).
Filter   Size ,   Δ / h C Φ 1 C Φ 2 C Φ 3 C Φ 4 Correlation Coefficient
642.0201.3640.4490.7910.590
1281.8501.2830.3840.6810.664
2561.6301.2220.3240.5870.747
Table 3. Coefficients in the extended model (7), (13).
Table 3. Coefficients in the extended model (7), (13).
Filter   Size ,   Δ / h C Φ 1 C Φ 2 ( 0 ) C Φ 2 ( 1 ) C Φ 2 ( 2 ) C Φ 2 ( 3 ) C Φ 3 C Φ 4 Correlation Coefficient
643.3051.060–0.134–0.044–0.2690.3561.3350.604
1283.0461.030–0.137–0.061–0.2330.3591.2020.684
2562.7401.029–0.137–0.069–0.1910.3131.0140.760
Table 4. Coefficients in the constant coefficient model (8).
Table 4. Coefficients in the constant coefficient model (8).
Filter   Size ,   Δ / h C T 1 C T 2 C T 3 Correlation Coefficient
640.01780.00840.06150.869
1280.01760.00420.06170.868
2560.02070.00940.06650.844
Table 5. Coefficients in the simplified model (8) without C T 2 .
Table 5. Coefficients in the simplified model (8) without C T 2 .
Filter   Size ,   Δ / h C T 1 C T 3 Correlation Coefficient
640.01850.06190.869
1280.01790.06190.868
2560.02110.06690.844
Table 6. Burke–Schumann solution.
Table 6. Burke–Schumann solution.
z > z s t z < z s t
Y F = z z s t 1 z s t Y F = 0
Y O = 0 Y O = 1 z z s t
Y P = 1 Y F Y P = 1 Y O
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Balabanov, R.; Usov, L.; Troshin, A.; Vlasenko, V.; Sabelnikov, V. A Differential Subgrid Stress Model and Its Assessment in Large Eddy Simulations of Non-Premixed Turbulent Combustion. Appl. Sci. 2022, 12, 8491. https://doi.org/10.3390/app12178491

AMA Style

Balabanov R, Usov L, Troshin A, Vlasenko V, Sabelnikov V. A Differential Subgrid Stress Model and Its Assessment in Large Eddy Simulations of Non-Premixed Turbulent Combustion. Applied Sciences. 2022; 12(17):8491. https://doi.org/10.3390/app12178491

Chicago/Turabian Style

Balabanov, Roman, Lev Usov, Alexei Troshin, Vladimir Vlasenko, and Vladimir Sabelnikov. 2022. "A Differential Subgrid Stress Model and Its Assessment in Large Eddy Simulations of Non-Premixed Turbulent Combustion" Applied Sciences 12, no. 17: 8491. https://doi.org/10.3390/app12178491

APA Style

Balabanov, R., Usov, L., Troshin, A., Vlasenko, V., & Sabelnikov, V. (2022). A Differential Subgrid Stress Model and Its Assessment in Large Eddy Simulations of Non-Premixed Turbulent Combustion. Applied Sciences, 12(17), 8491. https://doi.org/10.3390/app12178491

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop