Subgrid-Scale Models for Compressible Large-Eddy Simulations^{⋆}

M. Pino Martín

Department of Aerospace Engineering and Mechanics, University of Minnesota, 110 Union St. SE, Minneapolis, MN 55455, USA pino@aem.umn.edu

Ugo Piomelli

Department of Mechanical Engineering, University of Maryland, College Park, MD 20742, USA ugo@eng.umd.edu

Graham V. Candler

Department of Aerospace Engineering and Mechanics, University of Minnesota, 110 Union St. SE, Minneapolis, MN 55455, USA candler@aem.umn.edu

Communicated by M.Y. Hussaini

Received 12 March 1999 and accepted 11 August 1999

Abstract. An a priori study of subgrid-scale (SGS) models for the unclosed terms in the energy equation is carried out using the flow field obtained from the direct simulation of homogeneous isotropic turbulence. Scale-similar models involve multiple filtering operations to identify the smallest resolved scales that have been shown to be the most active in the interaction with the unresolved SGSs. In the present study these models are found to give more accurate prediction of the SGS stresses and heat fluxes than eddy-viscosity and eddy-diffusivity models, as well as improved predictions of the SGS turbulent diffusion, SGS viscous dissipation, and SGS viscous diffusion.

1. Introduction

Large-eddy simulation (LES) is a technique intermediate between the direct numerical simulation (DNS) of turbulent flows and the solution of the Reynolds-averaged equations. In LES the contribution of the large, energy-carrying structures to momentum and energy transfer is computed exactly, and only the effect of the smallest scales of turbulence is modeled. Since the small scales tend to be more homogeneous and universal, and less affected by the boundary conditions, than the large ones, there is hope that their models can be simpler and require fewer adjustments when applied to different flows than similar models for the Reynolds-averaged Navier–Stokes equations.

⋆ The authors gratefully acknowledge the support from the Air Force Office of Scientific Research, under Grant Nos.

AF/F49620-98-1-0035 (MPM and GVC) and AF/F49620-97-1-0244 (UP), monitored by D.L. Sakell. This work was also sponsored by the Army High Performance Computing Research Center under the auspices of the Department of the Army, Army Research Laboratory cooperative agreement number DAAH04-95-2-0003/contract number DAAH04-95-C-0008, the content of which does not necessarily reflect the position or the policy of the government, and no official endorsement should be inferred. A portion of the computer time was provided by the University of Minnesota Supercomputing Institute.

While a substantial amount of research has been carried out into modeling for the LES of incompressible flows, applications to compressible flows have been significantly fewer, due to the increased complexity introduced by the need to solve an energy equation, which introduces extra unclosed terms in addition to the subgrid-scale (SGS) stresses that must be modeled in incompressible flows. Furthermore, the form of the unclosed terms depends on the energy equation chosen (internal or total energy, total energy of the resolved field, or enthalpy).

Early applications of LES to compressible flows used a transport equation for the internal energy per unit mass ε (Moin et al., 1991; El-Hady et al., 1994) or for the enthalpy per unit mass h (Speciale et al., 1988; Erlebacher et al., 1992). In these equations the SGS heat flux was modeled in a manner similar to that used for the SGS stresses, while two additional terms, the SGS pressure-dilatation Π_{dil} and the SGS contribution to the viscous dissipation ε_{v}, were neglected.

Vreman et al. (1995a, b) performed a priori tests using DNS data obtained from the calculation of a mixing layer at Mach numbers from 0.2 to 0.6. They found that the SGS pressure-dilatation π_{dil} and SGS viscous dissipation ε_{v} are of the same order as the divergence of the SGS heat flux Q_{j}, and that modeling ε_{v} improves the results, especially at moderate or high Mach numbers. They also proposed the use of a transport equation for the total energy of the filtered field, rather than either the enthalpy or the internal energy equations; the same unclosed terms that appear in the internal energy and enthalpy equations are also present in this equation.

Very few calculations have been carried out using the transport equation for the total energy, despite the desirable feature that it is a conserved quantity, and that all the SGS terms in this equation can be cast in conservative form. This equation has a different set of unclosed terms, whose modeling is not very advanced yet. Normand and Lesieur (1992) performed calculations of a transitional boundary layer, and modeled only the SGS heat flux, neglecting all the other terms. Knight et al. (1998) performed the LES of isotropic homogeneous turbulence on unstructured grids and compared the results obtained with the Smagorinsky (1963) model with those obtained when the energy dissipation was provided only by the dissipation inherent in the numerical algorithm. They modeled the SGS heat flux and an SGS turbulent diffusion term, and neglected the SGS viscous diffusion. Comte and Lesieur (1998) proposed the use of an eddy-diffusivity model for the sum of the SGS heat flux and SGS turbulent diffusion, neglecting the SGS viscous diffusion.

In this paper the flow field obtained from a DNS of homogeneous isotropic turbulence is used to compute the terms in the energy equations, and evaluate eddy-viscosity and scale-similar models for their parametrization. We place emphasis on the total energy equation, both because of the lack of previous studies in the terms that appear in it, and because of the desirability of solving a transport equation for a conserved quantity. In the remainder of the paper the governing equations are presented and the unclosed terms are defined. The DNS database used for the a priori tests is described. Finally, several models for the unclosed terms are presented and tested.

2. Governing Equations

To obtain the equations governing the motion of the resolved eddies, we must separate the large from the small scales. LES is based on the definition of a filtering operation: a resolved variable, denoted by an overbar, is defined as (Leonard, 1974)

\overline{f}(x) = ∫_{D} f(x′)G(x, x′; \overline{Δ})dx′,  (1)

where D is the entire domain, G is the filter function, and \overline{Δ} is the filter-width associated with the wavelength of the smallest scale retained by the filtering operation. Thus, the filter function determines the size and structure of the small scales.

In compressible flows it is convenient to use Favre-filtering (Favre, 1965a, b) to avoid the introduction of SGS terms in the equation of conservation of mass. A Favre-filtered variable is defined as f̃ = \overline{ρf}/\overline{ρ}. In addition to the mass and momentum equations, one can choose solving an equation for the internal energy, enthalpy, or total energy. Applying the Favre-filtering operation, we obtain the resolved transport equations

\frac{∂\overline{ρ}}{∂t} + \frac{∂}{∂x_{j}} (\overline{ρ} ũ_{j}) = 0,  (2)
\frac{∂\overline{ρ} ũ_{i}}{∂t} + \frac{∂}{∂x_{j}} (\overline{ρ} ũ_{i}ũ_{j} + \overline{p}δ_{ij} − σ̃_{ji}) = −\frac{∂τ_{ji}}{∂x_{j}} ,  (3)
\frac{∂(\overline{ρ} ε̃)}{∂t} + \frac{∂}{∂x_{j}} (\overline{ρ} ũ_{j}ε̃) + \frac{∂q̃_{j}}{∂x_{j}} + \overline{p}S̃_{kk} − σ̃_{ji}S̃_{ij} = −C_{v} \frac{∂Q_{j}}{∂x_{j}} − Π_{dil} + ε_{v},  (4)
\frac{∂(\overline{ρ} h̃)}{∂t} + \frac{∂}{∂x_{j}} (\overline{ρ} ũ_{j}h̃) + \frac{∂q̃_{j}}{∂x_{j}} − \frac{∂\overline{p}}{∂t} − ũ_{j} \frac{∂\overline{p}}{∂x_{j}} − σ̃_{ji}S̃_{ij} = −C_{v} \frac{∂Q_{j}}{∂x_{j}} − Π_{dil} + ε_{v},  (5)
\frac{∂}{∂t}(\overline{ρ} Ẽ) + \frac{∂}{∂x_{j}} [(\overline{ρ} Ẽ + \overline{p})ũ_{j} + q̃_{j} − σ̃_{ij}ũ_{i}] = −\frac{∂}{∂x_{j}} (γC_{v}Q_{j} + \frac{1}{2}J_{j} − D_{j}) .  (6)

Here ρ is the density, u_{j} is the velocity in the x_{j} direction, p is the pressure, ε = c_{v}T is the internal energy per unit mass, T is the temperature; h = ε + p/ρ is the enthalpy per unit mass; E = ε + u_{i}u_{i}/2 is the total energy per unit mass, and the diffusive fluxes are given by

σ̃_{ij} = 2μ̃S̃_{ij} − \frac{2}{3} μ̃δ_{ij}S̃_{kk},  q̃_{j} = −k̃\frac{∂T̃}{∂x_{j}} ,  (7)

where S_{ij} = \frac{1}{2}(∂u_{i}/∂x_{i} + ∂u_{j}/∂x_{i}) is the strain rate tensor, and μ̃ and k̃ are the viscosity and thermal conductivity corresponding to the filtered temperature T̃.

The effect of the SGSs appears on the right-hand side of the governing equations through the SGS stresses τ_{ij}; SGS heat flux Q_{j}; SGS pressure-dilatation Π_{dil}; SGS viscous dissipation ε_{v}; SGS turbulent diffusion ∂J_{j}/∂x_{j}; and SGS viscous diffusion ∂D_{j}/∂x_{j}. These quantities are defined as

τ_{ij} = \overline{ρ} (ũ_{i}u_{j} − ũ_{i}ũ_{j}),  (8)
Q_{j} = \overline{ρ} (ũ_{j}T − ũ_{j}T̃) ,  (9)
Π_{dil} = \overline{pS_{kk}} − \overline{p}S̃_{kk},  (10)
ε_{v} = \overline{σ_{ji}S_{ij}} − σ̃_{ji}S̃_{ij},  (11)
J_{j} = \overline{ρ} (ũ_{j}u_{k}u_{k} − ũ_{j}ũ_{k}u_{k}) ,  (12)
D_{j} = \overline{σ_{ij}u_{i}} − σ̃_{ij}ũ_{i}.  (13)

The equation of state has been used to express pressure-gradient and pressure-diffusion correlations in terms of Q_{j} and Π_{dil}. It is also assumed that \overline{μ(T)S_{ij}} ≃ μ(T̃)S̃_{ij}, and that an equivalent equality involving the thermal conductivity applies. Vreman et al. (1995b) performed a priori tests using DNS data obtained from the calculation of a mixing layer at Mach numbers in the range 0.2–0.6, and concluded that neglecting the nonlinearities of the diffusion terms in the momentum and energy equations is acceptable.

3. A priori Test

One method to evaluate the performance of models for LES or RANS calculations is the a priori test, in which the velocity fields obtained from a direct simulation are filtered to yield the exact SGS terms, and the filtered quantities are used to assess the accuracy of the parametrization. The database used in this study was obtained from the calculation of homogeneous isotropic turbulence decay.

The Navier–Stokes equations were integrated in time using a fourth-order Runge–Kutta method. The spatial derivatives were computed using an eighth-order accurate central finite-difference scheme. The results have been validated by comparison with the DNS data of Martín and Candler (1998, 1999).

Figure 1. Energy spectrum; ⋄ , location of the filter-widths used in the a priori test;  k^{−5/3} slope;  , DNS.

q^{2} = u_{i}u_{i}, and η is the Kolmogorov length scale.

The simulations were performed on grids with 256^{3} points. The computational domain is a periodic box with length 2π in each dimension. The fluctuating fields were initialized as in Ristorcelli and Blaisdell (1997).

The calculation was performed at a Reynolds number Re_{λ} = u′λ/ν = 50, where λ = ⟨u^{2}⟩^{1/2}/ ⟨(∂u/∂x)^{2}⟩^{1/2} is the Taylor microscale and u′ = (u_{i}u_{i})^{1/2} is the turbulence intensity, and at a turbulent Mach number M_{t} = u′/a = 0.52, where a is the speed of sound. The initial flow field is allowed to evolve for four dimensionless time units τ_{t} = λ/u′, so that the energy spectrum may develop an inertial range that decays as k^{−5/3}, where k is the nondimensional wave number.

The filtered fields were obtained using a top-hat filter, which is defined in one dimension as

\overline{f}_{i} = \frac{1}{2n} (f_{i−n/2} + 2 ∑_{i−n/2+1}^{i+n/2−1} f_{i} + f_{i+n/2}) .  (14)

Various filter-widths \overline{Δ} = nΔ (where Δ is the grid size and n = 4, 8, 16, and 32) were used. Note that the grid resolution is high enough that n = 2 would correspond to a DNS. The location of the various filter cutoffs along the energy spectrum at t/τ_{t} = 6.5 are shown in Figure 1; they cover the decaying range of the spectrum (n = 4), the inertial range (n = 8 and 16), and the energy-containing range (n = 32). With the filters used, respectively, 5%, 15%, 40%, and 70% of the total turbulent kinetic energy resides in the SGSs. The two intermediate values are representative of actual LES calculations, in which the SGS kinetic energy is typically between 15% and 30% of the total energy. A higher percentage of SGS energy in general indicates an under-resolved calculation. In the following, results will be shown for \overline{Δ} = 8Δ, except when evaluating the effect of filter-width.

The accuracy of a model is evaluated by computing the exact term R and its model representation M and comparing the two using the correlation coefficient C(R) and the root-mean-square (rms) amplitudes ⟨(R − ⟨R⟩)^{2}⟩^{1/2} and ⟨(M − ⟨M⟩)^{2}⟩^{1/2}. The correlation coefficient is given by

C(R) = \frac{⟨(R − ⟨R⟩)(M − ⟨M⟩)⟩}{(⟨(R − ⟨R⟩)^{2}⟩⟨(M − ⟨M⟩)^{2}⟩)^{1/2}} ,  (15)

where the brackets ⟨·⟩ denote averaging over the computational volume. A “perfect” model would give a correlation coefficient of 1. In the following, the quantities plotted are made nondimensional using the initial values of ρ, u′, and λ.

4. Models for the Momentum Equation

The SGS stresses (8) are the only unclosed term that appears in the momentum equation. Various types of models have been devised to represent the SGS stresses. Eddy-viscosity models try to reproduce the global exchange of energy between the resolved and unresolved stresses by mimicking the drain of energy associated with the turbulence energy cascade. Yoshizawa (1986) proposed an eddy-viscosity model for weakly compressible turbulent flows using a multiscale direct-interaction approximation method. The anisotropic part of the SGS stresses is parametrized using the Smagorinsky (1963) model, while the SGS energy τ_{kk} is modeled separately:

τ_{ij} − \frac{δ_{ij}}{3} τ_{kk} = −C_{s}^{2}2\overline{Δ}^{2}\overline{ρ} |S̃| (S̃_{ij} − \frac{δ_{ij}}{3} S̃_{kk}) = C_{s}^{2}α_{ij},  τ_{kk} = C_{I}2\overline{ρ} \overline{Δ}^{2}|S̃|^{2} = C_{I}α, (16)

with C_{s} = 0.16, C_{I} = 0.09, and |S̃| = (2S̃_{ij}S̃_{ij})^{1/2}.

Moin et al. (1991) proposed a modification of the eddy-viscosity model (16) in which the two model coefficients were determined dynamically, rather than input a priori, using the Germano identity L_{ij} = T_{ij} − τ̂_{ij} (Germano, 1992), which relates the SGS stresses τ_{ij} to the “resolved turbulent stresses”

L_{ij} = \binom{}{\overline{ρu_{i}} \overline{ρu_{j}}/\overline{ρ}} − ˆ \overline{ρu_{i}} ˆ \overline{ρu_{j}}/ˆ\overline{ρ}, and the subtest stresses T_{ij} = ˆ \overline{ρ} ˘ ũ_{i}u_{j} − ˆ \overline{ρ} ˘ ũ_{i}˘ũ_{j} (where ˘ f̃ = ˆ \overline{ρf}/ˆ\overline{ρ}, and

the hat represents the application of the test filter Ĝ of characteristic width Δ̂ = 2\overline{Δ}) that appear if the filter Ĝ is applied to (3). Moin et al. (1991) determined the model coefficients by substituting (16) into the Germano identity and contracting with S̃_{ij}. In the present paper the contraction proposed by Lilly (1992) to minimize the error in a least-squares sense are used instead. Accordingly, the two model coefficients for the dynamic eddy-viscosity (DEV) model will be given by

C = C_{s}^{2} = \frac{⟨L_{ij}M_{ij}⟩}{⟨M_{kl}M_{kl}⟩},  C_{I} = \frac{⟨L_{kk}⟩}{⟨β − α̂⟩},  (17)

where β_{ij} = −2Δ̂^{2} ˆ\overline{ρ} |˘S̃|(˘S̃_{ij} − δ_{ij} ˘S̃_{kk}/3), M_{ij} = β_{ij} − α̂_{ij}, and β = 2Δ̂^{2} ˆ\overline{ρ} |˘S̃ |^{2}.

Scale-similar models are based on the assumption that the most active SGSs are those closer to the cutoff, and that the scales with which they interact are those immediately above the cutoff wave number (Bardina et al., 1980). Thus, scale-similar models employ multiple operations to identify the smallest resolved scales and use the smallest “resolved” stresses to represent the SGS stresses. Although these models account for the local energy events, they underestimate the dissipation.

Speciale et al. (1988) proposed the addition of a scale-similar part to the eddy-viscosity model of Yoshizawa (1986) introducing the mixed model. In this way, the eddy-viscosity contribution provides the dissipation that is underestimated by purely scale-similar models. This mixed model was also used by Erlebacher et al. (1992) and Zang et al. (1992), and is given by

τ_{ij} − \frac{δ_{ij}}{3} τ_{kk} = C_{s}α_{ij} + A_{ij} − \frac{δ_{ij}}{3} A_{kk},  τ_{kk} = C_{I}α + A_{kk},  (18)

where A_{ij} = \overline{ρ} (˜ũ_{i}ũ_{j} − ˜ ũ_{i}˜ũ_{j}).

Erlebacher et al. (1992) tested the constant coefficient model a priori by comparing DNS and LES results of compressible isotropic turbulence and found good agreement in the dilatational statistics of the flow, as well as high correlation between the exact and the modeled stresses. Zang et al. (1992) compared the DNS and LES results of isotropic turbulence with various initial ratios of compressible to total kinetic energy. They obtained good agreement for the evolution of quantities such as compressible kinetic energy and fluctuations of the thermodynamic variables.

Dynamic model adjustment can be also applied to the mixed model (18), to yield the dynamic mixed model (DMM)

C = \frac{⟨L_{ij}M_{ij}⟩ − ⟨N_{ij}M_{ij}⟩}{⟨M_{lk}M_{lk}⟩},  C_{I} = \frac{⟨L_{kk} − N_{kk}⟩}{⟨β − α̂⟩},  (19)

with B_{ij} = ˆ \overline{ρ} ( ˘ ˜˘ũ_{i}˘ũ_{j} − ˘ ˜ ˘ ũ_{i}˘˜˘ũ_{j}), and N_{ij} = B_{ij} − Â_{ij}.

An issue that requires some attention is the necessity to model separately the trace of the SGS stresses τ_{kk}. Yoshizawa (1986), Moin et al. (1991), and Speciale et al. (1988) proposed a separate model for this term. Erlebacher et al. (1992) conjectured that, for turbulent Mach numbers M_{t} < 0.4 this term is negligible; their DNS of isotropic turbulence confirm this conjecture. Zang et al. (1992) confirmed these results a posteriori: they ran calculations with 0 ≤ C_{I} ≤ 0.066 (the latter value is ten times higher than that predicted by the theory) and observed little difference in the results.

Comte and Lesieur (1998) proposed incorporating this term into a modified pressure P. This leads to the presence of an additional term in the equation of state, which takes the form

P = \overline{ρ} RT̃ + \frac{3γ − 5}{6}τ_{kk};  (20)

for γ = \frac{5}{3} the additional term is zero, and for γ = \frac{7}{5} it might be negligible, unless M_{t} is very large. This observation can be used to explain a posteriori the insensitivity of the LES results to the value of C_{I} discussed by Zang et al. (1992): the SGS stress trace can be approximately incorporated in the pressure with no modification to the equation of state. Another factor may be that both the calculations by Erlebacher et al. (1992) and those by Zang et al. (1992) used mixed models, in which the scale-similar part gave a contribution to the normal SGS stresses. Thus, τ_{kk} is taken into account, at least partially, by the scale-similar contribution.

If the mixed model is used, the trace of the SGS stresses can be parameterized without requiring a separate term. A one-coefficient dynamic mixed model (DMM-1) would be of the form

τ_{ij} = Cα_{ij} + A_{ij},  (21)

with

C = \frac{⟨L_{ij}M_{ij}⟩ − ⟨N_{ij}M_{ij}⟩}{⟨M_{lk}M_{lk}⟩}.  (22)

The models DEV (16)–(17), DMM (18)–(19), and DMM-1 (21)–(22) are evaluated in Figures 2–4. Figure 2(a) shows that the DMM-1 model gives the highest correlation for the diagonal components of the SGS stress tensor; Figure 2(b) shows that neither the eddy-viscosity model nor the two-coefficient mixed model DMM predict the rms of the SGS stresses accurately. The DMM-1 model gives the most accurate prediction among those tested.

Figure 3(a) shows the correlation coefficient for the off-diagonal components of the SGS stress. As in incompressible flows, the eddy-viscosity model gives very poor correlation (near 0.2), while much improved results are obtained with the mixed models. Note that the correlation coefficient for DMM and

Figure 2. A priori comparison of the normal SGS stresses τ_{11}. (a) Correlation coefficient and (b) nondimensional rms magnitude. , Eddy-viscosity model DEV (16)–(17); , two-coefficient mixed model DMM (18)–(19); , one-coefficient mixed model DMM-1 (21)–(22); △ , DNS.

Figure 3. A priori comparison of the off-diagonal SGS stresses τ_{12}. (a) Correlation coefficient and (b) nondimensional rms magnitude. , Eddy-viscosity model DEV (16)–(17); , two-coefficient mixed model DMM (18)–(19) and one-coefficient mixed model DMM-1 (21)–(22); △ , DNS.

Figure 4. Nondimensional rms magnitude of τ_{11} versus filter-width at t/τ_{t} = 6.5. , Eddy-viscosity model DEV (16)–(17); , two-coefficient mixed model DMM (18)–(19); , one-coefficient mixed model DMM-1 (21)–(22); △ , DNS.

DMM-1 overlap in the figure. Figure 3(b) shows the rms of τ_{12}. DEV underpredicts the rms magnitude of the exact term, while DMM and DMM-1 slightly overpredict it.

The coefficient C_{s} remained nearly constant at a value of 0.15 throughout the calculation, consistent with the theoretical arguments (Yoshizawa, 1986). The coefficient of the SGS energy, C_{I}, on the other hand, has a value three times higher than predicted by the theory, consistent with the results of Moin et al. (1991).

Figure 4 shows the rms magnitude of τ_{11} versus the filter-width, at time t/τ_{t} = 6.5. For very small filter-widths (\overline{Δ}/Δ = 4), all the models are accurate, reflecting the capability of dynamic models to turn off the model contribution when the grid becomes sufficiently fine to resolve all the turbulent structures (models with constants assigned a priori, such as the Smagorinsky (1963) model, do not have this characteristic). For \overline{Δ}/Δ = 8, consistent with the results shown above, the one-coefficient mixed model DMM gives the most accurate predictions. For intermediate filter-widths, up to \overline{Δ}/Δ = 16, the best prediction is given by the DMM-1 model; when this filter-width is used the unresolved scales contain a considerable amount of energy, 40%. For \overline{Δ}/Δ = 32, it appears that the DMM model predicts the rms

Figure 5. Comparison of unclosed terms in the energy equations. (a) Nondimensional terms in the internal energy or enthalpy equations and (b) nondimensional terms in the total energy equation. , Divergence of the SGS heat flux,

C_{v} ∂Q_{j}/∂x_{j};  , SGS viscous dissipation ε_{v};  , pressure dilatation Π_{dil};  , divergence of the SGS heat
flux, γC_{v} ∂Q_{j}/∂x_{j};  , SGS turbulent diffusion ∂J_{j}/∂x_{j};  , SGS viscous diffusion ∂D_{j}/∂x_{j}.

magnitude accurately. However, since the DMM model overpredicts the rms significantly for \overline{Δ}/Δ = 8 and 16, the accurate prediction given by DMM for \overline{Δ}/Δ = 32 is a coincidence. When \overline{Δ}/Δ = 32 the SGSs contain a large contribution from the energy-containing eddies (70% of the energy is in the SGS); since \overline{Δ}/Δ = 32 is not in the inertial range the assumptions on which LES modeling is based fail. The same results are found for τ_{12} (not shown).

5. Models for the Energy Equations

Figure 5 compares the magnitude of the unclosed terms appearing in the internal-energy and enthalpy equations (4) and (5), respectively (Figure 5(a)) and in the total energy equation (Figure 5(b)). Unlike in the mixing layer studied by Vreman et al. (1995b), in this flow the pressure dilatation Π_{dil} is negligible, and the viscous dissipation ε_{v} is one order of magnitude smaller than the divergence of the SGS heat flux. In the total energy equation (6), the SGS turbulent diffusion ∂J_{j}/∂x_{j} is comparable with the divergence of the SGS heat flux and the SGS viscous diffusion is one order of magnitude smaller than the other terms. In this section several models for the more significant terms are examined.

5.1. SGS Heat Flux

The simplest approach to modeling the SGS heat flux Q_{j} is to use an eddy-diffusivity model of the form

Q_{j} = −\frac{\overline{ρ} ν_{T}}{Pr_{T}} \frac{∂T̃}{∂x_{j}} = −C \frac{\overline{Δ}^{2}\overline{ρ} |S̃|}{Pr_{T}} \frac{∂T̃}{∂x_{j}} ,  (23)

where C is the eddy-viscosity coefficient that can be either assigned if a model of the form (16) is used, or computed dynamically as in (17). The turbulent Prandtl number Pr_{T} can be also fixed or calculated dynamically according to

Figure 6. A priori comparison of the SGS heat flux Q_{j}. (a) Correlation coefficient and (b) nondimensional rms magnitude.

, Eddy-diffusivity model (23), P r_{T} = 0.7; , eddy-diffusivity model (23), Prandtl number adjusted according to (24); , mixed model (26)–(27); △, DNS.

Pr_{T} = \frac{C⟨T_{k}T_{k}⟩}{⟨K_{j}T_{j}⟩},  (24)

where

T_{j} = −Δ̂^{2} ˆ\overline{ρ} |˘S̃|\frac{∂ ˘T̃}{∂x_{j}} + \overline{Δ}^{2} ˆ\overline{ρ} |S̃|\frac{∂T̃}{∂x_{j}} ,  K_{j} = (\frac{ˆ\overline{ρu_{j}} \overline{ρT}}{\overline{ρ}} ) − \frac{ˆ\overline{ρu_{j}} ˆ \overline{ρT}}{ˆ\overline{ρ}} .  (25)

A mixed model of the form

Q_{j} = −C \frac{\overline{Δ}^{2}\overline{ρ} |S̃|}{Pr_{T}} \frac{∂T̃}{∂x_{j}} + \overline{ρ} ( ˜ũ_{j}T̃ − ˜ ũ_{j} ˜T̃)  (26)

was proposed by Speciale et al. (1988). The model coefficients C and Pr_{T} can again be assigned or adjusted dynamically according to (19) and

Pr_{T} = C \frac{⟨T_{k}T_{k}⟩}{⟨K_{j}T_{j}⟩ − ⟨V_{j}T_{j}⟩},  (27)

with

V_{j} = ˆ \overline{ρ} ( ˘ ˜˘ũ_{j} ˘T̃ − ˘ ˜ ˘ ũ_{j} ˘ ˜ ˘T̃) − ˆ \overline{ρ} ( ˜ũ_{j}T̃ − ˜ ũ_{j} ˜T̃).  (28)

Figure 6(a) shows the correlation coefficient for the three models described above. Both eddy-viscosity models overlap on the plot giving a poor correlation factor, roughly 0.2, whereas the mixed model gives a correlation above 0.6. Both eddy viscosity models under-predict the rms of the exact Q_{j}, shown in Figure 6(b), while the mixed model is more accurate. The mixed model maintains accuracy for all filter-widths \overline{Δ}/Δ ≤ 16 (Figure 7).

5.2. SGS Viscous Dissipation

The other term in the enthalpy or internal energy equations that was found to be significant in the present flow is the viscous dissipation ε_{v}. In this section the three models proposed by Vreman et al. (1995b) are tested:

Figure 7. Nondimensional rms magnitude of Q_{j} versus filter-width at t/τ_{t} = 6.5. , eddy-diffusivity model (23), P r_{T} = 0.7; , eddy-diffusivity model (23), Prandtl number adjusted according to (24); , mixed model (26)–(27); △ , DNS.

ε_{v}^{(1)} = C_{ε1} ( ˜σ̃_{ji}S̃_{ij} − ˜ σ̃_{ij} ˜S̃_{ij});  (29)
ε_{v}^{(2)} = C_{ε2}\overline{ρ} q^{3}/\overline{Δ},  q^{2} ∼ \overline{Δ}^{2}|S̃|^{2};  (30)
ε_{v}^{(3)} = C_{ε3}\overline{ρ} q^{3}\overline{Δ},  q^{2} ∼ ˜ ũ_{k}ũ_{k} − ˜ ũ_{k}˜ũ_{k}.  (31)

The first is a scale-similar model; the second and third represent the SGS dissipation as the ratio between the cube of the SGS velocity scale, q̃, and the length scale. The velocity scale can be obtained using either the Yoshizawa (1986) model, as in (30), or the scale-similar model as in (31). Vreman et al. (1995b) fixed the values of the coefficients by matching the rms magnitude of the modeled and exact terms obtained from the a priori test, and obtained C_{ε1} = 8, C_{ε2} = 1.6, and C_{ε3} = 0.6. In the present study the dynamic procedure will be used instead to determine the coefficients. The analog of the Germano identity for this term reads

⟨ ˆσ̃_{ji}S̃_{ij} − ˆ \overline{ρσ_{ij}} ˆ \overline{ρS_{ij}}/ˆ\overline{ρ}^{2}⟩ = ⟨E_{v}^{(n)} − ε̂_{v}^{(n)}⟩ ,  (32)

and the modeled terms ε_{v}^{(n)} can be given respectively by (29)–(31), while the E_{v}^{(n)} are

E_{v}^{(1)} = C_{ε1} ( ˘ ˜˘σ̃_{ji} ˘S̃_{ij} − ˘ ˜ ˘ σ̃_{ij} ˘ ˜ ˘S̃_{ij});  (33)
E_{v}^{(2)} = C_{ε2} ˆ\overline{ρ} ˘ q̃^{3}/Δ̂,  ˘  q̃^{2} ∼ Δ̂^{2}|˘S̃|^{2};  (34)
E_{v}^{(3)} = C_{ε3} ˆ\overline{ρ} ˘ q̃^{3}/Δ̂,  ˘  q̃^{2} ∼ ˘ ˜ ˘ ũ_{i}˘ũ_{j} − ˘ ˜ ˘ ũ_{i}˘˜˘ũ_{j}.  (35)

Figure 8(a) shows the correlation coefficient for the three models. The scale-similar model gives the highest correlation. The use of a velocity scale obtained from the scale-similar assumption, however, results in improved prediction of the rms magnitude; using q ∼ \overline{Δ}|S̃| yields a significant overprediction of the rms. The values of the coefficients obtained from the dynamic adjustment in this flow are significantly lower than those obtained in the mixing layer by Vreman et al. (1995b). For the particular filter-width shown, we obtained C_{ε1} = 2.4, and C_{ε2} = 0.03, while C_{ε3} increased monotonically in time from 0.25 to 0.4. The fact that with these values the first and third models match the rms magnitude of the exact term indicates a lack of universality of these constants. Dynamic adjustment of the model coefficient appears to be beneficial for this term.

The modeling of the viscous dissipation is more sensitive than the other terms to the filter-width. The prediction accuracy deteriorates with increasing filter-width, and in this case even for \overline{Δ}/Δ = 16 none of the models is particularly accurate (Figure 9).

Figure 8. A priori comparison of the SGS viscous dissipation ε_{v}. (a) Correlation coefficient and (b) nondimensional rms

magnitude.  , scale-similar model (33);  , dynamic model (34);  , dynamic model (35); △ , DNS.

Figure 9. Nondimensional rms magnitude of ε_{v} versus filter-width at t/τ_{t} = 6.5.  , scale-similar model (33);  ,

dynamic model (34); , dynamic model (35); △ , DNS.

5.3. SGS Turbulent Diffusion

The SGS turbulent diffusion ∂J_{j}/∂x_{j} appears in the total energy equation (6). Comte and Lesieur (1998) did not model this term explicitly, but added it to the SGS heat flux by using an eddy-diffusivity model to parametrize

(ρ̃Eu_{j} + p̃u_{j}) − (\overline{ρ} Ẽũ_{j} + \overline{p}ũ_{j}) = γ\overline{ρ} (ũ_{j}T − ũ_{j}T̃) + J_{j} ≃ −\frac{ν_{T}}{Pr_{T}} \frac{∂T̃}{∂x_{j}};  (36)

with this model, however, the SGS turbulent diffusion J_{j}, which depends mostly on the unresolved velocity fluctuations, is modeled in terms of the temperature gradient. In an isothermal flow, J_{j} may be nonzero, and, even if the temperature is not constant, there is no reason to couple a term due to mechanical energy gradients to the temperature. A model of the form (36) effectively neglects J_{j}.

The only attempt to model the SGS turbulent diffusion was that by Knight et al. (1998). They argue that ũ_{i} ≃ ˜ ũ_{i} and propose a model of the form

J_{j} ≃ ũ_{k}τ_{jk}.  (37)

A dynamic scale-similar model can be obtained using the generalized central moments (Germano, 1992)

τ(u_{i}, u_{j}) = \overline{ρ} [ũ_{i}u_{j} − ũ_{i}ũ_{j}] ,  (38)
τ(u_{i}, u_{j}, u_{k}) = \overline{ρ} ũ_{i}u_{j}u_{k} − ũ_{i}τ(u_{j}, u_{k}) − ũ_{j}τ(u_{i}, u_{k}) − ũ_{k}τ(u_{i}, u_{j}) − \overline{ρ} ũ_{i}ũ_{j}ũ_{k}.  (39)

Using this notation the turbulent diffusion term can be written as

2J_{j} = τ(u_{j}, u_{k}, u_{k}) + 2ũ_{k}τ(u_{j}, u_{k}),  (40)

since τ(u_{j}, u_{k}) = τ_{jk}. Using this formalism, scale-similar models can be derived by approximating the quadratic terms using the filtered velocities ũ_{j} to replace the velocities u_{j}; for instance, one can write

τ(u_{i}, u_{j}) ∼ τ(ũ_{i}, ũ_{j}) ⇒ \overline{ρ} (ũ_{i}u_{j} − ũ_{i}ũ_{j}) ∼ \overline{ρ} ( ˜ũ_{i}ũ_{j} − ˜ ũ_{i}˜ũ_{j}) .  (41)

If the proportionality constant in (41) is set to one, the scale-similar part of the mixed model (18) is obtained. Analogously, the triple product can be written as

2J_{j} = τ(u_{j}, u_{k}, u_{k}) + 2ũ_{k}τ(u_{j}, u_{k})
≃ C_{J}τ(ũ_{j}, ũ_{k}, ũ_{k}) + 2ũ_{k}τ(u_{j}, u_{k})
= C_{J} [\overline{ρ} ˜ ũ_{j}ũ_{k}ũ_{k} − \overline{ρ} ˜ ũ_{j} ˜ũ_{k}˜ũ_{k} − ˜ ũ_{j}A_{kk} − 2˜ũ_{k}A_{jk}] + 2ũ_{k}τ_{jk},  (42)

the last term is parametrized by the same model used in the momentum equation. The coefficient C_{J} can be set using the identity

ˆ\overline{ρ} ũ_{j}ũ_{k}ũ_{k} − ˆ \overline{ρ} ˘ ũ_{j} ˘ũ_{k} ˘ũ_{k} = 2J_{j} − 2Ĵ_{j},  (43)

where

2J_{j} = C_{J} [ˆ\overline{ρ} ˘ ˜ ˘ ũ_{j} ˘ũ_{k} ˘ũ_{k} − ˆ \overline{ρ} ˘ ˜ ˘ ũ_{j} ˘˜ ˘ũ_{k} ˘˜ ˘ũ_{k} − ˘ ˜ ˘ ũ_{j}B_{kk} − 2˘˜˘ũ_{k}B_{jk}] + 2˘ũ_{k}T_{jk},
B_{jk} = ˆ \overline{ρ} ( ˘ ˜˘ũ_{j} ˘ũ_{k} − ˘ ˜ ˘ ũ_{j} ˘˜ ˘ũ_{k}) ,  (44)

to yield

C_{J} = \frac{⟨( ˆ\overline{ρ} ũ_{j}ũ_{k}ũ_{k} − ˆ \overline{ρ} ˘ ũ_{j} ˘ũ_{k} ˘ũ_{k}) P_{j} − Q_{j}P_{j}⟩}{⟨P_{k}P_{k}⟩},  (45)

where

P_{j} = [ˆ\overline{ρ} ˘ ˜ ˘ ũ_{j} ˘ũ_{k} ˘ũ_{k} − ˆ \overline{ρ} ˘ ˜ ˘ ũ_{j} ˘˜ ˘ũ_{k} ˘˜ ˘ũ_{k} − ˘ ˜ ˘ ũ_{j}B_{kk} − 2˘˜˘ũ_{k}B_{jk}]
− [ ˆ\overline{ρ} ˜ ũ_{j}ũ_{k}ũ_{k} − ˆ \overline{ρ} ˜ ũ_{j} ˜ũ_{k}˜ũ_{k} − ˆ ˜ ũ_{j}A_{kk} − 2 ˆ˜ũ_{k}A_{jk}] ,  (46)
Q_{j} = 2 (˘ũ_{k}T_{jk} − ˆ ũ_{k}τ_{jk}) .  (47)

Figure 10(a) shows the correlation coefficient for the two models (37) and (42) and using (21)–(22) to model τ_{jk}. The correlation factor is greater than 0.7 for both models, and both models overpredict slightly the rms magnitude of J_{j} (Figure 10(b)). When the one-coefficient, scale-similar model is used this overprediction is significantly reduced. Both models perform equally well for \overline{Δ}/Δ ≤ 16, while neither is accurate for \overline{Δ}/Δ = 32.

Subgrid-Scale Models for Compressible Large-Eddy Simulations  13

Figure 10. A priori comparison of the SGS turbulent diffusion J_{j}. (a) Correlation coefficient and (b) nondimensional rms

magnitude.  , knight et al.(1998);  , scale-similar, one-coefficient model; △ , DNS.

Figure 11. A priori comparison of the SGS viscous diffusion D_{j}. (a) Correlation coefficient and (b) nondimensional rms magnitude. , scale-similar model; △ , DNS.

5.4. SGS Viscous Diffusion

The SGS viscous diffusion ∂D_{j}/∂x_{j} is the smallest of the terms in the total energy equation, and is about 5% of the divergence of Q_{j}. No model for this term has been proposed in the literature to date. One possibility is to parametrize it using a scale-similar model of the form

Figure 12. A priori comparison of the sum of the SGS terms in the total energy equation (6). (a) Correlation coefficient and (b) nondimensional rms magnitude. , Model; △ , DNS.

D_{j} = C_{D}( ˜σ̃_{ij}ũ_{i} − ˜ σ̃_{ij} ˜ũ_{i}),  (48)

in which the coefficient can be obtained from

C_{D} = \frac{⟨[ ˆ\overline{ρσ_{ij}} \overline{ρu_{i}}/\overline{ρ}^{2} − ˆ \overline{ρσ_{ij}} ˆ \overline{ρu_{i}}/ˆ\overline{ρ}^{2}] R_{j}⟩}{⟨R_{k}R_{k}⟩},  (49)

where

R_{l} = ( ˘ ˜˘σ̃_{lk} ˘ũ_{k} − ˘ ˜ ˘ σ̃_{lk} ˘˜ ˘ũ_{k}) − ( ˆ ˜σ̃_{lk}ũ_{k} − ˆ ˜ σ̃_{lk}˜ũ_{k}) .  (50)

As can be seen from Figure 11, this model gives a poor correlation and poor agreement for the prediction of the rms magnitude. However, since the viscous diffusion is relatively small, its contribution to the total energy spectrum does not go to the inertial range, but rather to the decaying range. In this situation the accuracy of the model is degraded, as shown by Meneveau and Lund (1997). Thus, the scale-similar approach may still give good predictions when this term is significant. In this particular flow, the error given by the model (or by not using a model) may be tolerable given the small contribution that the term gives to the energy budget.

5.5. General Considerations

In addition to the term-by-term comparisons shown before, it is possible to evaluate the global accuracy of the models by comparing the sum of the exact SGS terms and the modeled quantity, namely,

E_{SGS} = γC_{v}Q_{j} + \frac{1}{2}J_{j} − D_{j}.  (51)

The mixed model (26)–(27) was used for the SGS heat flux, the scale-similar model (44)–(45) for the SGS turbulent diffusion, and the SGS viscous diffusion has been neglected. Figure 12(a) shows the correlation coefficient for the exact and modeled quantities. While the individual correlations were roughly 0.6 and 0.7 for the SGS heat flux model and the SGS turbulent diffusion, respectively, the global correlation drops just below 0.6 when considering the sum of the terms. Figure 12(b) shows the rms for both quantities. The agreement between the exact and modeled quantities is slightly less accurate than for the SGS heat flux alone, Figure 6(b), but more accurate than for the SGS turbulent diffusion alone, Figure 10(b). Figure 12 shows that the overall performance is very good.

6. Conclusions

Several mixed and eddy-viscosity models for the momentum and energy equations have been tested. The velocity, pressure, density, and temperature fields obtained from the DNS of homogeneous isotropic turbulence at Re_{λ} = 50, M_{t} = 0.52 were filtered and the unclosed terms in the momentum, internal energy, and total energy equations were computed.

In the momentum equation, mixed models were found to give better prediction, in terms of both correlation and rms amplitude, than the pure eddy-viscosity models. The dynamic adjustment of the model coefficient was beneficial, as already observed by Moin et al. (1991).

In the internal energy and enthalpy equations only the divergence of the SGS heat flux was significant in this flow; the SGS pressure dilatation Π_{dil} and viscous dissipation ε_{v}, which were significant in the mixing layer studied by Vreman et al. (1995b), were found to be negligible here. Once again, mixed dynamic models gave the most accurate results.

In the total energy equation two additional terms are present, one of which, the turbulent diffusion ∂J_{j}/∂x_{j}, is significant. The model proposed by Knight et al. (1998) and a new scale-similar model proposed here correlate well with the actual SGS turbulent diffusion, and predict the correct rms amplitude. However, the new scale-similar model was found to be more accurate. A mixed model for the SGS viscous diffusion was also proposed and tested, although this term is much smaller than the others. The accuracy of the models for the sum of the terms was also evaluated, and it was found that the models proposed still predict nearly the correct rms amplitude, and an acceptable value of the correlation coefficient.

The results obtained in this investigation are promising and indicate that it is possible to model accurately the terms in the energy equations. Further work may extend these results to cases in which the pressure-dilatation is significant, as well as to inhomogeneous flows, and evaluate these models a posteriori.

References

Bardina, J., Ferziger, J.H., and Reynolds, W.C. (1980). Improved subgrid-scale models for large eddy simulation. AIAA

Paper 80-1357. Comte, P., and Lesieur, M. (1998). Large-eddy simulation of compressible turbulent flows. In: Advances in Turbulence

Modeling, edited by D. Olivari. Von Karman Institute for Fluid Dynamics, Rhode-Ste-Genèse, 4:1–4:133. El-Hady, N., Zang, T.A. , and Piomelli, U. (1994). Application of the dynamic subgrid-scale model to axisymmetric

transitional boundary layer at high speed. Phys. Fluids A, 6, 1299–1309. Erlebacher, G., Hussaini, M.Y., Speziale, C.G., and Zang, T.A. (1992). Toward the large-eddy simulation of compressible

turbulent flows. J. Fluid Mech., 238, 155–185. Favre, A. (1965a). Équations des gaz turbulents compressible. I. Formes générales. J. Méc., 4, 361–390. Favre, A. (1965b). Équations des gaz turbulents compressible. II. Méthode des vitesses moyennes; méthode des vitesses

macroscopiques pondérées par la masse volumique. J. Méc., 4, 391–421. Germano, M. (1992). Turbulence: the filtering approach. J. Fluid Mech., 238, 325–336. Knight, D., Zhou, G., Okong’o, N., and Shukla, V. (1998). Compressible large eddy simulation using unstructured grids.

AIAA Paper 98-0535. Leonard, A. (1974). Energy cascade in large-eddy simulations of turbulent fluid flows. Adv. Geophys., 18A, 237–248. Lilly, D.K. (1992). A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids A, 4, 633–635. Martín, M.P., and Candler, G.V. (1998). Effect of chemical reactions on decaying isotropic turbulence. Phys. Fluids, 10,

1715–1724. Martín, M.P., and Candler, G.V. (1999). Subgrid-scale model for the temperature fluctuations in reacting hypersonic

turbulent flows. Phys. Fluids, 11, 2765–2771. Meneveau, C., and Lund, T.S. (1997). The dynamic Smagorinsky model and scale-dependent coefficients in the viscous

range of turbulence. Phys. Fluids, 9, 3932–3934. Moin, P., Squires, K.D., Cabot, W.H., and Lee, S. (1991). A dynamic subgrid-scale model for compressible turbulence

and scalar transport. Phys. Fluids A, 3, 2746–2757.

16  M.P. Martín, U. Piomelli, and G.V. Candler

Normand, X., and Lesieur, M. (1992). Direct and large-eddy simulation of laminar breakdown in high-speed axisymmetric

boundary layers. Theoret. Comput. Fluid Dynamics, 3, 231–252. Ristorcelli, J.R., and Blaisdell, G.A. (1997). Consistent initial conditions for the DNS of compressible turbulence. Phys. Flu-

ids, 9, 4–6. Smagorinsky, J. (1963). General circulation experiments with the primitive equations. I. The basic experiment. Mon.

Weather Rev. 91, 99–164. Speziale, C.G., Erlebacher, G., Zang, T.A., and Hussaini, M.Y. (1988) The subgrid-scale modeling of compressible tur-

bulence. Phys. Fluids A, 31, 940–942. Vreman, B., Geurts, B., and Kuerten, H. (1995a). A priori tests of large eddy simulation of the compressible mixing layer.

J. Engrg. Math., 29, 299–327. Vreman, B., Geurts, B., and Kuerten, H. (1995b). Subgrid-modeling in LES of compressible flow. Appl. Sci. Res., 54,

191–203. Yoshizawa, A. (1986). Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling.

Phys. Fluids A, 29, 2152–2164. Zang, T.A., Dahlburg, R.B., and Dahlburg, J.P. (1992). Direct and large-eddy simulations of three-dimensional compress-

ible Navier–Stokes turbulence. Phys. Fluids A, 4, 127–140.