© ESO, 2011

1. Introduction

The effects of numerically unresolved turbulence have recently met increasing attention in a variety of astrophysical simulations (see, for instance, Scannapieco & Brüggen 2008; Maier et al. 2009; Joung et al. 2009; Oppenheimer & Davé 2009). Some approaches comprise subgrid scale (SGS) models, although these models are basically phenomenological parameterizations of astrophysical processes on length scales smaller than the grid scale. The full multi-scale dynamics of turbulence, however, is not embraced. The essence of an SGS model in the fluid-dynamical sense is that, at high Reynolds numbers, energy is transported through a turbulent cascade from larger, numerically resolved length scales to the subgrid scales. The energy of the unresolved turbulent velocity fluctuations is eventually dissipated into heat. Numerical simulations, in which an explicit closure for the turbulent cascade is applied on the grid scale, are called large eddy simulations (LES). A closure is an approximation to an SGS quantity in terms of resolved flow quantities.

If the unresolved turbulent velocity fluctuations reach a non-negligible fraction of the speed of sound, they give rise to a turbulent pressure in addition to the thermal pressure of the gas. This contribution to the pressure is proportional to the energy density of the SGS turbulence. Turbulent pressure effects in the baryonic gas component of star-forming galaxies are discussed in Burkert et al. (2010). In contemporary numerical simulations of disk galaxies (e.g., Dobbs et al. 2008; Tasker & Tan 2009; Agertz et al. 2009, 2010) or in galactic-scale simulations of the interstellar medium (e.g., Joung & Mac Low 2006; Joung et al. 2009), the minimal grid scale (or the SPH smoothing length) Δ ≳ 1   pc. Since this length scale is comparable to the size of molecular clouds, the unresolved turbulent velocity fluctuations can exceed the speed of sound in the cold-gas phase. Consequently, it can be expected that a significant turbulent pressure is caused by turbulence below the grid scale. To compute the turbulent pressure, which has an impact on the star formation rate through the stability of the gas against gravitational collapse, a model for the highly compressible regime is indispensable. Bonazzola et al. (1987) and Bonazzola et al. (1992) formulate an analytic theory to calculate the turbulent pressure of isotropic compressible turbulence on the integral length scale. By applying an SGS model, on the other hand, the turbulent pressure can be computed for any length scale within the inertial subrange. Joung et al. (2009) present an SGS model that is based on the equation for the kinetic energy of the unresolved turbulent velocity fluctuations, the so-called SGS turbulence energy, where energy is solely supplied by supernova feedback. Since the non-diagonal SGS turbulence stresses are neglected, the model of Joung et al. (2009) reduces the effects of SGS turbulence to the turbulent pressure alone, and the turbulence energy cascade, i.e., the production of SGS turbulence by the shear of the numerically resolved flow, is not considered.

Another example for this type of SGS models is the model of Scannapieco & Brüggen (2008) for the simulation of Rayleigh-Taylor-driven turbulence in active galactic nuclei, where it is assumed that SGS turbulence is produced by buoyancy processes only on unresolved length scales. These processes are modelled by an equation for the characteristic length scale of the Rayleigh-Taylor instability. In contrast, Schmidt et al. (2006b) incorporate unresolved buoyancy effects into an SGS model that includes the production by shear for the treatment of turbulent combustion in thermonuclear supernovae.

In the cosmological simulations of galaxy clusters by Maier et al. (2009), the role of SGS turbulence has been explored with a numerical technique that combines adaptive mesh refinement and LES. They apply the SGS turbulence energy model of Schmidt et al. (2006a). The main effect of the SGS model is an enhancement of the turbulent heating in the cluster core. The SGS turbulence energy also serves as a tracer of turbulence production in the intergalactic medium (Iapichino et al. 2011). However, a deficiency in these simulations is that the employed SGS model is only applicable to moderately compressible turbulence. Shocks are treated tentatively, i.e., SGS turbulence production is suppressed in the vicinity of shock fronts. While this is not a severe constraint for the bulk of the intracluster medium, in which the Mach numbers of the turbulent flow are small compared to unity, an erroneous production of SGS turbulence energy is likely to occur near accretion shocks in the outer regions of the cluster.

In this article, we improve on the previous approaches to SGS modelling by addressing the closure problem for highly compressible turbulence. In Sect. 2, we discuss the meaning of the compressible Euler equations in the context of computational fluid dynamics. The verification of the proposed closure and the calibration of the closure coefficients are presented in Sect. 3. Data from several high-resolution simulations of forced turbulence (Schmidt et al. 2007, 2009; Federrath et al. 2010b) allow us to compute the rate at which energy is transferred from length scales greater than the filter length to smaller length scales and to test the correlation with different closures. As a result, we propose a combination of the eddy-viscosity closure, which has been used successfully in LES of incompressible turbulence, and a nonlinear closure put forward by Woodward et al. (2006). Then we show that physically reasonable statistics of the SGS turbulence energy and the rate of dissipation are obtained for varying grid resolutions and forcing in LES of supersonic turbulence (Sect. 4). Furthermore, we investigate correlations of the SGS quantities with quantities derived from the numerically resolved flow. We demonstrate that the turbulent pressure and energy dissipation cannot be predicted in a straightforward way on the basis of the resolved turbulent flow, as proposed, for instance, by Pan et al. (2009) and Zhu et al. (2010). Instead, a full SGS model is needed to estimate unresolved turbulence effects. In the last section, we summarize the results and discuss potential astrophysical applications of the closure for the highly compressible turbulent cascade in combination with the phenomenological approaches described above.

2. The compressible Euler equations with subgrid-scale dynamics

The Reynolds number of turbulent flows in astrophysics is usually considered to be high enough so that the approximation of an inviscid fluid can be applied on numerically resolvable length scales. For a physically complete picture, we begin with the compressible Navier-Stokes equations encompassing all physical length scales. This acknowledges that perfect fluids do not exist in nature and that the notion of viscous dissipation is essential for turbulence. The fluid-dynamical variables determined by this set of equations are denoted by $\stackrel{\mathrm{\infty }}{\mathit{\varrho }}$ for the mass density of the gas, $\stackrel{\mathrm{\infty }}{{u}}$ for the velocity of the flow, etc.¹. The resolution of a numerical simulation is given by the size of the grid cells Δ, which is called the cutoff scale or the grid scale. A consistent formulation of the equations of fluid dynamics with a cutoff scale Δ can by derived from the Navier-Stokes equations by means of the filter formalism introduced by Germano (1992). Generalizing this formalism to compressible fluid dynamics is straightforward (see Schmidt et al. 2006a). The basic idea is to identify the numerically computed solution with filtered variables $\mathit{\varrho }\mathrm{:}\mathrm{=}\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}{\mathrm{⟩}}_{\mathrm{\Delta }}$, ${u}\mathrm{:}\mathrm{=}\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathit{/}\mathit{\varrho }$, etc. The filter operator  ⟨ · ⟩ _Δ smooths the physical variables that are given by the Navier-Stokes equations over the length scale Δ. In LES, the filtering corresponds to the discretization of the equations of fluid dynamics. The dynamical equations for the computable, filtered quantities are similar to the unfiltered equations, with additional terms that are related to the subgrid-scale dynamics on length scales ℓ < Δ.

Let us consider the dynamical equation for the momentum density of the fluid, which is given by the partial differential equation (PDE) $\frac{\mathit{\partial }}{\mathit{\partial t}}\left(\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\right)\mathrm{+}\mathbf{\nabla }\mathrm{·}\left(\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\mathrm{\otimes }\stackrel{\mathrm{\infty }}{{u}}\right)\mathrm{=}\stackrel{\mathrm{\infty }}{\mathit{\rho }}\left(\stackrel{\mathrm{\infty }}{{g}}\mathrm{+}\stackrel{\mathrm{\infty }}{{f}}\right)\mathrm{-}\mathbf{\nabla }\stackrel{\mathrm{\infty }}{\mathit{P}}\mathrm{+}\mathbf{\nabla }\mathrm{·}\stackrel{\mathrm{\infty }}{{\sigma }}\mathit{,}$(1)where $\stackrel{\mathrm{\infty }}{\mathit{g}}$ and $\stackrel{\mathrm{\infty }}{\mathit{f}}$ are the accelerations due to gravity and other mechanical forces acting on the fluid, and $\stackrel{\mathrm{\infty }}{\mathit{P}}$ is the thermal pressure. The viscous dissipation tensor $\stackrel{\mathrm{\infty }}{\mathit{\sigma }}$ is defined by ${\stackrel{\mathrm{\infty }}{\mathit{\sigma }}}_{\mathit{ij}}\mathrm{=}\mathrm{2}\mathit{\nu }\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\left({\stackrel{\mathrm{\infty }}{\mathit{S}}}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\stackrel{\mathrm{\infty }}{\mathit{d}}{\mathit{\delta }}_{\mathit{ij}}\right)\mathit{,}$(2)where ν is the microscopic viscosity of the fluid², the rate of strain ${\stackrel{\mathrm{\infty }}{\mathit{S}}}_{\mathit{ij}}$ is the symmetic part of the Jacobian matrix ${\stackrel{\mathrm{\infty }}{\mathit{u}}}_{\mathit{i,j}}\mathrm{=}{\mathit{\partial }}_{\mathit{j}}{\stackrel{\mathrm{\infty }}{\mathit{u}}}_{\mathit{i}}$, and $\stackrel{\mathrm{\infty }}{\mathit{d}}\mathrm{=}{\mathit{u}}_{\mathit{i,i}}$. After applying a homogeneous filter operator that is uniform in time, Eq. (1) is converted into an equation for the filtered momentum, $\mathit{\varrho }{u}\mathrm{=}\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\mathrm{⟩}$. This equation has the same form as the original equation, except for one term. Because of the nonlinear advection term, the filtering introduces a stress term that accounts for the interaction between the numerically resolved flow and velocity fluctuations on the subgrid scales: $\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathit{\varrho }{u}\mathrm{\right)}\mathrm{+}\mathbf{\nabla }\mathrm{·}\mathrm{(}\mathit{\varrho }{u}\mathrm{\otimes }{u}\mathrm{)}\mathrm{=}\mathit{\rho }\mathrm{(}{g}\mathrm{+}{f}\mathrm{)}\mathrm{-}\mathbf{\nabla }\mathit{P}\mathrm{+}\mathbf{\nabla }\mathrm{·}\mathrm{(}{\sigma }\mathrm{+}{{\tau }}_{\mathrm{sgs}}\mathrm{)}\mathit{,}$(3)where the SGS turbulence stress tensor is defined by ${{\tau }}_{\mathrm{sgs}}\mathrm{=}\mathrm{-}\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\mathrm{\otimes }\stackrel{\mathrm{\infty }}{{u}}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathrm{+}\mathit{\varrho }{u}\mathrm{\otimes }{u}\mathit{.}$(4)In the following, the components of τ_sgs are simply denoted by τ_ij. The second-order moment $\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\mathrm{\otimes }\stackrel{\mathrm{\infty }}{{u}}{\mathrm{⟩}}_{\mathrm{\Delta }}$ is not explicitly computable in LES because the variations in the mass density $\stackrel{\mathrm{\infty }}{\mathit{\varrho }}$ and the velocity $\stackrel{\mathrm{\infty }}{{u}}$ below the grid scale are unknown. For this reason, an approximation in terms of filtered quantities has to be devised. This is the closure problem³.

The SGS turbulence energy density is defined by the difference between the resolved kinetic energy and the filtered kinetic energy: ${\mathit{K}}_{\mathrm{sgs}}\mathrm{:}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\stackrel{\mathrm{\infty }}{{u}}\mathrm{·}\stackrel{\mathrm{\infty }}{{u}}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\varrho }\mathrm{|}{u}{\mathrm{|}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{tr}\hspace{0.17em}{{\tau }}_{\mathrm{sgs}}\mathit{,}$(5)where tr   τ_sgs = τ_ii is the trace of the SGS turbulence stress tensor. One can see that the trace of τ_sgs gives rise to the term $\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathbf{\nabla }{\mathit{K}}_{\mathrm{sgs}}$ on the right hand side of Eq. (3). This term can be absorbed into the pressure gradient if the thermal pressure P is replaced by the effective pressure ${\mathit{P}}_{\mathrm{eff}}\mathrm{=}\mathit{P}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{\mathrm{sgs}}\mathrm{=}\mathit{P}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{tr}\hspace{0.17em}{{\tau }}_{\mathrm{sgs}}\mathit{.}$(6)The relative contribution of the turbulent pressure ${\mathit{P}}_{\mathrm{sgs}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{\mathrm{sgs}}$ compared to the thermal pressure P is characterized by the SGS turbulence Mach number ${\mathrm{ℳ}}_{\mathrm{sgs}}\mathrm{=}\mathrm{(}\mathrm{2}{\mathit{K}}_{\mathrm{sgs}}\mathit{/}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, where c_s is the thermal speed of sound. ℳ_sgs depends on the temperature of the fluid and the cutoff scale Δ. The dependence on Δ is investigated in Sect. 4.2. Joung et al. (2009) define the turbulent pressure by P_sgs = (γ − 1)K_sgs, where γ is the adiabatic coefficient of the gas. We emphasize that, except for γ = 5/3, this definition is inconsistent with the decomposition of the fluid-dynamical equations, which fixes the coefficient to 2/3 (see also Chandrasekhar 1951). This is reasonable because the turbulent pressure is solely a property of the turbulent flow of a gas on a given length scale, whereas γ is a microscopic property of the gas that is related to the thermal motions of the atoms or molecules.

The SGS turbulence energy is an intermediate reservoir of energy that exchanges energy with the resolved flow and loses energy by dissipation into heat. For the computation of K_sgs, a PDE has to be solved in addition to the filtered equations for the resolved gas dynamics: $\frac{\mathit{\partial }}{\mathit{\partial t}}{\mathit{K}}_{\mathrm{sgs}}\mathrm{+}\mathbf{\nabla }\mathrm{·}\mathrm{\left(}{u}{\mathit{K}}_{\mathrm{sgs}}\mathrm{\right)}\mathrm{=}\mathrm{\Gamma }\mathrm{+}\mathrm{\Sigma }\mathrm{-}\mathit{\rho }\mathrm{(}\mathit{ϵ}\mathrm{+}\mathit{\lambda }\mathrm{)}\mathrm{+}D\mathit{.}$(7)While Σ = τ_ijS_ij is the rate of SGS turbulence energy production by the turbulent cascade through the cutoff scale Δ (also called the turbulence energy flux) and ρϵ is the viscous dissipation rate smoothed over Δ, effects caused by SGS fluctuations of the gravitational potential and the thermal pressure are given by Γ and ρλ, respectively. The term D accounts for SGS transport effects. We refer to Schmidt et al. (2006a), Eqs. (33)–(37), for the exact definitions of these terms. For our purpose it is sufficient to discuss the closures of these terms, which are approximations in terms of the numerically resolved variables and K_sgs.

To compute the SGS turbulence stress tensor (4), we propose the following closure for the highly compressible regime: ${\mathit{\tau }}_{\mathit{ij}}\mathrm{=}\mathrm{2}{\mathit{C}}_{\mathrm{1}}\mathrm{\Delta }\mathrm{(}\mathrm{2}\mathit{\varrho }{\mathit{K}}_{\mathrm{sgs}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{S}}_{\mathit{ij}}^{\mathrm{\ast }}\mathrm{-}\mathrm{2}{\mathit{C}}_{\mathrm{2}}{\mathit{K}}_{\mathrm{sgs}}\frac{\mathrm{2}{\mathit{u}}_{\mathit{i,k}}{\mathit{u}}_{\mathit{j,k}}}{\mathrm{|}\mathbf{\nabla }\mathrm{\otimes }{u}{\mathrm{|}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{C}}_{\mathrm{2}}\mathrm{)}{\mathit{K}}_{\mathrm{sgs}}{\mathit{\delta }}_{\mathit{ij}}\mathit{,}$(8)where |∇ ⊗ u|: = (2u_i,ku_i,k)^1/2 is the norm of the resolved velocity derivative, ${\mathit{S}}_{\mathit{ij}}^{\mathrm{\ast }}\mathrm{=}{\mathit{S}}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}\mathit{d}{\mathit{\delta }}_{\mathit{ij}}$ is the trace-free part of ${\mathit{S}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{u}}_{\mathit{i,j}}\mathrm{+}{\mathit{u}}_{\mathit{j,i}}\mathrm{)}$, and d = u_i,i. While the first term in Eq. (8) corresponds to the eddy-viscosity closure that is commonly used in incompressible LES, the second, nonlinear term was investigated by Woodward et al. (2006) for transonic decaying turbulence. The standard eddy-viscosity closure follows if C₂ = 0. In general, the linear eddy-viscosity term dominates if (K_sgs/ρ)^1/2 is small compared to Δ|S^∗| ≲ Δ|∇ ⊗ u|. On the other hand, for strong turbulence intensity, i.e., (K_sgs/ρ)^1/2 ≳ Δ|∇ ⊗ u|, the nonlinear term contributes significantly. This particularly applies to intermittent events in supersonic turbulence, for which Δ|∇ ⊗ u| ≳ c_s. In moderately compressible turbulence, nonlinear contributions affect the high-intermittency tails of the turbulent energy distribution. Independent of the values of C₁ and C₂, τ_ii = −2K_sgs, as required by the identity (5). We denote the trace-free part of the SGS turbulence stress tensor by ${\mathit{\tau }}_{\mathit{ij}}^{\mathrm{\ast }}$. The verification of the generalized closure (8) for τ_ij and the determination of the coefficients C₁ and C₂ for supersonic turbulence is the key to computing the turbulent pressure ${\mathit{P}}_{\mathrm{sgs}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{\mathrm{sgs}}$, as K_sgs first and foremost depends on the production rate Σ = τ_ijS_ij in Eq. (7).

Due to the microscopic viscosity ν of the fluid, the viscous stresses ${\stackrel{\mathrm{\infty }}{\mathit{\sigma }}}_{\mathit{ij}}$ dissipate kinetic energy on the smallest dynamical length scales ℓ ~ η of the physical flow $\stackrel{\mathrm{\infty }}{{u}}$. The length scale η is called the Kolmogorov scale. In the filtered momentum Eq. (3), viscous dissipation effects are given by the divergence of the filtered tensor ${\mathit{\sigma }}_{\mathit{ij}}\mathrm{=}\mathrm{⟨}{\stackrel{\mathrm{\infty }}{\mathit{\sigma }}}_{\mathit{ij}}{\mathrm{⟩}}_{\mathrm{\Delta }}$. The corresponding rate of energy dissipation, filtered on the grid scale, is given by $\mathit{\varrho ϵ}\mathrm{=}\mathrm{⟨}{\stackrel{\mathrm{\infty }}{\mathit{\sigma }}}_{\mathit{ij}}{\stackrel{\mathrm{\infty }}{\mathit{u}}}_{\mathit{i,j}}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathrm{=}\mathrm{⟨}\mathrm{2}\mathit{\nu }\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\hspace{0.17em}\stackrel{\mathrm{\infty }}{\mathit{S}}\begin{array}{}\hspace{0.17em}\mathrm{\ast }\\ \mathit{ij}\end{array}\stackrel{\mathrm{\infty }}{\mathit{S}}\begin{array}{}\hspace{0.17em}\mathrm{\ast }\\ \mathit{ij}\end{array}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathrm{=}\mathrm{⟨}\mathit{\nu }\stackrel{\mathrm{\infty }}{\mathit{\varrho }}\mathrm{|}\stackrel{\mathrm{\infty }}{\mathit{S}}\begin{array}{}\hspace{0.17em}\mathrm{\ast }\\ \mathit{ij}\end{array}{\mathrm{|}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathrm{\Delta }}\mathit{.}$(9)It is important to note that ϱϵ ≠ σ_iju_i,j, where σ_ij and u_i,j are the filtered viscous stress tensor and the filtered velocity gradient, respectively.

For fully developed incompressible turbulence, the Kolmogorov scale can be related to the Reynolds number: η/L ~ Re^3/4, where Re: = VL/ν for an integral length L and characteristic velocity V of the flow. As pointed out at the beginning of this section, Re is assumed to be very high in astrophysical systems. In this case, η is much smaller than any feasible grid resolution Δ, and simple scaling arguments show that the viscous stress term in the filtered momentum Eq. (3) is negligible (Röpke & Schmidt 2009), i.e., |σ| ≪ |τ_sgs|. Consequently, the physical energy dissipation occurs entirely on subgrid scales ℓ ≪ Δ. As η decreases in comparison to Δ, the velocity fluctuations on ever smaller length scales give rise to arbitrarily steep velocity gradients, which add up to a non-vanishing product of the viscosity times the squared rate of strain on the right hand side of Eq. (9), regardless of how small the viscosity is. This results in a non-zero, asymptotically constant mean rate of energy dissipation in the limit η → 0 (ν → 0), which is supported by experimental and numerical evidences (see Frisch 1995; Ishihara et al. 2009). We may reasonably conjecture that the viscous dissipation tensor is negligible in the filtered momentum equation and the energy dissipation rate does not vanish in the limit of infinite Reynolds numbers also in the case of compressible turbulence. A posteriori tests imply that this conjecture is fulfiled for driven supersonic turbulence (see Sect. 4.2). However, the question of energy dissipation in inhomogeneous turbulence is more difficult. For example, it is known from boundary layers of terrestrial turbulence that viscous effects can affect the flow on relatively large scales near a wall. For this reason, the microscopic viscosity cannot be neglected in LES of such flows. Although solid walls are not encountered in astrophysics, many relevant problems exhibit pronounced inhomogeneities, and we cannot entirely exclude the possibility that viscous effects might become noticeable on resolved length scales in certain cases.

A closure for ϵ follows from simple dimensional reasoning: $\mathit{\rho ϵ}\mathrm{=}{\mathit{C}}_{\mathit{ϵ}}\frac{{\mathit{K}}_{\mathrm{sgs}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{\rho }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\Delta }}\mathrm{·}$(10)Here, it is assumed that the SGS turbulence energy is dissipated into heat at a rate proportional to K_sgs divided by the time scale Δ(K_sgs/ρ)^−1/2. For the pressure-dilatation term ρλ several closures have been proposed (e.g., Sarkar 1992; Fureby et al. 1997). However, applying a priori tests (see Sect. 3), we find that these closures clearly fail in the case of supersonic turbulence. The simplest solution is to neglect pressure dilatation entirely (Woodward et al. 2006). In this article, we also set ρλ = 0, although we are aware that pressure-dilatation effects have potential significance, particularly in the case of adiabatic turbulence. The transport term in Eq. (7) can be modelled by a gradient-diffusion approximation (see Sagaut 2006): $D\mathrm{=}\mathbf{\nabla }\mathrm{·}\left[{\mathit{\kappa }}_{\mathrm{sgs}}\mathbf{\nabla }\left(\frac{{\mathit{K}}_{\mathrm{sgs}}}{\mathit{\rho }}\right)\right]\mathit{,}$(11)where the SGS turbulent diffusivity is approximated by κ_sgs ≈ 0.65Δ(ρK_sgs)^1/2, as shown by Schmidt et al. (2006a).

In this work, we assume that self-gravity has no significant effects on length scales ℓ ≲ Δ. This corresponds to the condition that the local Jeans length λ = c_s(π/Gρ)^1/2, where G is the gravitational constant, is sufficiently large compared to the grid scale Δ (Truelove et al. 1997; Federrath et al. 2010a). Thus, setting Γ = 0, the filtered equations resulting from the compressible Navier-Stokes equations in the limit of η ≪ Δ (Schmidt et al. 2006a) read as $\begin{array}{cc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{\varrho }\mathrm{+}\mathbf{\nabla }\mathrm{·}\mathrm{\left(}{u}\mathit{\varrho }\mathrm{\right)}\hspace{0.17em}& \mathrm{=}\mathrm{0}\mathit{,}\\ \mathrm{=}\mathit{\rho }\mathrm{(}{g}\mathrm{+}{f}\mathrm{)}\mathrm{-}\mathbf{\nabla }\mathrm{(}\mathit{P}\mathrm{+}{\mathit{P}}_{\mathrm{sgs}}\mathrm{)}\mathrm{+}\mathbf{\nabla }\mathrm{·}{{\tau }}_{\mathrm{sgs}}^{\mathrm{\ast }}\mathit{,}& \end{array}$where $\mathit{E}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{u}}^{\mathrm{2}}\mathrm{+}{\mathit{E}}_{\mathrm{int}}$ is the sum of the resolved kinetic and internal energy density, and −ℒ accounts for sources and sinks of the internal energy due to heating and cooling, respectively. Since the resolved fluid dynamics on length scales ℓ ≥ Δ is unaffected by the viscosity of the fluid, the above set of equations defines the compressible Euler equations for computational fluid dynamics in a physically meaningful and consistent way. These equations are supplemented by an equation of state, the SGS turbulence energy Eq. (7), and the Poisson equation for the gravitational potential. The pure compressible Euler equations without SGS terms, on the other hand, do not follow from the compressible Navier-Stokes equation in the limit of infinite Reynolds number. In this case, there is no viscous dissipation at all, and, by definition, ϵ vanishes identically. This is a mathematical idealization that does not describe turbulent flows in nature.

As a special case, implicit large eddy simulations (ILES) follow from the above approach. ILES is the most commonly used method in astrophysical fluid dynamics. It is based on two assumptions, which are usually not stated in the literature. First, the discretization of the compressible Euler equations introduces a dissipative leading error term D_num in the momentum Eq. (11). Implicitly, this term is assumed to be equivalent to the SGS turbulence stress term ∇·τ_sgs. The second assumption in ILES is that u·D_num = −ρϵ, i.e., kinetic energy on the resolved scales is directly dissipated into heat at a rate that approximates the viscous dissipation on unresolved length scales. This is referred to as numerical viscosity or numerical dissipation. Actually, the following equations are solved in ILES: $\begin{array}{cc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{\varrho }\mathrm{+}\mathbf{\nabla }\mathrm{·}\mathrm{\left(}{u}\mathit{\varrho }\mathrm{\right)}\hspace{0.17em}& \mathrm{=}\mathrm{0}\mathit{,}\\ \mathrm{=}\mathit{\rho }\mathrm{(}{g}\mathrm{+}{f}\mathrm{)}\mathrm{-}\mathbf{\nabla }\mathit{P}\mathrm{+}{D}_{\mathrm{num}}\mathit{,}& \end{array}$Despite the lack of a mathematical justification, ILES serves as an approximation of turbulent compressible fluid dynamics that has proven its utility in numerous astrophysical applications. Benzi et al. (2008) demonstrate that the ILES approach closely reproduces two-point statistics of weakly compressible turbulence in the inertial subrange in comparison to direct numerical simulations that solve the Navier-Stokes equations. In this article, we make use of ILES to compute high-resolution data for the explicit verification of SGS closures. In contrast, LES treat the energy dissipation explicitly. However, it cannot be avoided that numerical schemes for compressible fluid dynamics such as the piecewise parabolic method (PPM, Colella & Woodward 1984) introduce some numerical dissipation. Thus, by running an LES with an explicit SGS model, there is inevitably a numerical dissipation channel that competes with the transfer of energy to the subgrid scales and the subsequent dissipation of SGS turbulence energy into heat. Notwithstanding this caveat, we demonstrate in this article that physically sensible predictions can by made by an explicit SGS model, which are not possible on the basis of ILES.

3. Closure verification

To test different closures for the turbulence energy flux Σ, we apply the method described in Schmidt et al. (2006a). The basic idea is to use data from ILES of non-selfgravitating isothermal and adiabatic turbulence with high numerical resolution and to apply an explicit filter to these data on a length scale that is in between the forcing and the dissipative range. The applied filters are Gaussian with filter lengths Δ_G = 24Δ and 32Δ for 768³ and 1024³ grids, respectively. Although the bottleneck effect has some influence on the chosen length scales (Schmidt et al. 2009; Federrath et al. 2010b), we show that the results change only slightly if the filter length decreases or increases by a factor of two. Moreover, comparing to box filters, the results turn out to be rather insensitive to the filter type.

3.1. Single-coefficient closures for supersonic isothermal turbulence

The turbulence energy flux on the filter scale Δ_G can be computed explicitly from the unfiltered numerical data by the formula ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{=}\left(\mathrm{-}\overline{\mathit{\varrho }{\mathit{u}}_{\mathit{i}}{\mathit{u}}_{\mathit{j}}}\mathrm{+}\frac{\overline{\mathit{\varrho }{\mathit{u}}_{\mathit{i}}}\hspace{0.17em}\overline{\mathit{\varrho }{\mathit{u}}_{\mathit{i}}}}{\overline{\mathit{\varrho }}}\right){\overline{\mathit{u}}}_{\mathit{i,j}}\hspace{0.17em}\mathit{,}$(18)where the first factor on the right hand side is the turbulence stress tensor on the filter scale (defined analogous to Eq. (4)) and the second factor is the derivative of the filtered velocity. As a shorthand notation, we denote the explicitly filtered quantities by an overline, for instance, $\overline{\mathit{\rho }}\mathrm{=}\mathrm{⟨}\mathit{\varrho }{\mathrm{⟩}}_{{\mathrm{\Delta }}_{\mathrm{G}}}$, and ${\overline{\mathit{u}}}_{\mathit{i}}\mathrm{=}\mathrm{⟨}\mathit{\varrho }{\mathit{u}}_{\mathit{i}}{\mathrm{⟩}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathit{/}\overline{\mathit{\rho }}$.

Fig. 1 Correlation diagrams for the SGS turbulence energy flux in the case of isothermal supersonic turbulence with solenoidal forcing (ℳrms ≈ 5.3). The applied filter length is 32Δ. The blue dots indicate the average prediction of the closure for a given value of Σ32Δ.

The above expression for Σ_ΔG can be compared to closures. For the eddy-viscosity closure, Σ_ΔG is given by ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{=}{\mathit{C}}_{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{G}}\mathrm{\left(}\mathrm{2}\overline{\mathit{\rho }}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right|}{\overline{\mathit{S}}}^{\hspace{0.17em}\mathrm{\ast }}{\mathrm{|}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}\mathit{,}$(19)where ${\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\overline{\mathit{\varrho }{\mathit{u}}^{\mathrm{2}}}\mathrm{-}\frac{{\overline{\mathit{\varrho u}}}^{\mathrm{2}}}{\overline{\mathit{\varrho }}}\right)$(20)is the turbulence energy on length scales smaller than Δ_G. Strictly speaking, K_ΔG is the turbulence energy for the length scales ranging from the grid resolution Δ to the smoothing length Δ_G of the Gaussian filter. If Δ_G is sufficiently larger than Δ, this distinction can be neglected (see Schmidt et al. 2006a).

When defining ${\mathit{C}}_{\mathrm{1}}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{:}\mathrm{=}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}\mathit{,}$(21)the squared error function of the closure can be written as ${\mathrm{err}}^{\mathrm{2}}\mathrm{\left(}{\mathit{C}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{𝒱}}{\left|{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}\mathrm{-}{\mathit{C}}_{\mathrm{1}}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\right|}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{3}}\mathit{x,}$(22)where Σ_ΔG and K_ΔG are given by Eqs. (18) and (20), respectively, and the volume integral extends over the whole domain . The minimum of err²(C₁) yields the least squares error solution for the closure coefficient, ${\mathit{C}}_{\mathrm{1}}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{𝒱}}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\overline{\mathit{d}}}^{\mathrm{]}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{x}}{{\mathrm{\int }}_{\mathrm{𝒱}}\mathrm{|}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathrm{|}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{3}}\mathit{x}}\mathit{,}$(23)where ${\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{=}{\mathrm{\Delta }}_{\mathrm{G}}\mathrm{\left(}\mathrm{2}\overline{\mathit{\rho }}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right|}{\overline{\mathit{S}}}^{\hspace{0.17em}\mathrm{\ast }}{\mathrm{|}}^{\mathrm{2}}$ for the eddy-viscosity closure. For statistically stationary and isotropic turbulence, the closure coefficient C₁ is independent of the filter length scale because of the local equilibrium of the transfer of turbulence energy in the inertial subrange. Thus, the value of C₁ inferred from Eq. (23) is an approximation to the coefficient of the SGS closure for Σ in LES.

To calculate C₁, we use data from two 1024³ simulations of supersonic isothermal turbulence with a root-mean-square (RMS) Mach number around 5.5 (Federrath et al. 2010b). Statistically stationary and isotropic turbulence is produced by stochastic forcing. Solenoidal (divergence-free) forcing is applied in one simulation, while the forcing is compressive (rotation-free) in the other simulation. We choose Δ_G = 32Δ for the filtering of the simulation data. Figures 1 and 2 show the correlation between ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}$ and Σ_ΔG by means of two-dimensional probability density functions. For the eddy-viscosity closure (19), the correlation is quite good (the spacing of the contour lines in Figs. 1 and 2 is logarithmic in the two-dimensional probability density), but there is a problem with negative flux values. The values of the closure coefficient C₁ following from Eq. (23) are listed in Table 1. Also listed are the correlation coefficients $\begin{array}{ccc}\mathrm{corr}& \mathrm{\left[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathit{,}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{\right]}\mathrm{=}& \\ & \frac{{\mathrm{\int }}_{\mathrm{𝒱}}\left[{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{-}\mathrm{⟨}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{⟩}\right]\left[{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{-}\mathrm{⟨}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{⟩}\right]\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{x}}{\mathrm{std}\mathrm{\left[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{\right]}\mathrm{std}\mathrm{\left[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{\right]}}\mathit{,}& \end{array}$(24)where std [·]  denotes the standard deviation and the angle brackets indicate an average over the whole domain.

Fig. 2 Correlation diagrams for isothermal supersonic turbulence with compressive forcing (ℳrms ≈ 5.6) as in Fig. 1.

An important question for the LES of supersonic turbulence is whether shocks can be accommodated in closures for the turbulence energy. Since the eddy-viscosity closure originates from incompressible turbulence, Maier et al. (2009) suggests setting Σ equal to zero in the vicinity of shock fronts. This should suppress the spurious production of SGS turbulence energy by the large strain at shock fronts. Thus, we tested whether excluding shocks in the computation of the eddy-viscosity closure for the turbulence energy flux would improve the correlations. However, panels (b) in Figs. 1 and 2 make it clear that such a cutoff deteriorates the correlations and implies a significant underestimate of large positive fluxes. Although the applied shock detection criterion $\overline{\mathit{d}}\mathit{<}\mathrm{-}{\overline{\mathit{c}}}_{\mathrm{s}}\mathit{/}{\mathrm{\Delta }}_{\mathrm{G}}$ is rather crude, we interpret this trend as an indication that shocks must not be separated from the supersonic turbulent cascade.

Fig. 3 Slices of the two-dimensional probability density functions plotted in Figs. 1 and 2, showing the predictions of different closures for the values of the explicitly computed energy flux Σ32Δ that are indicated by the vertical dashed lines. The top and bottom rows of panels show the results for solenoidal and compressive forcing, respectively.

In addition to the conventional eddy-viscosity closure, we investigate a closure that is based on the determinant of the velocity gradient (Woodward et al. 2001). In this case, the trace-free part of the SGS turbulence stress tensor is still given by the expression ${\mathit{\tau }}_{\mathit{ij}}^{\mathrm{\ast }}\mathrm{=}\mathrm{2}\mathit{\rho }{\mathit{\nu }}_{\mathrm{sgs}}{\mathit{S}}_{\mathit{ij}}^{\mathrm{\ast }}$. The eddy viscosity, however, does not depend on K_sgs, but is defined by ν_sgs = −C₁Δ²|S^∗|⁻²detS^∗. Hence, the turbulence energy flux on the filter scale is given by ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathit{C}}_{\mathrm{1}}\overline{\mathit{\rho }}{\mathrm{\Delta }}_{\mathrm{G}}^{\mathrm{2}}\begin{array}{}\det \end{array}{\overline{\mathrm{S}}}^{\hspace{0.17em}\mathrm{\ast }}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}$(25)for this closure. Woodward et al. (2001) employ the same method as we do to test the correlation between their closure and the turbulence energy flux. A particularly interesting feature of the determinant is that it switches signs, thereby describing two different flow topologies. In one case, the determinant is negative. This corresponds to the forward turbulent cascade transporting energy from large to smaller eddies. In the other case, the flow is contracting in one dimension and expanding in the other two. Then the determinant is positive, corresponding to a backscattering of energy from small to larger eddies. This phenomenon can be explained by the alignment of vortices along a single stretching direction (the “tornado" topology). While an energy flux of the form (19) fails to describe the reverse cascade, we see in panels (c) of Figs. 1 and 2 that the determinant closure yields a good correlation for negative energy flux. However, the overall correlation does not significantly improve (see Table 1), because of the relatively large scatter in the forward cascade.

In Woodward et al. (2006), a nonlinear expression for the turbulence stress tensor is investigated, which depends on the full Jacobian ∇ ⊗ u of the velocity: ${\mathit{\tau }}_{\mathit{ij}}\mathrm{=}\mathrm{-}\mathrm{2}{\mathit{C}}_{\mathrm{1}}{\mathit{K}}_{\mathrm{sgs}}\frac{\mathrm{2}{\mathit{u}}_{\mathit{i,k}}{\mathit{u}}_{\mathit{j,k}}}{\mathrm{|}\mathbf{\nabla }\mathrm{\otimes }{u}{\mathrm{|}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{C}}_{\mathrm{1}}\mathrm{)}{\mathit{K}}_{\mathrm{sgs}}{\mathit{\delta }}_{\mathit{ij}}\mathit{.}$(26)Since |∇ ⊗ u| = (2u_i,ku_i,k)^1/2, the above expression fulfills the identity τ_ii = −2K_sgs. The corresponding turbulence energy flux on the filter length scale Δ_G is given by ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{=}\mathrm{-}\mathrm{4}{\mathit{C}}_{\mathrm{1}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\frac{{\overline{\mathit{u}}}_{\mathit{i,k}}{\overline{\mathit{u}}}_{\mathit{j,k}}\overline{\mathit{S}}\begin{array}{}\hspace{0.17em}\mathrm{\ast }\\ \mathit{ij}\end{array}}{\mathrm{|}\mathbf{\nabla }\mathrm{\otimes }\overline{{u}}{\mathrm{|}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{C}}_{\mathrm{1}}\mathrm{)}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}\mathit{.}$(27)Figures 1d and 2d show that the correlation is excellent for the above closure, with correlation coefficients above 0.99, as listed in Table 1. Like the determinant closure discussed above, the trace-free part of the nonlinear closure for the SGS turbulence stress switches signs and, thus, allows for a backward energy cascade. To compare the quality of the different closures, slices of the two-dimensional probability density functions are plotted in Fig. 3. While both the determinant and the nonlinear closures yield good approximations to negative values of the energy flux, the nonlinear closure is cleary superior for large positive fluxes. However, the tails are flatter for compressively driven turbulence, which implies that large deviations are more frequent in this case. The eddy-viscosity closure closely reproduces low positive values of the energy flux. For high positive values, the deviations are less pronounced in comparison to the determinant closure, but greater than for the nonlinear closure.

3.2. Mixed closure for supersonic isothermal turbulence

Even though the correlation of the turbulence energy flux is very good, the purely nonlinear closure (26) is generally not adequate as a model for the turbulence stress tensor for the following reasons. Most importantly, rotation invariance is violated because of the antisymmetric part of ∇ ⊗ u. As a consequence, spurious turbulence energy would be produced for a uniformly rotating fluid. Apart from that, the application of this closure in LES of forced turbulence show that the growth of turbulence energy during the transition from laminar to turbulent flow is insufficient for this closure. This is because of the linear dependence on K_sgs⁴. As pointed out by Woodward et al. (2006), a seed term has to be included in order to trigger the production of turbulence energy. If the seed term is constructed from the symmetric part of the velocity gradient, then turbulence energy production vanishes for a uniformly rotating fluid, and, consequently, the problem of rotation invariance is also resolved. Woodward et al. (2006) consider a linear combination of the closure (26) with the determinant closure. Because of the relatively large scatter of the determinant closure for large positive energy fluxes, however, we propose a combination of the nonlinear closure with the linear eddy-viscosity closure. Conceptually, this combination has the advantage that the eddy-viscosity closure, which is well established for LES of incompressible turbulence, follows as a limiting case.

The least-squares-error approach can be generalized to a mixed closure with two coefficients, C₁ and C₂. For ${\mathit{C}}_{\mathrm{1}}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{=}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\overline{\mathit{d}}\mathit{,}$(28)the closure coefficients are given by the linear system of equations $\begin{array}{ccc}& {\mathrm{(}}^{\mathrm{\int }}\mathrm{𝒱}\mathrm{|}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathrm{|}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{)}}\hspace{0.17em}{\mathit{C}}_{\mathrm{1}}\mathrm{+}{\mathrm{(}}^{\mathrm{\int }}\mathrm{𝒱}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{)}}\hspace{0.17em}{\mathit{C}}_{\mathrm{2}}& \\ & \mathrm{=}{\mathrm{\int }}_{\mathrm{𝒱}}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\overline{\mathit{d}}}^{\mathrm{]}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{x,}& \\ & {\mathrm{(}}^{\mathrm{\int }}\mathrm{𝒱}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{)}}\hspace{0.17em}{\mathit{C}}_{\mathrm{1}}\mathrm{+}{\mathrm{(}}^{\mathrm{\int }}\mathrm{𝒱}\mathrm{|}{\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}{\mathrm{|}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{)}}\hspace{0.17em}{\mathit{C}}_{\mathrm{2}}& \\ & \mathrm{=}{\mathrm{\int }}_{\mathrm{𝒱}}{\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{[}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\overline{\mathit{d}}}^{\mathrm{]}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{x,}& \end{array}$where $\begin{array}{ccc}{\mathit{f}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}& \mathrm{=}{\mathrm{\Delta }}_{\mathrm{G}}\mathrm{\left(}\mathrm{2}\overline{\mathit{\rho }}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\right|}{\overline{\mathit{S}}}^{\hspace{0.17em}\mathrm{\ast }}{\mathrm{|}}^{\mathrm{2}}\mathit{,}& \\ {\mathit{g}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}& \mathrm{=}\mathrm{-}\mathrm{4}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\frac{{\overline{\mathit{u}}}_{\mathit{i,k}}{\overline{\mathit{u}}}_{\mathit{j,k}}\overline{\mathit{S}}\begin{array}{}\hspace{0.17em}\mathrm{\ast }\\ \mathit{ij}\end{array}}{\mathrm{|}\mathbf{\nabla }\mathrm{\otimes }\overline{{u}}{\mathrm{|}}^{\mathrm{2}}}\mathrm{·}& \end{array}$The solutions for C₁ and C₂ that are obtained from our numerical data are listed in Table 2. As one can see, the correlation coefficients are about as high as for the purely nonlinear closure and there is only little variation with the forcing and the filter length scale. For Δ_G = 64Δ the ratio of Δ_G to the integral scale L is 8, which is quite small. As a consequence, the driving force marginally influences this length scale. The filter length Δ_G = 16Δ, on the other hand, is significantly affected by numerical dissipation. The correlation diagrams for the mixed closure are plotted in Fig. 4. Although the relation between ${\mathrm{\Sigma }}_{\mathrm{32}\mathrm{\Delta }}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}$ and Σ_32Δ is slightly tilted for negative fluxes, the results are comparable to Figs. 1d and 2d. This can also be seen in Fig. 3. Therefore, we base our SGS model on the mixed nonlinear closure (8) with the averaged coefficients C₁ = 0.02 and C₂ = 0.7.

Fig. 4 Correlation diagrams for the mixed closures for isothermal supersonic turbulence with different forcing.

Fig. 5 Correlation diagrams for the mixed closure in the case of isothermal a) and adiabatic b) turbulence with lower rms Mach numbers as in Fig. 4.

3.3. Supplementary tests for different Mach numbers

For the data listed in Table 2, the average Mach numbers associated with the filter scale, ${\mathrm{ℳ}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{=}\mathrm{⟨}\mathrm{(}\mathrm{2}{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathit{/}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{⟩}$, assume values around the speed of sound. Thus, the question arises whether the closure coefficients calculated above are applicable to subsonic velocity fluctuations. To test the mixed closure for a different Mach number, we calculated the turbulence energy flux ${\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}$ with fixed values C₁ = 0.02 and C₂ = 0.7 for data from a simulation of isothermal turbulence with ℳ_rms ≈ 2.2 (Schmidt et al. 2009). The resulting correlation diagram is plotted in Fig. 5a. There is a small bias to overestimate the turbulence energy flux, but the prediction of the mixed closure is still very good. Indeed, a correlation coefficient $\mathrm{corr}\mathrm{\left[}{\mathrm{\Sigma }}_{\mathrm{G}}\mathit{,}{\mathrm{\Sigma }}_{\mathrm{G}}^{\mathrm{\left(}\mathrm{cls}\mathrm{\right)}}\mathrm{\right]}\mathrm{=}\mathrm{0.990}$ is obtained in this case (see Table 3). In addition, we investigated data from an adiabatic turbulence simulation (Schmidt et al. 2007), in which the rms Mach number gradually decreases with time because of the dissipative heating of the gas. The results are summarized in Table 3, and the correlation diagram for the final snapshot of the simulation is shown in Fig. 5b. Our results suggest that the closure coefficients are not very sensitive to the Mach number. Nevertheless, we cannot exclude that the optimal values of C₁ and C₂ differ significantly if the SGS turbulence Mach number ℳ_sgs (see Sect. 2) is only a tiny fraction of the speed of sound. Answering this question is left for future studies.

Fig. 6 Visualization of the SGS turbulence energy density Ksgs in a 5123 LES with solenoidal forcing.

3.4. Energy dissipation

For the turbulence energy on the length scale Δ_G, which is defined by Eq. (20), a dynamical equation analogous to Eq. (7) can be formulated. By averaging this equation over the whole periodic domain and assuming statistical equilibrium, i.e., ∂_t ⟨ K_ΔG ⟩  ≃ 0, the following global balance equation is obtained: $\mathrm{⟨}{\mathrm{\Sigma }}_{{\mathrm{\Delta }}_{\mathrm{G}}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\overline{\mathit{Pd}}\mathrm{-}\overline{\mathit{P}}\hspace{0.17em}\overline{\mathit{d}}\mathrm{⟩}\mathrm{-}{\mathit{C}}_{\mathit{ϵ}}⟨\frac{{\mathit{K}}_{{\mathrm{\Delta }}_{\mathrm{G}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\overline{\mathit{\rho }}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{\Delta }}_{\mathrm{G}}}⟩\mathrm{\simeq }\mathrm{0.}$(35)The first term is the mean turbulence energy flux, the second term is the mean pressure dilatation (see Schmidt et al. 2006a), and the third term is the mean dissipation rate expressed in terms of K_ΔG. Substituting Eq. (18) for Σ_ΔG yields the coefficient of turbulence energy dissipation, C_ϵ. From the supersonic isothermal turbulence data, we find a value C_ϵ ≈ 1.5, which is somewhat higher yet still comparable to typical values calculated for incompressible turbulence (see Sagaut 2006).

4. Large eddy simulations of forced supersonic turbulence

To investigate statistical properties of the SGS turbulence energy and related quantities, we run LES of forced supersonic isothermal turbulence with the SGS model defined in Sects. 2 and 3. For the implementation, we use the code Enzo 1.5 developed by the Laboratory for Computational Astrophysics at the University of California in San Diego (http://lca.ucsd.edu). In these simulations, we apply solenoidal, compressive and mixed force fields to produce statistically stationary and homogeneous turbulence with different rms Mach numbers (see Schmidt et al. 2006c; Schmidt et al. 2009; Federrath et al. 2010b). The forcing acts on length scales around the integral length L, where L is one half of the box size. The autocorrelation time of the force field is given by the time scale T = L/V, where the characteristic velocity V specifies the magnitude of the turbulent velocity fluctuations on the integral scale. The mixture of solenoidal (divergence-free) and compressive (rotation-free) modes of the force field is adjusted by means of a Helmholtz decomposition with weighing parameter 0 ≤ ζ ≤ 1. Purely solenoidal forcing results for ζ = 1. When setting the adiabatic exponent γ = 1.001, the energy dissipated per integral time is much lower than the internal energy for ℳ_rms up to about 10. For this reason, the gas is pseudo-isothermal. This approximate treatment of isothermality enables us to monitor energy conservation. With our implementation of the SGS model, the sum of resolved kinetic energy, SGS turbulence energy, and internal energy minus the power of the forcing integrated over time is conserved for the whole computational domain to a relative precision better than 10^-8. The fraction of computational time consumed by the SGS model is in the percent range.

Fig. 7 Visualization of the denstrophy Ω1/2 for the same snapshot as in Fig. 6.

4.1. Correlations with resolved flow quantities and the effective pressure

As an example, Fig. 6 shows a visualization of K_sgs prepared from an LES with 512³ grid cells. The parameters of this simulation were chosen to match the ILES with solenoidal forcing in Sect. 3.1. The rms Mach number of the flow is about 5.5 in the statistically stationary regime. In the reddish regions, K_sgs is higher than the spatial mean, while it is lower in the bluish regions. For comparison, Fig. 7 shows the local denstrophy ${\mathrm{\Omega }}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\left|\mathbf{\nabla }\mathrm{\times }{}^{\mathrm{(}}\mathit{\varrho }^{\mathrm{1}\mathit{/}\mathrm{2}}{u}{}^{\mathrm{)}}\right|}^{\mathrm{2}}$, which is an indicator of compressible turbulent velocity fluctuations (Kritsuk et al. 2007). It appears that high SGS turbulence energy is concentrated in regions of intense denstrophy. On average, K_sgs ~ 0.1Δ²Ω_1/2 for high denstrophy values, as one can see in the correlation diagram of K_sgs vs. Δ²Ω_1/2 in Fig. 8a. The same relation is found for compressive forcing (see panel (b) of Fig. 8). Nevertheless, the local values of K_sgs and Δ²Ω_1/2 deviate substantially from the average relation. This is a consequence of the various processes contributing to the SGS dynamics, which are not fully encompassed by the derivative of the resolved velocity field. For this reason, derived quantities such as the rate of strain or the denstrophy are only of limited utility for estimating the effects of turbulence on unresolved length scales. Since ${\mathit{P}}_{\mathrm{sgs}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{K}}_{\mathrm{sgs}}$, this applies also to the turbulent pressure.

The phase diagrams of the effective pressure (6) vs. the mass density are plotted in Fig. 9 for both LES. One can see that the average of the effective pressure for a given mass density closely follows the isothermal relation P ∝ ρ. This is because the mean turbulent pressure is small compared to the thermal pressure for the resolution Δ = L/256 (see Sect. 4.2). Locally, however, the intermittency of turbulent velocity fluctuations can give rise to an effective pressure that exceeds the thermal pressure by one order of magnitude. For this reason, the contribution of the turbulent pressure P_sgs is locally not negligible. This effect becomes stronger as the cutoff scale Δ increases in comparison to the integral scale of turbulence.

In Sect. 2, we argue that the viscous stress term in the filtered momentum Eq. (3) vanishes in the limit of infinite Reynolds number, and the rate of energy dissipation on the grid scale, ϵ, is determined by the SGS turbulence energy (see Eq. (10)). In contrast, an extrapolation of the expression for the microscopic dissipation rate on length scales ℓ ~ η to the grid scale Δ was proposed by Pan et al. (2009): $\mathit{\varrho ϵ}\mathrm{=}\mathit{\varrho }{\mathit{\nu }}_{\mathrm{\Delta }}\mathrm{|}{\mathit{S}}^{\mathrm{\ast }}{\mathrm{|}}^{\mathrm{2}}\mathit{.}$(36)The grid-scale viscosity ν_Δ = const. in the above expression is treated as a constant coefficient that is determined by the mean numerical dissipation of PPM. Defining the compressible Reynolds number of the resolved flow by $\mathrm{R}{\mathrm{e}}_{\mathrm{\Delta }}\mathrm{=}\mathrm{2}{\mathit{L}}^{\mathrm{2}}\mathrm{⟨}\mathrm{|}{\mathit{S}}^{\mathrm{\ast }}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}{\mathit{u}}_{\mathrm{rms}}^{\mathrm{2}}$, where u_rms is the root mean square velocity⁵, the viscosity can be evaluated from ν_Δ = VL/Re_Δ. The problem with this approach is that the viscosity on the grid scale, which corresponds to the SGS eddy-viscosity, cannot be assumed to be constant.

Fig. 8 Correlation diagrams of the SGS turbulence energy vs. the denstrophy, normalized by the cutoff scale Δ, for 5123 LES with solenoidal and compressive forcing. The contours are logarithmic. The average relation between both quantities is indicated by the dotted lines, and the dashed line shows the relation Ksgs ~ 0.1Δ2Ω1/2.

Fig. 9 Phase diagrams of the effective pressure defined by Eq. (6) vs. the mass density for 5123 LES with solenoidal and compressive forcing. The contours are logarithmic. The averages of the SGS turbulence energy for particular values of the denstrophy are indicated by the dotted lines.

Fig. 10 Correlation diagrams of the normalized rate of energy dissipation defined by Eq. (36), where νΔ assumes a constant value that is given by the numerical Reynolds number, vs. the rate of energy dissipation (10) that is predicted by the SGS model. The averages of expression  (36) for given values of ϱϵsgs are indicated by the dotted lines.

When neglecting diffusion, compressibility, and the nonlinear term in the SGS turbulence stress (8), the equilibrium between production and dissipation of SGS turbulence energy in Eq. (7) implies K_sgs ~ (C₁/C_ϵ)ϱΔ²|S^∗|², hence, ϵ ~ (Δ/C_ϵ)²(C₁|S^∗|)³ according to Eq. (10). We emphasize that a relation of the form ϵ ~ Δ²|S^∗|³ follows from any common SGS model under the assumption of local equilibrium (Sagaut 2006). Comparing to Eq. (36), we see that that ν_Δ ~ Δ²|S^∗|, which is not a constant. This is because ϱϵ ≠ σ_iju_i,j ∝ |S^∗|², as explained in Sect. 2. The discrepancy becomes apparent in Fig. 10, which shows the correlation diagrams of the rate of energy dissipation calculated via Eq. (36) vs. ϵ following from the SGS model. Toward low values of ϵ, we find an average relation close to |S^∗|² ∝ ϵ^2/3, which is just the relation that follows from the above estimate of the equilibrium dissipation rate. This behaviour is reasonable because the contribution of the nonlinear term in the closure (8), which is neglected in the estimate, is relatively small for low values of K_sgs (corresponding to low energy dissipation). Moreover, the unresolved velocity fluctuations tend to be significantly smaller than the speed of sound in this limit, which corresponds to low compressibility. Consequently, the results from the LES support the theoretical arguments against Eq. (36) as an approximation to the dissipation rate on the grid scale. Although Eqs. (10) and (36) yield about the same mean dissipation rate, the former determines the local rate of energy dissipation on the footing of a physically well motivated scale-separation of fluid dynamics, while the latter is based on a putative analogy between the numerical and the microscopic viscosity.

4.2. Dependence on the cutoff scale

The scaling of the turbulent velocity fluctuations in supersonic hydrodynamic turbulence has been inferred from energy spectrum functions and structure functions (e.g., Kritsuk et al. 2007; Schmidt et al. 2008, 2009; Federrath et al. 2010b; Price & Federrath 2010). The pure velocity scaling in the supersonic regime is stiffer than Kolmogorov scaling, and it appears that the scaling exponent depends on the forcing. For example, Federrath et al. (2010b) find indices of the turbulence energy spectra β = −1.86 ± 0.05 and −1.94 ± 0.05 for solenoidal and compressive forcing, respectively. For incompressible turbulence, β = −5/3. Velocity variables with fractional massweighing, in particular ρ^1/3u, exhibit similar scaling laws, which can be interpreted as indicating universality (Kritsuk et al. 2007; Schmidt et al. 2008).

The SGS turbulence energy is given by the fluctuations in the velocity and density fields on length scales ℓ ≲ Δ, as defined by Eqs. (4) and (5). However, there is no obvious relation to the known scaling laws for turbulence, because the decomposition of the fluid dynamical variables cannot be related to the two-point statistics (structure functions) or the Fourier modes (energy spectra) in a straightforward manner. To determine the scaling of SGS turbulence as a function of Δ, we run several LES with Δ ranging from L/256 to L/32. The mean values of the rms Mach number and the SGS turbulence Mach number are plotted as functions of time in Fig. 11. The flow approaches a statistically stationary state after about 2 integral time scales (Schmidt et al. 2009; Federrath et al. 2010b), for which ℳ_rms settles at values between 5 and 6. The temporal variation in ℳ_rms is caused by the stochastic forcing. As expected, $\mathrm{⟨}{\mathrm{ℳ}}_{\mathrm{sgs}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ decreases with the cutoff scale. Averaging the spatial means from t = 2T to 10T, we find the time-averaged mean values listed in Table 4. As one can see in Fig. 12a, the time averages of $\mathrm{⟨}{\mathrm{ℳ}}_{\mathrm{sgs}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ closely follow power laws, $\mathrm{⟨}{\mathrm{ℳ}}_{\mathrm{sgs}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\propto }{\mathrm{\Delta }}^{{\mathit{\alpha }}_{\mathrm{ℳ}}}\mathit{,}$(37)with α_ℳ = 0.475 ± 0.004 for solenoidal and 0.451 ± 0.026 for compressive forcing.

Fig. 11 Temporal evolution of the rms Mach number (top) and the mean SGS turbulence Mach number (bottom) for solenoidal (left column) and compressive forcing (right column). The cutoff length Δ decreases from L/32 (light colour) to L/256 (full colour).

Fig. 12 Scaling laws for the mean SGS turbulence Mach number a) and energy b) as functions of the numerical resolution Δ.

Fig. 13 Temporal evolution of the SGS turbulence energy (top) and the dissipation rate (bottom) for solenoidal (left column) and compressive forcing (right column). The cutoff length Δ decreases from L/32 (light colour) to L/256 (full colour).

The behaviour of the mean SGS turbulence energy is similar, although the intermittent fluctuations of  ⟨ K_sgs ⟩  are more pronounced than the mean SGS turbulence Mach number (see top panels of Fig. 13). The higher degree of intermittency stems from the mass density that is included in K_sgs. For compressive forcing,  ⟨ K_sgs ⟩  is systematically lower in comparison to the LES with solenoidal forcing. This indicates that the total amount of energy in the turbulent structures on a given length scale is smaller in the compressive forcing case. The ratio of the mean values of K_sgs for compressive and solenoidal forcing in Table 4 approximately agrees with the ratio 0.38 that is inferred from the filtered high-resolution data. On the other hand, the scaling laws $\mathrm{⟨}{\mathit{K}}_{\mathrm{sgs}}\mathrm{⟩}\mathrm{\propto }{\mathrm{\Delta }}^{{\mathit{\alpha }}_{\mathit{K}}}\mathit{,}$(38)are nearly the same for solenodial and compressive forcing (see Fig. 12b). We find the slopes α_K = 0.799 ± 0.009 and 0.769 ± 0.029, which agree within the error bars. This result is remarkable, because it suggests that the scaling properties of turbulence on small length scales are independent of the forcing.

As can be seen in the bottom panels of Fig. 13, the above scaling law of  ⟨ K_sgs ⟩  results in a mean dissipation rate  ⟨ ρϵ ⟩  that is independent of the cutoff scale, which is an essential property of the energy dissipation predicted by the SGS model. The time-averaged mean values are listed in Table 4. Since turbulence is statistically stationary in these simulations, the constant rate of energy dissipation corresponds to a mean turbulence energy flux that is also constant. Theoretical arguments implying the existence of a scale-free energy flux in the compressible turbulent cascade have recently been presented by Aluie (2011). The significantly lower mean dissipation rate in the case of compressive forcing is consistent with the energy spectra of ρ^1/3u (see Fig. A.1, in Federrath et al. 2010b). This mass-weighted velocity variable is related to the energy dissipation rate (Kritsuk et al. 2007). Moreover, Fig. 14 (left panel) shows that the growth of the mean internal energy in time becomes smaller as the weighing parameter ζ decreases from 1 (solenodial forcing) to 0 (compressive forcing). After subtracting the contribution from numerically resolved compression effects (see Schmidt et al. 2006c), we find that, independent of the cutoff length Δ, about 3/4 of the change of the internal energy stems from SGS turbulence energy dissipation. The remainder is caused by numerical dissipation. This does not imply that the total rate of energy dissipation is much higher in LES than in ILES, because the total energy dissipation is always determined by the energy injection due to the forcing. In conclusion, the greater part of kinetic energy is dissipated through the SGS turbulence energy reservoir at a scale-free rate.

To quantify the relative importance of high values of ℳ_sgs, we determined the volume fractions of cells with an SGS turbulence Mach number greater than a particular value. This fraction is given by 1 − cdf(ℳ_sgs), where cdf(ℳ_sgs) is the cumulative distribution function of ℳ_sgs. In Fig. 15, the resulting functions are plotted for the LES with different cutoff lengths. As expected, the fraction with ℳ_sgs > 1 decreases with the cutoff length Δ. However, the tails toward high ℳ_sgs demonstrate that even at relatively high resolution there are supersonic velocity fluctuations on unresolved length scales, and the corresponding turbulent pressure decreases only a little with the cutoff scale.

Fig. 14 Left panel: time evolution of the mean internal energy for forcing with varying ζ and about the same rms Mach number. Right panel: volume fractions of zones, in which the SGS turbulence Mach number is greater than any given value for the same forcing parameters as in the left panel (the curves for ζ = 1 and ζ = 2/3 almost coincide).

Fig. 15 Volume fractions of zones, in which the SGS turbulence Mach number is greater than a certain value for different numerical resolutions. The dashed lines follows from the filtering of 10243 ILES data with filter length 16Δ. This corresponds to the 643 LES, for which Δ/L = 1/32.

For the lowest resolution LES, we can compare the distribution of ℳ_sgs to the distribution inferred from the corresponding filtered 1024³ data (see Sect. 3.1). The filter length Δ_G = 16Δ = L/32 is equivalent to the cutoff length in the 64³ LES. Choosing an even lower resolution of the LES, corresponding to a larger filter length for the ILES, turned out not to be feasible. Even for the LES with Δ/L = 1/32, the forcing range and the range of length scales that are directly affected by numerical dissipation overlap. For the filtering of the ILES, on the other hand, the filter length cannot be lowered (corresponding to a higher resolution of the LES), because the dynamical range of fluctuations between the grid scale and the filter length would become insufficient and the numerical smoothing would be too strong. Nevertheless, Fig. 15a demonstrates that the distributions agree remarkably well in the case of solenoidal forcing. For compressive forcing, there are larger discrepancies. However, given that Gaussian filtering only corresponds roughly to the implicit filter in an LES and that the SGS model is based on various approximations, the match is quite satisfactory.

The larger deviations in the case of compressive forcing suggest that it is not possible to calibrate the SGS model coefficients in such a way that an optimal match is obtained both for solenoidal and for compressive forcing at the same time. The different shape of the distribution that is obtained from the high-resolution simulation with compressive forcing points toward a missing physical effect such as the pressure-dilatation, which is entirely neglected in our SGS model. Either way, purely compressive forcing is a limiting case. In nature, some mixture of solenoidal and compressive forcing is more likely to occur. In Fig. 14 (right panel), we compare the distributions of ℳ_sgs for force fields with ζ varying from 1 (solenoidal) to 0 (compressive). High SGS turbulence Mach numbers become more frequent as the contribution of compressive forcing modes increases.

4.3. Dependence on the Mach number

At fixed resolution, the SGS turbulence energy increases with the resolved kinetic energy of the flow. For isothermal turbulence, this also implies an increase in the SGS turbulence Mach number with rising rms Mach number. To investigate this dependence, we varied the magnitude of mixed forcing with ζ = 2/3 and ζ = 1/3. In the case ζ = 1/2, a forcing field with two solenoidal and one longitudinal degrees of freedom is obtained. Figure 8 in Federrath et al. (2010b) demonstrates that the ratio of the energy that is contained in transversal and longitudinal modes approaches a constant value for ζ > 1/2 (see also Kritsuk et al. 2010). Correspondingly, the cumulative distributions of ℳ_sgs plotted in Fig. 14b show that there is almost no difference between forcing with ζ = 2/3 and purely solenoidal forcing (ζ = 1). On the other hand, one can see that there is a noticeable influence of compressive modes for ζ = 1/3, but the distribution differs from the purely compressive case (ζ = 0).

The dependence of the RMS SGS turbulence Mach number on the forcing amplitude is shown in the top panels of Fig. 16 for 256³ LES. The corresponding characteristic Mach numbers and time-averaged statistics are listed in Table 5. Independent of the Mach number, the ratio of $\mathrm{⟨}{\mathrm{ℳ}}_{\mathrm{sgs}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ to ℳ_rms is nearly equal for ζ = 2/3 and 1/3. This ratio is 0.102 for purely solenoidal forcing (ζ = 1), and 0.138 for purely compressive forcing (ζ = 0). Consequently, $\mathrm{⟨}{\mathrm{ℳ}}_{\mathrm{sgs}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{/}{\mathrm{ℳ}}_{\mathrm{rms}}$ is mostly determined by the grid scale Δ, except for small ζ. As one can see in the middle panels of Fig. 16, the normalized mean SGS turbulence energy,  ⟨ K_sgs ⟩ /ρ₀V², is about the same for the different Mach numbers, with a weak trend to decrease toward low Mach numbers (see also Table 5). The same behaviour is found for the mean dissipation rate (Fig. 16, bottom panels), which further supports the validity of the SGS model in the supersonic regime. Following the trend discussed in Sect 4.2, clearly the mixture of solenoidal and compressive forcing modes has an influence.

Figure 17 shows the distributions of ℳ_sgs for ζ = 2/3 (a) and 1/3 (b), as explained in Sect. 4.2. For both forcing types, the volume fractions with ℳ_sgs > 1 increase with the forcing magnitude. For rms Mach numbers greater than 5, supersonic turbulent velocity fluctuations at the cutoff scale fill more then 10% of the total volume. If the ratio between the integral scale of turbulence and the cutoff scale is less, the volume filling factor increases further.

Fig. 16 Temporal evolution of the SGS turbulence Mach number (top) and energy (middle) and the dissipation rate (bottom) for ζ = 2/3 (left column) and ζ = 1/3 (right column). The different lines in each plot correspond to different forcing magnitudes (see Table 5), which are specified by the values of the characteristic Mach number Ma = V/c0 (c0 is the initial speed of sound).

5. Conclusion

In formulating a mixed closure for the flux of energy from the numerically resolved to the unresolved scales, we generalized the subgrid-scale (SGS) turbulence energy model to the regime of highly compressible turbulence. This closure is based on ideas of Woodward et al. (2006) and features a nonlinear term in addition to a linear eddy-viscosity term. In general, the turbulence energy cascade is an important source of SGS turbulence energy production in turbulent flows, and it should not be neglected even if other subresolution sources are emphasized in particular astrophysical applications (e.g., Scannapieco & Brüggen 2008; Joung et al. 2009). Our proposed closure for the transfer of energy by the turbulent cascade complements theses models. We verified this closure by means of explicit filtering of high-resolution data from various simulations of supersonic isothermal and adiabatic turbulence (Schmidt et al. 2007, 2009; Federrath et al. 2010b). Tests in large eddy simulations of forced supersonic turbulence show that the SGS model meets several important requirements that should be satisfied by any sound SGS model:

• for statistically stationary turbulence, an equilibrium betweenthe production and the dissipation of SGS turbulence is reached.The mean SGS turbulence energy depends on the grid scale via apower law;

• the SGS turbulence Mach number, which specifies the importance of the turbulent pressure on the grid scale relative to the thermal pressure of the gas, depends linearly on the rms Mach number of the resolved turbulence;

• the SGS turbulence energy dissipation is independent of the grid scale.

Forced turbulence simulations in a periodic box are very suitable to testing the properties listed above, because of the well-defined statistics of the isotropic, homogeneous and stationary turbulence that is produced. In addition to these properties, we find a dependence of the SGS turbulence energy and the rate of energy dissipation on the mixture of solenoidal (divergence-free) and compressive (rotation-free) forcing modes. However, the scaling laws for the SGS turbulence energy are very similar for solenoidal and compressive forcing. This indicates that the SGS model describes the dynamics in the inertial subrange, although the length scales close to the cutoff scale are affected by numerical dissipation. The differences in the mean values result from the substantial differences in the turbulent flow structure on larger scales (see Federrath et al. 2008, 2009, 2010b).

Fig. 17 Volume fractions of zones, in which the SGS turbulence Mach number is greater than any given value for different rms Mach numbers and mixtures of solenoidal and compressive modes, as in Fig. 16.

The implementation of the SGS model into a fluid-dynamical code such as Enzo is rather straightforward. The SGS turbulence energy can be treated as a passive scalar with various source terms. To evaluate the mixed closure for the turbulence energy cascade, derivatives of the resolved velocity field are computed by means of centred differences. Care must be taken to ensure energy conservation, particularly if the net change of the SGS turbulence energy in a certain grid cell exhausts the available energy over a time step or, vice-versa, if too much energy is drained from the resolved scales. However, these are exceptions that can be handled numerically. The effective increase in computing time is less than 10%. Moreover, the hydrodynamic coupling of the SGS model to the resolved flow introduces a reduction of the bottleneck effect in the turbulence energy spectra (Woodward et al. 2006; Schmidt 2010).

In our large eddy simulations, we find correlations between resolved flow quantities, such as the rate of strain or the denstrophy and SGS quantities, but the scatter is large. This is problematic if one intends to estimate unresolved flow properties on the basis of such correlations. In particular, this applies to the calculation of the dissipation rate from a constant numerical viscosity and the rate of strain, as proposed by Pan et al. (2009). We have shown that the assumption of a constant numerical dissipation coefficient is inconsistent with the equilibrium relation between the dissipation rate and the rate of strain on the grid scale in the limit of large Reynolds numbers. This relation, which follows from the SGS turbulence energy model with the linear eddy-viscosity closure (and also from the Smagorinsky model), is verified by our LES data for low turbulence intensity, while deviations become apparent for strong turbulent dissipation. This can be understood as a consequence of the nonlinear term in the closure for the turbulence energy flux. Also the estimate of turbulent pressure effects using the rate of strain and the vorticity of the resolved flow that is put forward by Zhu et al. (2010) is incomplete because they do not distinguish between the contributions from the resolved flow and from the subgrid scales. The predictions from both approaches with regard to turbulence in the intergalactic medium are compared in ongoing work (Iapichino et al. 2011).

From the probability distributions of the SGS turbulence Mach number, it follows that the turbulent pressure locally exceeds the thermal pressure even at moderate rms Mach numbers and for relatively small grid scales. Since the grid scale in contemporary galactic disk simulations (e.g., Agertz et al. 2009; Tasker & Tan 2009) is close to molecular cloud scales (a few pc), unresolved supersonic velocity fluctuations are quite likely and the turbulent pressure plays an important role. This has implications for the treatment of collapsing gas regions. The criterion for gravitational stability, which influences the grid resolution in adaptive mesh refinement simulations and controls the production of sink particles to capture the collapsing gas, is usually based on the thermal Jeans mass (among other criteria; see Federrath et al. 2010a). To account for the effects of turbulence below the grid scale, we suggest including the turbulent pressure in the definition of the Jeans mass, in analogy to the magnetic pressure in self-gravitating MHD turbulence.

To model the fragmentation below the grid scale in more detail, a possible approach could be based on the assumption that the local star formation efficiency is regulated by turbulence on the grid scale. Then the star formation rate can be parameterized in terms of the turbulent Mach number that is calculated from the SGS model (see Krumholz & McKee 2005; Padoan & Nordlund 2009). On the other hand, star formation acts back on the SGS turbulence energy via stellar feedback. As suggested by Joung et al. (2009), a stellar feedback term can be included in the SGS turbulence energy equation. Statistically, we have

$\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{K}}_{\mathrm{sgs}}\mathrm{~}\mathrm{\Sigma }\mathrm{+}{\mathrm{\Sigma }}_{\mathit{\star }}\mathrm{-}\mathit{\rho ϵ,}$where Σ_⋆ ∝ ρe _⋆ /τ_ff accounts for the energy injection per unit mass, e_⋆, by supernovae. The associated time scale is the free-fall time scale τ_ff =  [3π/(32Gρ)] ^1/2, which is the fundamental time scale of star formation. Neglecting the fluctuations of the gas density and setting the mean production rate Σ ~ V²/T, where V and T are the typical velocity and the turn-over time scale, respectively, of the resolved turbulent flow, it follows that the turbulent pressure in equilibrium is of the order ${\mathit{P}}_{\mathrm{sgs}}\mathrm{~}\mathit{\rho }{\mathrm{\Delta }}^{\mathrm{2}\mathit{/}\mathrm{3}}{\left[\frac{{\mathit{V}}^{\mathrm{2}}}{\mathit{T}}\mathrm{+}{\mathit{C}}_{\mathit{\star }}\frac{{\mathit{e}}_{\mathit{\star }}}{{\mathit{\tau }}_{\mathrm{ff}}}\right]}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{·}$(39)Kolmogorov scaling becomes manifest in the factor Δ^2/3 in Eq. (39). For highly compressible turbulence, however, the scaling of the SGS turbulence energy deviates from the Kolmogorov law (see Sect. 4.2). Depending on the ratios V²/e _⋆  and T/τ_ff, the production of SGS turbulence energy by the turbulent cascade or by supernovae dominates. The model of Joung et al. (2009) follows in the limit (V²/e _⋆ )(τ_ff/T) ≪ 1. In general, shear instabilities, gravitational instabilities, cooling instabilities, etc. above the grid scale feed energy to smaller scales. For disk galaxies, a simple estimate can be obtained from the velocity dispersion of atomic hydrogen, which is about 10   km   s^-1. Agertz et al. (2009) show that gravitational instabilities grow on length scales ranging from 0.1 to about 2   kpc. Setting V = 10⁶   cm   s^-1 and assuming L > 0.1   kpc ≈ 3 × 10²⁰   cm, the turbulence energy flux to smaller length scales is V²/T = V³/L ≲ 0.003   erg   g^-1   s^-1. On the other hand, ${\mathit{e}}_{\mathit{\star }}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{48}}\hspace{0.17em}\mathit{\eta }\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}^{-1}\mathrm{\approx }\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{15}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{g}}^{-1}$ and C _⋆  ≈ 0.025 imply Σ _⋆  ~ 0.03η   (n/1   cm^-3)^1/2   erg   g^-1   s^-1 (see Joung et al. 2009). Since the efficiency of the energy transfer from supernova blast waves to the interstellar gas is roughly η ≃ 0.1 (Mac Low & Klessen 2004), the energy injection by stellar feedback is comparable to the turbulence energy flux for atomic hydrogen with density n ~ 1   cm^-3. The dependence of Σ _⋆  on the gas density implies a greater contribution from supernova feedback in the cold-gas phase, but the above estimate does not account for the intermittency of turbulent velocity fluctuations, which entails large deviations from the mean. As a consequence, the assumption of Joung et al. (2009) to consider the energy injection by supernova as the main source of the turbulent pressure is marginally fulfiled in cosmological simulations, in which the internal structure of galaxies is very poorly resolved. In galactic-scale simulations with high resolution, on the other hand, turbulence is not uniformly produced, and including the turbulence energy cascade improves the description of numerically unresolved processes. In particular, it will be useful to attempt a further generalization of the SGS model to multiphase turbulence. A very simple ansatz has recently been presented by Murante et al. (2010). A complete SGS model that treats a warm and a cold gas phase is presently under development (Braun & Schmidt, in prep.).

Including stellar feedback and cooling into our SGS model will be of further utility for the numerical treatment of turbulence in the intergalactic medium (see Springel & Hernquist 2003), where turbulence is produced by different processes (see, for instance, Cen & Ostriker 1999; Subramanian et al. 2006; Ryu et al. 2008; Iapichino et al. 2008, 2011). Oppenheimer & Davé (2009) show that a significant amount of the line broadening of O VI in cosmological simulations stems from numerically unresolved turbulence. They apply a heuristic model in the postprocessing of the simulation data. When using an SGS model, on the other hand, the effect on the line broadening can be computed on the fly. Moreover, metals are mixed into the intergalactic medium by turbulence that is driven by galactic outflows. SGS turbulence enhances the turbulent mixing. Following an approach that is quite similar to the treatment of stellar feedback in galaxy simulations, our SGS model can be used in combination with phenomenological models for supernova-driven outflows (Joung et al. 2009; Evoli & Ferrara 2011). In both cases, the use of adaptive mesh refinement is mandatory for achieving a sufficient dynamical range. Therefore, an essential objective for future work will be to incorporate the new closure for the highly compressible turbulent cascade into fluid mechanics with adaptively refined large eddy simulations (FEARLESS; Maier et al. 2009).

Acknowledgments

We thank Jens Niemeyer and Ralf Klessen for valuable discussions. Computations described in this work were performed using the Enzo code developed by the Laboratory for Computational Astrophysics at the University of California in San Diego (http://lca.ucsd.edu). The computational resources were provided by the HLRBII project h0972 at the Leibniz Supercomputer Centre in Garching, Germany. C.F. acknowledges funding from the Landesstiftung Baden-Würrtemberg via their programme International Collaboration II (grant P-LS-SPII/18), from the German Bundesministerium für Bildung und Forschung via the ASTRONET project STAR FORMAT (grant 05A09VHA), from the International Max Planck Research School for Astronomy and Cosmic Physics (imprs-a) and from the Heidelberg Graduate School of Fundamental Physics (hgsfp), which is funded by the Excellence Initiative of the Deutsche Forschungsgemeinschaft (dfg) gsc 129/1. Furthermore, C.F. received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 247060) for the research presented in this work.

References

All Tables

All Figures

	Fig. 1 Correlation diagrams for the SGS turbulence energy flux in the case of isothermal supersonic turbulence with solenoidal forcing (ℳrms ≈ 5.3). The applied filter length is 32Δ. The blue dots indicate the average prediction of the closure for a given value of Σ32Δ.
In the text

	Fig. 2 Correlation diagrams for isothermal supersonic turbulence with compressive forcing (ℳrms ≈ 5.6) as in Fig. 1.
In the text

	Fig. 3 Slices of the two-dimensional probability density functions plotted in Figs. 1 and 2, showing the predictions of different closures for the values of the explicitly computed energy flux Σ32Δ that are indicated by the vertical dashed lines. The top and bottom rows of panels show the results for solenoidal and compressive forcing, respectively.
In the text

	Fig. 4 Correlation diagrams for the mixed closures for isothermal supersonic turbulence with different forcing.
In the text

	Fig. 5 Correlation diagrams for the mixed closure in the case of isothermal a) and adiabatic b) turbulence with lower rms Mach numbers as in Fig. 4.
In the text

	Fig. 6 Visualization of the SGS turbulence energy density Ksgs in a 5123 LES with solenoidal forcing.
In the text

	Fig. 7 Visualization of the denstrophy Ω1/2 for the same snapshot as in Fig. 6.
In the text

	Fig. 8 Correlation diagrams of the SGS turbulence energy vs. the denstrophy, normalized by the cutoff scale Δ, for 5123 LES with solenoidal and compressive forcing. The contours are logarithmic. The average relation between both quantities is indicated by the dotted lines, and the dashed line shows the relation Ksgs ~ 0.1Δ2Ω1/2.
In the text

	Fig. 9 Phase diagrams of the effective pressure defined by Eq. (6) vs. the mass density for 5123 LES with solenoidal and compressive forcing. The contours are logarithmic. The averages of the SGS turbulence energy for particular values of the denstrophy are indicated by the dotted lines.
In the text

	Fig. 10 Correlation diagrams of the normalized rate of energy dissipation defined by Eq. (36), where νΔ assumes a constant value that is given by the numerical Reynolds number, vs. the rate of energy dissipation (10) that is predicted by the SGS model. The averages of expression  (36) for given values of ϱϵsgs are indicated by the dotted lines.
In the text

	Fig. 11 Temporal evolution of the rms Mach number (top) and the mean SGS turbulence Mach number (bottom) for solenoidal (left column) and compressive forcing (right column). The cutoff length Δ decreases from L/32 (light colour) to L/256 (full colour).
In the text

	Fig. 12 Scaling laws for the mean SGS turbulence Mach number a) and energy b) as functions of the numerical resolution Δ.
In the text

	Fig. 13 Temporal evolution of the SGS turbulence energy (top) and the dissipation rate (bottom) for solenoidal (left column) and compressive forcing (right column). The cutoff length Δ decreases from L/32 (light colour) to L/256 (full colour).
In the text

	Fig. 14 Left panel: time evolution of the mean internal energy for forcing with varying ζ and about the same rms Mach number. Right panel: volume fractions of zones, in which the SGS turbulence Mach number is greater than any given value for the same forcing parameters as in the left panel (the curves for ζ = 1 and ζ = 2/3 almost coincide).
In the text

	Fig. 15 Volume fractions of zones, in which the SGS turbulence Mach number is greater than a certain value for different numerical resolutions. The dashed lines follows from the filtering of 10243 ILES data with filter length 16Δ. This corresponds to the 643 LES, for which Δ/L = 1/32.
In the text

	Fig. 16 Temporal evolution of the SGS turbulence Mach number (top) and energy (middle) and the dissipation rate (bottom) for ζ = 2/3 (left column) and ζ = 1/3 (right column). The different lines in each plot correspond to different forcing magnitudes (see Table 5), which are specified by the values of the characteristic Mach number Ma = V/c0 (c0 is the initial speed of sound).
In the text

	Fig. 17 Volume fractions of zones, in which the SGS turbulence Mach number is greater than any given value for different rms Mach numbers and mixtures of solenoidal and compressive modes, as in Fig. 16.
In the text