FLUIDS 3 , 054602 ( 2018 ) Study of subgrid-scale velocity models for reacting and nonreacting flows

A study is conducted to identify advantages and limitations of existing large-eddy simulation (LES) closures for the subgrid-scale (SGS) kinetic energy using a database of direct numerical simulations (DNS). The analysis is conducted for both reacting and nonreacting flows, different turbulence conditions, and various filter sizes. A model, based on dissipation and diffusion of momentum (LD-D model), is proposed in this paper based on the observed behavior of four existing models. Our model shows the best overall agreements with DNS statistics. Two main investigations are conducted for both reacting and nonreacting flows: (i) an investigation on the robustness of the model constants, showing that commonly used constants lead to a severe underestimation of the SGS kinetic energy and enlightening their dependence on Reynolds number and filter size; and (ii) an investigation on the statistical behavior of the SGS closures, which suggests that the dissipation of momentum is the key parameter to be considered in such closures and that dilatation effect is important and must be captured correctly in reacting flows. Additional properties of SGS kinetic energy modeling are identified and discussed.


I. INTRODUCTION
Turbulence is ubiquitous in various applications from aerodynamics, environment, and astrophysics to energy production.As such, turbulence modeling is of primary importance.However, it remains a highly challenging topic as turbulent flows involve a wide range of length and time scales which are strongly coupled and thus yield complex structures.Over the past few decades, large eddy simulation (LES) has become a widely used tool to simulate turbulent flows thanks to the increase in computational power.In its general description, LES is based on the use of a low-pass spatial filter that dissociates the large turbulence scales down to scales associated with the filter size .Only scales down to are resolved, while the residual part requires some form of modeling.Governing equations for LES are derived by applying a spatial filter to the Navier-Stokes equations, which after some manipulations lead to the appearance of these additional, residual terms which are unclosed.The LES approach has been extensively described in the past both for nonreacting [1,2] and reacting flows [3].The residual scales are often referred to as subgrid scales (SGS) when the filter size is assumed to be the size of a mesh element.Although this is not always the case in a LES, this distinction is unnecessary for the purposes of this paper, and thus the terms unresolved and subgrid will be used indistinctly from now on.
The main role of the small scales of turbulence is the dissipation of energy which is often modeled through the use of a category of model based on the eddy-viscosity concept, among which are Both nonreacting and reacting time instants of DNS are considered for the present study.The nonreacting cases are DNS of forced or freely decaying homogeneous isotropic turbulence in a periodic tridimensional box of length Ł 3 = (2π ) 3 discretized using N 3 points on a uniform mesh.A typical configuration of these DNS is presented in Fig. 1(a) where enstrophy isosurfaces are presented for case NR2.These DNS cover a range of Reynolds number based on the Taylor microscale Re λ from 49 to 1131; see Refs.[20,21].Details on the nonreacting DNS conditions can be found in Table I and in the associated references.Cases NR1 and NR2 consider decaying turbulence while case NR3 involves forced turbulence.Note that due to the decaying turbulence, data in Table I for cases NR1 and NR2 have been computed directly from DNS for the particular time instant analyzed here.
Two different reacting cases are considered in the present study.These are both DNS of threedimensional (3D), statistically planar, freely propagating, turbulent, premixed flames in a rectangular parallelepiped.Characteristics of these cases, estimated upstream of the flame, are summarized in Table II.Figure 1(b) shows the configuration under study.The case R1 considers a premixed methaneair flame with an equivalence ratio of 0.8.The reactant temperature is set at 600 K.The chemical mechanism used is the Smooke mechanism [24], which involves 16 species and 25 elementary reactions.Temperature-dependent properties are considered and binary Fickian diffusion is used with constant Lewis numbers for each species.The case R2 uses simplified single-step Arrheniustype reaction.Constant Lewis number is considered with Le = 1.0.The heat-release parameter, τ = (T ad − T 0 )/T 0 , is 4.5, where T ad , T 0 are respectively the adiabatic flame temperature and the reactant temperature.
To assess the subgrid-scale quantities, all these DNS are explicitly filtered (using a densityweighting average for the reacting cases).This filtering is performed within a 3 filter box centered at the filtering point.The filtered quantity of interest, F , is then computed as where x is the location of the considered point, ρ and ρ are the instantaneous and filtered densities, and x is the sample point in the filter subspace.The domain is in theory unbounded.In practice, the filter kernel G assumes very small values or zero values very quickly as the distance from x increases, and thus the integration domain is usually truncated to the sphere centered in x with radius TABLE II.DNS conditions for the reacting cases.s L is the laminar flame speed, δ th = (T b − T u )/|∇T | max is the laminar flame thermal thickness, and Da = L 11 s L /(δ th u ) and Ka = (u /s L ) 3/2 (δ th /L 1 1) 1/2 are the Damköler and Karlovitz numbers respectively.Other symbols are as for Table I .The filter considered here, G, is a Gaussian filter kernel defined as where C is a normalization constant such that G(r) integrates to unity.Numerically, the convolution in Eq. ( 1) is performed in physical space for each point x, by summing over the three-dimensional stencil of size s = [int( / h)] 3 with specific weights g i,j,k .These weights are computed by evaluating G at all points (i,j,k) within the stencil, and they sum to unity.This operation is performed efficiently in MATLAB using the intrinsic function smooth3.
Using the above methodology, the subgrid-scale turbulent kinetic energy, k, is defined as where U is the velocity vector.Note that, by definition, k is not a turbulent kinetic energy since in general the filter of a space-varying quantity has nonzero statistical mean.

B. Subgrid scale velocity modelling
Models for LES of reacting and nonreacting flows often require the estimation of the SGS kinetic energy, k, or the subgrid velocity scale, u = √ 2k/3.This quantity can be computed using a transport equation for k; see, for example, Ref. [1].However, this leads to additional unclosed terms [13,25] and computational cost, and is thus often not the preferred choice in a LES.Models which directly estimate the SGS velocity have been developed in the past mainly for nonreacting flows, and their corresponding model constants were derived assuming homogeneity and isotropy of the flow.In this study, four commonly used models are considered, as follows.

Smallest resolved velocity (SRV) model
This model assumes scale similarity between the smallest resolved scale and the largest residual scale.The main idea behind scale similarity is that a field φ at scale behaves similarly, and thus is proportional, to the same field at another scale.In fluid dynamics, this concept is related to that of energy cascade within the inertial range, and scale similarity in this case implies that the energy transfer across scales is similar to that from the smallest resolved scale to the next smallest resolved scale [26].Information at scale , φ , can be extracted by using a filter of size as φ = φ − φ, but this information would not be available in a LES.Using the scale similarity assumption, it can be extrapolated from the resolved scales from φ = φ − φ and then using proportionality between φ and φ .If φ is the velocity field, the above concept can be expressed as [1] where C p is a coefficient assumed to be 1 here, and the hat is used to indicate a test-filtering operation with a filter width .The test filter used in this study is a Gaussian filter of width = 2 and is density weighted for reacting flows.

054602-4
STUDY OF SUBGRID-SCALE VELOCITY MODELS FOR …

Bardina's model
This model also assumes scale similarity, but of kinetic energy rather than of velocity.Its original formulation [8] was not Galilean invariant.In this study, its Galilean-invariant form, proposed in Ref. [27], is used: where C b = 0.126 [8].

Lilly's model
The SGS velocity scale here is estimated from the SGS viscosity, ν sgs , as [7] where c L = 0.094 and c L = 10.64.While the model as originally proposed involves the constant c L in the denominator, for consistency with the other models, it is convenient to re-express it with the constant c L = 1/c L in the numerator.The Smagorinsky model, , is used in this work to estimate the SGS viscosity, where S ij is the filtered strain tensor and c s is taken to be 0.15.It is worth noticing that, although the Smagorinsky closure is used in LES to model the anisotropic part of the SGS stress tensor, Lilly's model is only related to the SGS viscosity.Moreover, because of the previous substitution, this model as it is used here only assumes that u is proportional to the rate of strain.It will be seen later in this paper that this assumption is incomplete as also the rotation contribution is important for the SGS velocity scale.

Colin's model
This model is derived from an expansion in Taylor series with truncation at the second order of the SRV model of Eq. ( 4).The curl operator is then applied to remove the dilatation effect and thus this model only considers the turbulent part of the SGS kinetic energy.For cases where the filter size coincides with the LES grid size, the SGS velocity scale is then expressed as [28] while a correction factor is applied otherwise.The model constant, C 2 = 2, is derived analytically and optimized to match the mean level of SGS kinetic energy in homogeneous isotropic turbulence [28].
A thorough study showing the sensitivity of these models to factors such as the Reynolds number and the filter width is not present in the literature and thus it is not clear how these factors affect the SGS velocity.Moreover, in the presence of combustion, the dilatation produced by the flame is expected to influence the small scales and it is not clear which is its effect on existing models.Although a posteriori analyses were carried out for LES of jet flames in Refs.[12,13], these analyses could only assess the effect of the SGS velocity modeling indirectly, by comparing measured statistics with LES data.A direct assessment of the ability of the above models to predict the SGS statistical properties under different conditions was thus never performed before.This analysis is performed next using the DNS database described in Sec.II A by comparing various statistical properties of k estimated using the models above to those obtained directly from DNS obtained using Eq.(3).

III. RESULTS
The DNS database for reacting and nonreacting flows presented earlier has been postprocessed in MATLAB to compute the required filtered velocity fields and consequently the SGS kinetic energy.Two filter sizes have been used to filter each DNS case in Tables I and II in order to investigate the statistical behavior of the SGS kinetic energy modeling at different scales.According to Pope's turbulent kinetic energy criterion [1], the SGS kinetic energy should amount to no more than 20% of the turbulent kinetic energy in a well-resolved LES.This criterion is evaluated for all cases shown in this work using data directly from DNS, presented in Table III.The values of the two filter sizes used, 1 and 2 > 1 , are reported for each case in Table IV.As shown in the table, all cases satisfy the Pope's criterion except when the larger filter size is used for cases NR1 and NR2 of Table I, in which case the unresolved kinetic energy amounts to respectively 25% and 30% of the turbulent kinetic energy.These two out of ten cases do not affect the conclusions in this paper and the related results are useful for model performance assessment in the case of coarse LES.It is worth mentioning that for the reacting cases Pope's criterion is to be evaluated in the nonreacting regions of the domain, since the SGS kinetic energy within the flame is affected by dilatation other than turbulence, as will be discussed later.
Prior to analyzing the behavior of the SGS kinetic energy models of Sec.II B, some properties of the SGS kinetic energy from DNS are discussed first.Figure 2 shows PDFs of log k + for both reacting and nonreacting cases and different filter sizes, where k + = k/ k is the SGS kinetic energy normalized using its volume average and + = /δ th in Fig. 2(b) is the filter width normalized with the laminar flame thickness.Note that the PDFs for reacting flows in Fig. 2(b) are constructed for points with c 1 − , where c = ( T − T u )/(T b − T u ) is the Favre-filtered reaction progress variable, and T , T u , and T b are respectively the filtered temperature and its unburned and burned values.The constant = 0.01 is used to only account for the statistical behavior inside the flame region, as the nonreacting behavior is already shown in Fig. 2(a), and to avoid the PDFs being affected by the larger number of samples in the numerical domain associated with nonreacting regions.The PDFs in Fig. 2 are compared to Gaussian PDFs constructed using the same mean and variance.The PDFs from DNS follow closely the Gaussian behavior, in particular for the nonreacting flows, implying that the overall behavior of k is log-normal.Also, the PDF width for reacting flows is larger than for the nonreacting cases (larger standard deviation).This is because k varies across the flame, which will be discussed in more detail in the next paragraphs, and thus a broader range of values is considered for k + for reacting cases.Nevertheless, the PDFs of log k + retain the Gaussian shape, indicating that despite the overall variation across the flame, k is log-normal.
A more quantitative investigation of the nearly Gaussian behavior of log k + can be performed by computing higher order statistics such as skewness and kurtosis of the PDFs shown in Fig. 2. The  skewness of a distribution X is a third-order moment and it is defined as γ , where E[•] is the expectation, and μ and σ are mean and standard deviation of the distribution.The skewness gives a measure of the asymmetry of the distribution, and it is zero for a Gaussian distribution.The kurtosis of X is a fourth-order moment giving a measure of the frequency of occurrence of extreme events, and it is defined as Ku . The kurtosis of a Gaussian distribution is 3, and values higher than 3 imply that the frequency of occurrences of extreme events is higher than for a Gaussian distribution.Skewness and kurtosis for the PDFs of Fig. 2 are shown in Table IV along with mean and standard deviation for completeness.It is worth noticing that μ = 0 despite k + = 1 as log k + = log k + and because of the conditioning on c for the reacting cases, while the higher order statistics are unaffected by the choice of the normalization for k.Skewness and kurtosis for cases NR1 and NR2 are very close to the values 0 and 3 respectively for both filter sizes explored in this work, suggesting that log k + has an almost-perfect Gaussian distribution in homogenous isotropic conditions at least for relatively low turbulence.Some skewness and kurtosis above the Gaussian value are observed for the higher turbulence case, NR3.However, these values may be affected by the limited number of samples used for this case.Some larger deviation from the Gaussian behavior exist for the reacting cases, R1 and R2, as was also observed for Fig. 2(b).A non-negligible positive skewness is observed for the lower turbulence case, R1, which is caused by the "bump" observed for log k + ≈ 2 in Fig. 2(b).This bump is in turn produced by the flame preheat region, which could be demonstrated by excluding the bins with value of progress value below about 0.5.Increasing the filter width or the Reynolds number reduces this effect (the skewness is closer to zero).The kurtosis for the reacting cases also decreases with increasing Reynolds number and filter width.It is not clear at this stage if a further increase of Reynolds number or filter size would produce a further decrease of skewness and kurtosis because of the limited set of data used for this study.Thus, in order not to shift the focus of this paper, further analyses are postponed to future work.Different statistics of k, obtained using the models presented in Sec.II B, are compared next with each other and those obtained from DNS.These comparisons are shown for clarity in the next paragraphs for nonreacting and reacting cases separately.

A. Nonreacting cases
Two filter sizes, /L 11 = 0.2 and /L 11 = 0.5, are used for the DNS cases NR1 and NR2 of Table I.This choice guarantees that is at least four times the Kolmogorov length scale and is mainly dictated by the narrow range of scales available for case NR1; see Table V.Two different filter sizes, /L 11 = 0.04 and /L 11 = 0.1, are used instead for case NR3.This different choice is due to memory limitations in MATLAB but does not affect the conclusions to be presented in this section.

Universality
The mean subgrid kinetic energy from DNS, k DNS , computed as a volume average over the entire domain, is shown in Table V for cases NR1 to NR3 and different filter sizes.Ideal values of the model constants for the four models of Sec.I are also shown.These values are set so that the mean SGS kinetic energy from the model matches the one in the DNS.The additional model shown in the table, localized diffusion-dissipation (LD-D), is discussed below in Sec.III B. The data in Table V show that the original constants for the SRV, Bardina's, and Lilly's models would yield a significant underestimation of k.This underestimation becomes more pronounced as the Reynolds number and the filter size increase.Thus, it appears that the model constants depend on Reynolds number and filter size.The model of Colin et al. [28] yields instead a reasonable approximation of the mean SGS kinetic energy using its original model constant, at least for relatively low Reynolds numbers (cases NR1 and NR2), which is somewhat expected since this constant was derived in Ref. [28] ad hoc to match the mean SGS kinetic energy in homogenous isotropic turbulence.The model constant, however, still increases for the highest Reynolds number and filter size.In order to provide an assessment of the universality of these constant, the overall mean and rms are also shown in Table V.According to the DNS data used, even with an ad hoc tuning of the model constants, an uncertainty of 18% to 22% on the value of the constant exists depending on the model used, the filter size, and the Reynolds number.In Sec.III C, it will be seen that for reacting cases these model constants are further dependent on the temperature.

Statistical behavior
The statistical behavior of the four models of Sec.II B is investigated next for cases NR1 to NR3 of Table I.Scatter plots of SGS kinetic energy from DNS, k DNS , versus modeled SGS kinetic energy, k model , are shown in Figs. 3 and 4 for cases NR1 and NR2 respectively.Case NR3 yields similar statistics and is not shown here for conciseness.The Pearson linear correlation coefficient [29], r, is also provided on each scatter plot.This coefficient, summarized in Table VI, indicates how data are correlated independently of the value of the mean, and thus it provides additional information with respect to the data shown for Table V.
For all models, as seen in Fig. 6, the Pearson coefficient decreases with increasing filter size, indicating a weaker correlation as one would expect.The highest correlation is observed for Bardina's model, while the SRV model has the lowest.It is interesting to notice that Colin's model, which is derived from the SRV model by Taylor expansion and truncation at the second order, yields higher correlations than the latter, and this aspect will be further discussed in Sec.IV.Moreover, Colin's model, which gave the best match in terms of mean SGS kinetic energy and universality of the model constant, is not the best performing model in terms of correlations.Similar relative behaviors between the various models are observed for the higher turbulence condition cases of Table I, case NR2 in Fig. 4, and case NR3.Similar to what was observed for Table V, the Pearson coefficient decreases with increasing filter size, as one would expect, and with increasing Reynolds number for all models but the SRV model, where r increases from NR1 to NR2 case and remains roughly constant for case NR3.This behavior may be due to the assumption of scale similarity in the SRV model which is not fully justified for the lowest turbulence case (NR1).However, Bardina's model, which is also based on the scale similarity assumption, yields relatively high correlations for both cases NR1 and NR2.The difference between these two models, apparently similar, can be understood by performing a Taylor expansion and will be discussed further in Sec.IV.
An interesting method to investigate the statistical behavior of k is through the so called q-q plots.These plots are constructed by sorting two sets of samples of two random variables in ascending order and plotting one against the other in a Cartesian plot.A linear behavior of the resulting plot would indicate that the two random variables have the same distributions (to within a linear scaling), whereas nonlinear behavior indicates that their distribution are different.The q-q plots of k DNS versus k model are shown in Fig. 5 for cases NR1 to NR3 of Table I to assess the performance of the four models of Sec.I on a statistical basis.The additional curve shown for the LD-D model will be discussed later.The SGS kinetic energy for each curve in Fig. 5 is normalized using its mean SGS kinetic energy (averaged over the entire volume), k + = k/ k .All four models of Sec.I result in fairly straight lines around the mean (given by k + = 1) and at larger values of k + , but the curves for the SRV and in particular for Colin's model seems to deviate significantly from the linear behavior for lower values of k + .This effect is more marked when the filter size increases, and deviations from linear behavior are observed also around k + = 1 in this case for the SRV and Colin's models.Among the four models discussed here, Bardina's model seems to maintain a rather good linear behavior for all values of k + , which confirms the strong correlation with the DNS data observed earlier for the scatter plots.
It is interesting to observe that Colin's (red), SRV (black), Lilly (green), and LD-D (dashed blue) curves are very close to each other.This means that the corresponding models yield very similar distributions, as is verified in Fig. 6.
As further assessment for the statistical behavior of the modeled SGS kinetic energy, probability density functions (PDF) of normalized k are shown for cases NR1 to NR3 and two filter sizes in Fig. 6.These PDFs are multiplied by k + in the figure so that equal areas correspond to equal probabilities in the semilog scale, since k + P (k + )d log(k + ) = P (k + )dk + .Plots of k + P (k + ) versus log(k + ) are the same as a plots of the PDF of log(k + ) versus log(k + ).The PDFs from the SRV and in particular Colin's models show a long tail for high k + values compared to the PDF from DNS.The long tail in Colin's model PDF indicates that this model has a larger scatter, which was also observed in the q-q plot of Fig. 5.An excellent agreement with the DNS data is observed for Bardina's model, which is also the only model among the ones discussed here which captures the peak k + location for all conditions explored.

B. A dissipation-diffusion based model for k
The analysis presented in the previous section shows Bardina's model to be very suitable to represent the subgrid kinetic energy on nonreacting flows for a range of different conditions of Reynolds numbers and filter sizes, as long as its model constant is appropriately chosen.Moreover, the analysis reveals that, surprisingly, Colin's model performance in terms of Pearson correlation is somewhat better than in the SRV model, which is the model from which it is derived.Using this information, one might think to use the same approach used in Ref. [28] and apply it to Bardina's model, rather than the SRV model.From Eq. ( 1), by assuming constant density and one-dimensional field for simplicity, the Taylor expansion gives where α 2 is an analytical function of the filter G.We note that the odd terms in the Taylor expansion are zero because of the assumed symmetry of G with respect to x .Combining the equations above  with Eq. ( 5) (for one dimension), one obtains By extending to three dimensions and nonconstant density, the equation above becomes This equation is in the form of two contributions, one proportional to a dissipation of kinetic energy and another proportional to the diffusion.For this reason, in this study the model above will be referred to as LD-D model.The constant α 2 is the same for both contributions and is assigned to unity in this study for simplicity.Furthermore, a global constant C m is introduced for conformity with the other models presented in Sec.II B, and its magnitude will be found next directly using the DNS database.The final form of the LD-D model is then The ideal values of C m for the conditions of Table I are shown in Table V.For nonreacting flows, the ideal value of the LD-D model constant is C m ≈ 0.76, and its rms is about 20%, which is similar to the rms observed for the other models.As for Bardina's model, the value of the constant seems to increase with Reynolds number and filter size.However, some sort of convergence around 0.77 is observed for case NR3.The Pearson linear correlation coefficients for the LD-D model are shown in Table VI.These coefficients are very similar to those observed for Bardina's model as one may expect, as the LD-D model is an approximation of Bardina's model.However, the Pearson coefficient becomes higher than that for Bardina's model for the lower turbulence case (NR1).This behavior, which is similar to the relative behavior observed between the SRV and Colin's model, is interesting in particular for reacting cases, to be explored later, because the flame effect usually results in a decrease of turbulence level across the flame front.
Additional insights can be given by observing the q-q plots in Fig. 5.In particular, one shall notice that the LD-D model follows closely Bardina's model behavior only for the high turbulence case (NR3).For cases NR1 and NR2, the LD-D deviates significantly from Bardina's model and from a linear behavior for low values of k, and in particular its behavior is closer to that of Lilly's model.Similar considerations can be made for the PDFs behavior in Fig. 6, where for the lower turbulence cases (NR1 and NR2) the PDF obtained from the LD-D model is more similar to that from Lilly's model than that from Bardina's.However, it is worth noticing that the Pearson correlation for the LD-D model is significantly better than that for Lilly's model.These properties of the LD-D and Lilly's models are not casual and will be further discussed in Sec.IV.
Overall, the LD-D model performs well for nonreacting flows and shows to be the second best model after Bardina's.However, with respect to Bardina's model, the LD-D model has the advantage that only derivatives, and not test-filter operators, are needed to evaluate the SGS kinetic energy.This characteristic may be relevant for cases where the test filtering introduces additional complexities and errors, like on unstructured grids [30][31][32].

C. Reacting cases
The analysis performed in the previous sections is repeated for the reacting flow cases to further assess advantages and limitations of the SGS velocity modeling presented in Sec.II B. Unlike nonreacting flows where homogeneous isotropic conditions where used, for reacting flows the turbulence and thus the SGS kinetic energy are affected by the flame dilatation effect on the velocity and the flame-turbulence interaction.As a result, k statistically varies across the domain as the reaction progresses.In realistic configurations k varies in space, it is thus of interest to understand if models for k are able to predict accurately this variation.The amount of SGS kinetic energy in a certain cell is the result of three competing effects: the dilatation effect, which tends to increase the value of k (velocity acceleration at subgrid scales); the flame-induced turbulence, produced by the combination of dilatation effect and flame movement; and the suppression of turbulence, due to the increased viscosity in the hot gases.We note that the second effect is expected to be small for the relatively low Re configurations investigated here.
The DNS database for reacting flows used in this study is presented in Table II.Two filter sizes, + = 0.2 and + = 0.5, and + = 0.5 and + = 2.0, are used respectively for cases R1 and R2 of Table II, where + = /δ th is the normalized filter width and δ th is the laminar flame thermal thickness.

Universality
The mean subgrid kinetic energy from DNS, k DNS , computed as average over the entire domain, is shown in Table VII.Ideal values of the model constants for the five k models of Eqs. ( 4)-( 9) are also shown.These values are set so that the mean SGS kinetic energy from the model matches the one in the DNS.In order to assess the effect of combustion on the modeling, these constants are also shown for conditional averaged values for different ranges of the filtered progress variable, c.As for the nonreacting cases, the original value of the model constant in the SRV, Bardina's, and Lilly's models would underpredict k, and thus the ideal value must increase significantly to match k DNS .Colin's model yields more accurate predictions if one consider the average of k over the entire domain.However, within the flame region (identified here as the region where 0.5 < c < 0.85), 054602-14 Colin's model original constant would yield a severe underprediction of k.This suggests that the dilation effect due to combustion on k is important, since Colin's model is designed to exclude the dilation from the SGS kinetic energy modeling.This is further explored later in this paper.For all models, the ideal constant varies with Reynolds number and filter size as observed for the nonreacting cases.An additional variation is observed when moving across the flame, which is due to the effect of temperature and heat release on the turbulence level itself.In particular, for the flame investigated in this study, the Reynolds number is observed to decrease moving from the reactant to the products (the overall effect of the flame is to dissipate turbulence).The ideal value of the constants for all five models decreases from the reactants ( c < 0.05) to the products ( c > 0.95) side, which is consistent with the observation that these constants increase with the Reynolds number, as discussed for Table V in the case of cold flows.Within the flame, k is influenced also by additional effects as discussed earlier, which is why the values of the constants do not decrease monotonically moving from the reactants to the products, as can be observed in Table VII, except in the case of the LD-D model, which shows to have a monotonic variation for all conditions explored here.It is worth noticing that the LD-D model does not exhibit the same variation across the flame as that observed in Bardina's model, from which the LD-D model is derived.The overall mean and rms of the ideal constants of the five models are also shown in Table VII.The different domain regions in Table VII are taken to be independent for the computation of mean and rms quantities.These values are useful as they give an indication of what is the uncertainty associated with the model constant if one had to guess its value for a different configuration.In particular, Lilly's model seems to be the most robust, yielding an rms of about 18% of the mean.The LD-D model has an rms of about 20%, which is significantly lower than that of Bardina's model, from which it is derived.The SRV and in particular Colin's model have instead a relatively large rms.
By the observation of the data for the ideal model constants for both reacting and nonreacting flows (Tables V and VII respectively), it is clear that SGS kinetic energy model constants vary with the turbulence level itself.In this respect, one can say that existing models for k are not universal, as a single value of these constants which can be used for different configurations does not exist.Indeed, FIG. 7. Scatter plots of SGS kinetic energy from DNS vs modeled SGS kinetic energy for case R1 of Table II.The scatter plots are shown for the four models of Sec.I and the LD-D model of Eq. ( 9), and results are presented for two different filter sizes.The quantities r and r f , indicated in the plots, represent the Pearson linear correlation coefficients obtained respectively using the entire domain and the flame region only, indicated here as the region where 0.5 c 0.8.Colors represent the progress variable field, c.
the value itself changes with the level of turbulence which is being "measured" by the model.A model constant with a lower rms has less dependence on the turbulence level, and thus the respective model has less uncertainty.In principle, it is also possible to evaluate the model constant dynamically.However, dynamic computations add computational cost and are not always possible as they require a sufficient range of scales for scale similarity to hold.For this reason, it would be important to have a DNS database of such values for different ranges of Reynolds number, which would serve as guidance for a posteriori (LES) analyses.

Statistical behaviour
Scatter plots for cases R1 and R2 of Table VII, colored by c, are shown in Figs 7 and 8 respectively for the five models of Eqs. ( 4)-( 9) and two filter widths.The "slope" of the scattered data increases with c, which is particularly evident for Bardina's and Lilly's models, suggesting that their respective model constants increase monotonically with the progress variable.This slope increase is consistent with the discussion made for Table VII.The Pearson correlation coefficient, r, decreases with FIG. 8. Scatter plots of SGS kinetic energy from DNS vs modeled SGS kinetic energy for case R2 of Table II.The scatter plots are shown for the four models of Sec.I and the LD-D model of Eq. ( 9), and results are presented for two different filter sizes.The quantities r and r f , indicated in the plots, represent the Pearson linear correlation coefficients obtained respectively using the entire domain and the flame region only, indicated here as the region where 0.5 c 0.8.Colors represent the progress variable field, c. increasing filter size, as one would expect for all models except Colin's.The reason for which is the fact that this model excludes the dilatation effect, and this will become clearer later.The Pearson coefficient also decreases from the lower to the higher turbulence case, although it remains almost constant in the cases of Colin's and LD-D models.Surprisingly, the highest correlation is observed for the LD-D model, which performs even better than Bardina's model from which it is derived, reaching peaks as high as r = 0.99 for certain conditions.The same relative behavior is observed between the SRV and Colin's model and this will be further discussed in Sec.IV.The smallest correlations are observed for the SRV model which is similar to what was observed for the nonreacting cases.
The conditional Pearson coefficient, r f , is also shown in Figs.7 and 8 for only the flame region, approximated here as the region where 0.5 c 0.8.r f is smaller than r for all conditions in the cases of the SRV, Bardina's, and Colin's model.A significant decrease is observed in particular for the latter, consistently with the discussion made for Table VII, which suggests that the dilatation effect on k is non-negligible.This will be further explored in Fig. 11.r f increases instead for Lilly's and LD-D models, indicating that these models are particularly suitable for capturing the subgrid phenomena affecting the kinetic energy across the flame.Overall, the LD-D model is shown to be very robust in terms of correlation for all conditions explored here.Further investigations are presented next.
In Fig. 9, q-q plots (in log scale) of normalized k are shown for cases R1 and R2 of Table I.As for the nonreacting case, all models show a linear behavior around the mean, k + ≈ 1, and a tendency toward a nonlinear behaviour for small values of k + .It is worth remembering that a linear behavior indicates a stronger statistical correlation.The similarity of the SRV and Colin's models, on the one hand, and Bardina's and LD-D models, on the other hand, is now evident from these plots.However, the LD-D model seems to maintain the linear behavior down to lower values of k + and to be closer to Lilly's than to Bardina's model for this range.The reason for this similarity between Lilly's and LD-D will be explained in Sec.IV.Let us note that the behavior of the q-q plots for high values of k + is not statistical because of the limited number of samples, and thus it is not discussed here.
PDFs of k + are shown in Fig. 10.As for the cold flow, the SRV and in particular Colin's model yield a long tail for high k + values compared to the PDF from DNS.This is indicative of a significant high scatter and in particular it suggests that large, unphysical values occur relatively frequently in Colin's model irrespective of the value of the mean.Compared to the cold flow, Bardina's PDF shows significant deviation from the DNS PDF for intermediate values of k.A very good agreement is shown instead for Lilly's model and in particular for the LD-D model proposed in this work.The LD-D model deviates significantly from the behavior of Bardina's model from which it is derived, and seems to be closer to Lilly's model as already observed for the q-q plots.The Taylor's expansion performed for the LD-D model and the drop of the higher order terms has the effect of shifting the peak k + to the right as shown in the figure.The same relative behavior is observed between the SRV and Colin's model.However, the SRV model already overpredicts the peak k + position, resulting in Colin's model to further overpredict this value.It is worth mentioning that the PDFs in Fig. 10 are obtained for points with c 1 − , with = 0.01 to avoid for the larger number of bins associated with c = 0 and c = 1 to hide the behavior of the PDFs across the flame.
As discussed earlier, k varies across the domain due to thermal expansion and the different production mechanisms of kinetic energy in the reacting and nonreacting regions.It is thus enlightening to investigate the variation of k + across the flame using the DNS data, and to assess if SGS modeling is able to capture this variation.Conditional averages of k + , conditioned on c, are computed for the five SGS kinetic energy models of Eqs. ( 4)-( 9) and compared with the conditional average of k + from DNS in Fig. 11 for cases R1 and R2 of Table II.The DNS data shows that the SGS kinetic energy is drastically reduced in the products ( c = 1), which is due to the increased viscosity.At intermediate values of c, however, k + reincreases in case R1 or remains constant in case R2.This is because of the dilatalional effect from the flame which prevents the decay of kinetic energy.This is particularly evident for the lower turbulence condition (case R1), where the dilatation effect is strong compared to that of turbulence as one would expect, resulting in a "bump" with local maximum at about c = 0.4 for both filter sizes tested here.It is worth mentioning again here that k is by definition not a turbulent kinetic energy and the data in Fig. 11 clearly show that the nonturbulence effects must be taken into account.When the initial turbulence increases (case R2), the effect of dilatation on k is smaller compared to that of "pure" turbulence and thus the "bump" is no longer visible.It is interesting to notice that the value of k + at intermediate values of c seems to be independent of the filter size for both cases R1 and R2, suggesting that the proportion of k + produced by the flame dilatation is scale independent.
Comparisons between the conditional k from DNS and from models help to enlighten the SGS modeling behavior across the flame.The scale-similarity models (SRV and Bardina's) are able to well predict the turbulence level on the reactant side ( c = 0) and to capture the correct evolution of k + due to dilatation.However, they both overestimate the location of the peak k + and its magnitude for all conditions explored here.Colin's model fails completely to predict the behavior and magnitude of k + at intermediate values of c, as one would expect since this model excludes the dilatation effect.Moreover, the value of k + predicted for the reactant side is overestimated with respect to DNS.This overestimation somehow compensates for the underestimation at intermediate c values, resulting in the behavior for the ideal model constant seen for Table VII.According to Fig. 11, Colin's model seems justified only for relatively high turbulence and filter sizes, as in this case more turbulence is fed within one filter-width length and the dilatation effect on k + becomes small with respect to that produced by the small-scale eddies.In realistic configurations, one can expect to have both regions of small (quasilaminar regions) and high (shears, swirls, obstacles, etc.) turbulence, and SGS modeling must be able to capture well both conditions.Moreover, in order to satisfy Pope's 80% turbulent kinetic energy criterion [1], the filter size must be small enough and in general, for reacting flows, at least of the order of the laminar flame thickness, and results in Fig. 11 show that the dilatation effect becomes more important for small filter sizes.Among all models investigated here, only Lilly's and the LD-D models predict the correct shape of the conditional k + , with the LD-D model exhibiting an almost perfect match with DNS for small filter sizes and only small deviations at larger filter sizes.Moreover, the LD-D model deviates significantly from Bardina's model from which it is derived and performs instead very similarly to Lilly's model.This aspect is discussed in the next section.By analyzing the results presented in this work, it appears that the LD-D model consistently outperforms the other models presented in this work for all conditions explored, both reacting and nonreacting.This model does not rely on the assumption of scale similarity, although it is derived from a model relying on it.Thus, it has the advantage to be computationally inexpensive, as it uses velocity gradients often already available from the flow in common solvers and is easy to implement compared to scale-similarity based approaches, which require a test filtering.

IV. DISCUSSION
In the previous section, for the nonreacting cases, various observations were made.First, the correlation coefficients of Colin's model is higher than for the SRV model, which may be surprising as Colin's model is based on a truncation of the Taylor expansion of the SRV model.This can be explained by the nature of the scale-similarity assumption which requires an additional filtering operation used in the SRV model.Similar to the Taylor expansion performed in Sec.III B, the filtering operation in the SRV model can be seen as introducing high-order terms and derivatives when evaluating a filtered quantity.These high-order terms have less physical meanings and as the SRV model requires an additional test filtering; its increased deviation with regards to DNS quantities can be explained by the introduction of even more high order quantities.In contrast, Colin's model only conserves the physically meaningful term (which can be interpreted in this case as a diffusion term), ∇ 2 U .We notice that the additional curl operation performed in Eq. ( 7) is irrelevant for nonreacting incompressible flows.An additional problem which comes with test-filtering operation is for unstructured grids (nonexplored here), as in this case the explicit filtering introduces additional noise at the higher orders, as shown in previous studies [30][31][32].Moreover, scale-similarity-based models are justified, at least conceptually, only within the inertial subrange [1], although they are observed to be successful also for flow regimes where this assumption lacks of justification [33].Scale-similarity-based models are thus to be used with care and the LD-D model developed here shows to be an excellent alternative.
It is interesting to observe the better quality, for nonreacting flows, of Bardina's model over the other models, and in particular over the SRV model which is also based on a scale-similarity assumption.Looking at the Taylor expansion of Sec.III B, the key difference between Bardina's and the SRV models can be seen as the presence of a pseudodissipation term [first term on the right-hand side in Eq. ( 9)] when computing the SGS velocity.This term is of particular relevance for the computation of the SGS velocity as the main role of the small scales of turbulence is the dissipation of energy, and it is observed to be of an order of magnitude higher than the diffusion term [second term on the right-hand side in Eq. ( 9)].Thus, although the SRV and Bardina's model look similar, there is a key difference between them which ultimately drives the model performance.The importance of this pseudodissipation term becomes even more evident for reacting flows.In this case it was observed that the LD-D model performance is closer to that of Lilly's model than Bardina's from which it is derived.This similarity can be explained by writing the filtered velocity gradient as a sum of strain and rotation rates, ∇ U = S + , where S and are the rate of strain and rate of rotation tensors respectively.By substituting this into the dissipation term in Eq. ( 9), one obtains The first term on the right-hand side in the equation above is the same one appearing in Eq. ( 6) for Lilly's model after using a Smagorinsky-like model for ν sgs .Thus, these two models would be identical (apart for a model constant) if the rotational rate and the diffusion term in Eq. ( 9) are zero.For nonreacting flows, the rotational term is found to be important as one would expect and this leads to the relatively poor performance of Lilly's model compared to the LD-D model.In reacting flows, the fluid expansion by heat release in the flame has a large influence through its effect on the trace of S ij , S kk .Indeed, the transport equation of vorticity can be written as The first two terms on the RHS can be interpreted as the effect of strain rate on the vorticity.In the flame, the dilatation effect proportional to S kk has a preferential direction in the flame-normal direction and acts as a filter which selectively suppresses vorticity not aligned with the flame normal.As a result, there is a large suppression of turbulence in and across the flame and the turbulence on the product side shows a clear anisotropy [34].This effect becomes more predominant as the turbulence intensity decreases and explains why the strain contribution to the subgrid turbulence becomes much more dominant in the flame compared to the rotational rate.As for nonreacting flows, the additional diffusion term in Eq. ( 9) is observed to be significantly smaller than the dissipation term in reacting flows.Thus, for reacting flows Lilly's and LD-D models capture the same physics and for this reason they have similar statistics.We note that the contribution of the diffusion term in Eq. ( 9), which is small for the conditions investigated in this study, was observed to increase with Reynolds number and filter size; thus it may play an important role in a posteriori analyses both in terms of physical consistency and numerical stability.For this reason, it is advised to retain this term for LES.As for the relative behavior between the SRV and Colin's models discussed for nonreacting flows in the previous paragraph, the poor predictions observed in reacting flows for Bardina's model is explainable by the noise introduced by the higher order terms, which are truncated in the LD-D model.The flame may have an amplification effect on this higher order noise.The scale-similarity assumption may be less justified across the flame because of the narrowing of the kinetic energy spectrum as a consequence of the flame dissipating the turbulence.The LD-D model does not rely on scale similarity and this can explain why it performs significantly better than Bardina's model in reacting flows.For reacting flows, the importance of the dilatation effect during the reaction is also particularly highlighted in the previous section.Indeed, due to the heat release and the presence of reactions, this effect cannot be neglected in SGS velocity models.The LD-D model shows an extremely good correlation with the DNS data, which is highlighted in the conditional average in Fig. 11 where only Lilly's and the LD-D models show good agreement with the DNS data.This is readily explained by the fact that these two models have an explicit dependence on the dissipation which is an extremely important effect inside the flame [34].As Colin's model does not account for that effect, a poor agreement is found.Bardina's and the SRV models show to be able to capture the dilatation effect on k.However, this effect seems to be overestimated and overall not properly captured.The relatively poor performance of the scale-similarity-based models may arise again from the test filtering necessary to estimate these models.As mentioned previously, the test filtering introduces high-order terms which have less direct physical meaning and may actually increase the deviation from the actual physical data.

V. CONCLUSION
An analysis of four common models for SGS kinetic energy was performed in this paper using a DNS database of reacting and nonreacting flows.An additional model, based on dissipation and diffusion of the kinetic energy [LD-D model, Eq. ( 9)], was developed from the observed behavior of the other models for nonreacting flows, and this model is shown to outperform all other models for reacting flows for all the conditions and filter sizes explored in this study.In particular, it is found that only the LD-D and Lilly's models are able to capture the correct evolution of k across the flame.For nonreacting flows, the LD-D model is outperformed only by Bardina's model from which it is derived.However, the LD-D model still compares very well with DNS and has the advantage of not relying on the scale-similarity assumption.
An investigation on the equations for the LD-D model suggest that the dissipation of kinetic energy intrinsically adopted in this model is key for accurately predicting the correct level of SGS kinetic energy.Thus, the LD-D model outperforms the SRV and Colin's models, which are proportional only to the diffusion of momentum.Moreover, both strain and rotational components of the dissipation are important in nonreacting flow turbulence, which is why the LD-D model outperforms Lilly's model, which only accounts for the strain contribution.The reason why, in reacting flows, the LD-D model significantly outperforms Bardina's model from which it is derived is identified to be in the effect of the higher order terms, which are truncated in the LD-D model, and which introduce noise at the small scales in Bardina's model.Moreover, unlike the LD-D model, Bardina's model relies on the scale-similarity assumption which may not be fully justified across the flame.A similar relative behavior is also observed for nonreacting flows between the SRV and Colin's model, the latter also being derived from the former, but not for reacting flows, where Colin's model is observed to significantly underestimate the SGS kinetic energy in the flame.This difference is due to the effect of dilatation, which plays an important role in combustiong flows, and it is excluded in Colin's model but not in the SRV model.
The findings of this study can be summarized in three main points: (i) The model constants in common SGS kinetic energy models depend on both Reynolds number and filter size.The DNS data shown here show that values suggested in the literature lead to underestimation of k and thus these constants need to be recalibrated.A range for these constants with associated uncertainty is suggested in this paper and exploring additional Reynolds number ranges would be beneficial to provide guidance in LES analyses where estimation of these constants from dynamic procedures is not possible.(ii) The effect of dilatation is important in reacting flows and must be included in the modeling.Moreover, the dilatation effect seems to be scale independent and is not captured correctly by scale-similarity-based models.(iii) The dissipation of momentum is the key parameter to consider in SGS kinetic energy modeling and is the main reason why the LD-D model presented here outperforms the other models.
The present study is to be completed with further analyses and in particular two directions will be explored in future works: (i) analysis on unstructured grids, since filter size or gradients on unstructured grids may introduce high-order term noise which may affect the results, and (ii) a posteriori analysis with LES to test the newly developed model on real configurations.

FIG. 1 .
FIG. 1.(a) Isosurface of enstrophy (green) for case NR2 threshold at a value μ + 2σ and (b) isosurface of heat release (red) (threshold at 0.8 of the maximum) for case R2 with isosurface of enstrophy (green).

FIG. 2 .
FIG. 2. PDFs of log 10 k + (in standardized form) from DNS () are compared with standard Gaussian PDF () for nonreacting (a) and reacting (b) cases and different filter widths.

FIG. 3 .
FIG. 3. Scatter plots of SGS kinetic energy from DNS vs modeled SGS kinetic energy for case NR1 of Table I.The scatter plots are shown for the four models of Sec.I and two filter sizes.The Pearson linear correlation coefficient, r, is indicated for each scatter plot.

FIG. 4 .
FIG. 4. Scatter plots of SGS kinetic energy from DNS vs modeled SGS kinetic energy for case NR2 of Table I.The scatter plots are shown for the four models of Sec.I and two filter sizes.The Pearson linear correlation coefficient, r, is indicated for each scatter plot.

TABLE I .
DNS conditions for the nonreacting cases.L 11 is the longitudinal integral length scale, λ is the Taylor microscale, η is the Kolmogorov length scale, and h is the grid spacing. .

TABLE III .
Ratios between SGS kinetic energy and turbulent kinetic energy for cases of Tables I and II are shown for two different filter sizes, 1 and 2 > 1 .

TABLE IV .
First four order moments of the distributions shown in Fig.2for log k + .

TABLE V .
Original and ideal values of the model constants to match the mean SGS kinetic energy from DNS.

TABLE VII .
Original and ideal values of the model constants to match the mean SGS kinetic energy from DNS.