1 Introduction

Turbulent flows with high Reynolds numbers are often encountered in computational astrophysics. Examples are the solar wind, stellar convection zones, star-forming clouds, and probably the gas in galaxy clusters. This review concentrates on computational methods that treat turbulence in the limit of high Reynolds numbers by explicitly solving the compressible Euler equations for the large-scale dynamics of the flow, while incorporating small-scale effects such as viscous dissipation into a subgrid-scale model. Since the non-linear turbulent interactions between different scales are at least partially resolved, this type of simulation is called large eddy simulation (LES).

The relative importance of non-linear interactions and viscous damping is specified by the Reynolds number. It is determined by the characteristic velocity V of the flow, its integral length scale L, and the microscopic viscosity ν:

$$\operatorname{Re} = \frac{{VL}} {\nu }.$$
(1)

The flow becomes turbulent if the non-linear interactions are much stronger than viscous damping. Generally, this happens if Re reaches values greater than a few 103 for flows with boundaries. But Re in astrophysical flows are typically much greater than that. For instance, an estimate for the turbulent convection zone of the Sun is Re ∼ 1014 (Canuto, 1994).

In principle, we can also define a scale-dependent Reynolds number Re() = υ′()/ν, where ν′() is the typical magnitude of velocity fluctuations on the length scale . The length sale of strong viscous damping is then given by Re(K) ∼ 1. For incompressible turbulence, substitution of the Kolmogorov-Obukhov scaling law υ′() ∼ (εℓ)1/3 yields (Frisch, 1995)

$${{{\epsilon ^{1/3}}\ell _{\rm{K}}^{4/3}} \over \nu }\sim 1.$$

. Since the mean dissipation rate is εV3/L, it follows that

$$\frac{L} {{\ell K}} \sim \operatorname{Re} ^{3/4} .$$
(2)

The problem of high Re is thus a problem of largely different length scales or, equivalently, a high number of degrees of freedom.

In a numerical simulation of turbulence, the range of length scales is limited by the grid scale Δ, which is simply the linear size of the grid cells. Only if Δ ≲ K, turbulence can be fully resolved by a so-called direct numerical simulation (DNS).Footnote 1 However, DNS become infeasible for very large Re because the total amount of floating point operations (FLOPs) increases with (L/Δ)4 ≳ (L/K)4 ∼ Re3. The scaling may differ for highly compressible turbulence, but the basic problem remains the same. For a DNS of solar convection over one dynamical time scale, it would be necessary to perform very roughly 1042 FLOPs, which would take far longer than the current age of the Universe on the fastest existing computer.

In practice, however, it is neither feasible nor absolutely necessary to account for all degrees of freedom in a simulation of high-Re turbulence. To reproduce certain statistical properties, a much coarser sampling of the degrees of freedom can be quite sufficient for many applications. This is why LES encompass only the energy-containing scales and structures dominated by non-linear interactions, which are part of the turbulent cascade down to a cutoff scale much greater than the microscopic dissipation scale. The cutoff scale is given by grid scale Δ. The defining criterion for LES is thus L ≫ Δ ≫ K or, equivalently,

$$\operatorname{Re} \gg \operatorname{Re} (\Delta ) \gg 1.$$

Here, Re(Δ) ∼ υ′(Δ)Δ/ν is the Reynolds number of subgrid-scale turbulence. The product υ′(Δ)Δ can be interpreted as turbulent viscosity of the numerically unresolved eddies of size ≲ Δ. The effective Reynolds number of the numerically computed flow is therefore given by

$$\operatorname{Re} _{eff} = \frac{{\operatorname{Re} }} {{\operatorname{Re} (\Delta )}} \sim \frac{{VL}} {{v'(\Delta )\Delta }} \sim \left( {\frac{L} {\Delta }} \right)^{4/3} .$$
(3)

This means that LES reduces the number of degrees of freedom by replacing the microscopic viscosity ν by a turbulent viscosity of the order υ′(Δ)Δ ≫ ν. As a result, the purely non-linear turbulent dynamics of the “large eddies” is separated from microscopic dissipation.Footnote 2 The biggest challenge when implementing this concept is to find an appropriate model for the coupling between the small- and large-scale dynamics.

A mathematical framework for LES is based on the notion of a filter, which separates large-scale ( ≳ Δ) from small-scale ( ≲ Δ) fluctuations. Filters can be used to decompose the equations of fluid dynamics into equations for smoothed variables, which have a very similar mathematical structure as the unfiltered equations, and equations for second-order moments of the fluctuations. The latter are interpreted as subgrid-scale variables. In Section 2, we will carry out the decomposition of the compressible Navier-Stokes equation by applying the filter formalism of Germano (1992). This formalism comprises the so-called Reynolds-averaged Navier-Stokes (RANS) equations as limiting case if the filter length is comparable to the integral length scale of the flow. This method is equivalent to numerically solving a mean-field theory for turbulent flow. Simulations based on the RANS equations work with low Reeff, while LES have high Reeff. In principle, second-order moments can be expressed in terms of higher-order moments. Since this would entail an infinite hierarchy of moments, the set of variables is limited by introducing closures. Usually, one attempts to find closures for the second-order moments by expressing them in terms of the filtered variables. This is what is called a subgrid-scale (SGS) model.Footnote 3 For example, a complete second-order closure model for turbulent convection is formulated in Canuto (1994). Much simpler, yet often employed is the one-equation model for the SGS turbulence energy K, i.e., the local kinetic energy of numerically unresolved turbulent eddies. For this reason, it is sometimes called the K-equation model. Closures for the transport and source terms in the SGS turbulence energy equation are presented in some detail in Section 3, followed by a discussion of how the closure coefficients can be determined (Section 4). Of particular importance is the prediction of the local turbulent viscosity, which is is given by Δ-√K times a dimension-less coefficient. The turbulent viscosity is required to calculate the turbulent stresses, which enter the equations for the filtered variables analogous to the viscous stresses in the unfiltered Navier-Stokes equations (see Section 3.1).

Filtering the dynamical equations is usually considered to be equivalent to numerical discretization. The filter length can then be identified with the grid scale Δ. Since the numerical truncation errors of finite difference or finite volume schemes for the computation of compressible astrophysical flows are more or less diffusion-like terms, they produce a numerical viscosity that effectively reduces the Reynolds number to a value comparable to Eq. (3). It is actually a common assumption that numerical viscosity approximates the turbulent viscosity on the grid scale. This leads to the notion of an implicit large eddy simulation (ILES) (Garnier et al., 2009), which is widely used for simulating turbulent flows in astrophysics. Numerous numerical studies demonstrated that ILES is a very robust method, which reliably predicts scaling laws of compressible turbulence at sufficiently high resolution (Sytine et al., 2000; Kritsuk et al., 2007; Benzi et al., 2008; Schmidt et al., 2008; Federrath et al., 2010b; Kritsuk et al., 2013). This is a consequence of the independence of inertial-range scaling from the dissipation mechanism, be it microscopic, turbulent or numerical viscosity, provided that the dynamical range of the simulation is large enough. In simulations of statistically stationary isotropic turbulence, however, the inertial subrange is very narrow for computationally feasible resolutions because the bottleneck effect distorts the spectrum over a large range of high wave numbers below the Nyquist wavenumber (Falkovich, 1994; Dobler et al., 2003; Schmidt et al., 2006a). It appears that LES with an explicit SGS model, such as the K-equation model, can reduce the bottleneck effect to some degree and reproduce scalings from ILES or DNS at lower resolution (Haugen and Brandenburg, 2006; Woodward et al., 2006; Schmidt, 2010). However, more systematic studies covering the parameters space of forced compressible turbulence are necessary to confirm this effect.

There are, of course, alternative methods of scale separation and a large variety of SGS models (for a comprehensive overview, see the monographs Sagaut, 2006; Garnier et al., 2009). An example are the Camassa-Holm equations, which follow from the incompressible Navier-Stokes equations by decomposing the trajectories of fluid elements into mean and fluctuating parts in the Lagrangian framework (Chen et al., 1998). Since the filtered component of the velocity is defined by an inverse Helmholtz operator of the form (1 − α22)−1, which is explicitly applied to determine the turbulent stresses in the filtered velocity equation, the resulting model is called Lagrangian-averaged Navier-Stokes α-model (LANS-α). Depending on the choice of α, the variables computed in LES based on LANS-α are typically smoothed over length scales somewhat larger than the grid resolution. In other words, this type of simulation partially resolves the sub-filter scales, which improves the controllability of the model. While there is no handle on the competition between the SGS model and numerical truncation errors on the grid scale in convectional LES, LANS-α can, in principle, alleviate this problem by adjusting the balance between truncation and model errors (Pietarila Graham et al., 2007). Although the idea is very elegant, the numerical studies discussed in Pietarila Graham et al. (2007, 2008) show that the applicability of LANS-α and similar models is limited, particularly for very high Re. Moreover, the generalization to compressible turbulence is not straightforward. Models such as LANS-α are not further covered by this review, but they might be an option for magnetohydrodynamical LES (Pietarila Graham et al., 2009).

In astrophysics, LES are mainly applied to complex systems. In simulations of cosmological structure formation, which are discussed in Section 6.3, the length scales on which turbulence is driven by gravity are varying. Although adaptive mesh refinement is applied to track down collapsing structures, it is difficult to to resolve a wide range of length scales between the smallest driving scale and the grid scale at the highest refinement level. In this situation, SGS effects can become fairly large. However, the variable grid scale complicates the scale separation in AMR simulations because energy has to be transferred between the resolved and SGS energy variables if a region is refined or de-refined. Section 5.1 describes how to combine LES and AMR. This method, for which the acronym Fearless (Fluid mEchanics for Adaptively Refined Large Eddy SimulationS) was coined in Maier et al. (2009), has been applied to galaxy clusters, the intergalactic medium, and primordial atomic cooling halos. The results from these simulations indicate that the contribution of the numerically unresolved turbulent pressure to the support against gravity is non-negligible and the turbulent viscosity tends to stabilize disk-like structures around collapsed gas clouds. Moreover, the SGS model provides indicators of turbulence production and dissipation and allows for the computation of the turbulent velocity dispersion. A difficulty is that turbulence production by cosmological structure formation is highly inhomogeneous. This entails the problem that the SGS model should dynamically adapt to conditions ranging from laminar flow to developed turbulence. Inhomogeneous and non-stationary turbulence can be treated by dynamical procedures for the calculation of closure coefficients or shear-improved SGS models, which decompose the numerically resolved flow into mean and fluctuating components. These techniques are outlined in Sections 4.2 and 5.2.

Furthermore, SGS models offer unique possibilities for modeling physical processes that are influenced by turbulence. An example is turbulent deflagration, where the turbulent diffusivity predicted by the SGS model dominates the effective flame propagation speed in underresolved numerical simulations. Turbulent deflagration plays a role at least in the initial phase of thermonuclear explosions of white dwarfs (see Section 6.1), which is one of the scenarios that are thought to produce type Ia supernovae. A recent application along similar lines are LES of isolated disk galaxies, where the SGS turbulence energy is a crucial parameter for calculating the star formation rate and the feedback due to supernova blast wave (see Section 6.2). Since the impact of feedback processes on the formation of galaxies and their evolution leaves many questions unanswered, galaxies are a particularly promising field of application.

While great progress has been made for compressible hydrodynamics, magnetohydrodynamical LES are still in their infancy. Several SGS models have been proposed in the context of terrestrial plasma physics (Müller and Carati, 2002a, b; Haugen and Brandenburg, 2006; Chernyshov et al., 2007; Pietarila Graham et al., 2009; Sondak and Oberai, 2012), but their applicability to astrophysical plasmas is unclear. Astrophysical MHD turbulence, particularly in the interstellar medium, extends to the supersonic and super-Alfvénic regimes. Moreover, plasmas become collisionless for high temperatures and low densities. A typical example is the solar corona. It is also likely to be the case in the intracluster medium. Since the fluid-dynamical description is not applicable in this case, kinetic methods have to be employed. Nevertheless, MHD-LES could provide a reasonable approximation on length scales that are sufficiently large compared to the characteristic scales of kinetic processes. In any case, SGS models for MHD turbulence will be a very challenging problem because of the local anisotropy of turbulent fluctuations, the potentially strong back-reaction from smaller to larger scales, and complicated dissipative processes such as turbulent reconnection (Brandenburg and Subramanian, 2005; Büchner, 2007; Zweibel and Yamada, 2009). In this area, extensive fundamental studies will be necessary.

2 Scale Separation

Large eddy simulations are based on the notion of scale separation. Although turbulence is a multi-scale phenomenon, with interactions among different length scales, a separation into smoothed and fluctuating components can be rigorously defined by means of filter operators. Of course, the filtering of non-linear terms gives rise to interactions between these components. Filter operators were originally applied in the context of mean-field theories, but can be generalized to LES. For incompressible hydrodynamical turbulence, Germano (1992) introduced a general framework that encompasses mean field theories as limiting case.

The smoothed component of a generic field variable q(x,t) is defined by means of a spatial low-pass filter, which is a convolution of q with an appropriate filter kernel G (see Chapter 2 in Sagaut, 2006):

$${\langle q\rangle _G}(x) = \int {G(x - x^{\prime})q(x^{\prime},t)} {{\rm{d}}^3}x^{\prime}.$$
(4)

A homogeneous isotropic low-pass filter has the following properties:

  • The filter kernel is independent of direction:

    $$G(x - x^\prime ) = G(r),\quad {\text{where }}\,r = |x - x^\prime |.$$
  • Filtering smoothes out fluctuations on length scales smaller than the filter length Δ G . Length scales that are large in comparison to Δ G are not affected. This implies

    $$G(x - x') \sim \left\{ {\begin{array}{*{20}c} {1/\Delta _G^3 } & { if \left| {x - x'} \right| \ll \Delta G,} \\ 0 & {if \left| {x - x'} \right| \gg \Delta G.} \\ \end{array} } \right.$$
  • The filter operator is linear, conserves constants, and commutes with spatial derivatives:

    $${\langle \Delta q\rangle _G} = \nabla {\langle q\rangle _G}.$$

    The simplest low-pass filter is the box or top-hat filter. For Cartesian coordinates x i , the kernel of the box filter is defined by

    $$G_{box} (x - x') = \prod\limits_{i = 1}^3 {G_i (x_i - x'_i ), where } G_i (x_i - x'_i ) = \left\{ {\begin{array}{*{20}c} {1/\Delta _i } & {if \left| {x_i - x'_i } \right| \leqslant \Delta _i /2,} \\ 0 & {otherwise.} \\ \end{array} } \right.$$
    (5)

    Usually, Δ i is assumed to be equal for all spatial dimensions. The mean value of q in a rectangular domain with periodic boundary conditions follows in the limit that Δ i is the linear size of the domain in each dimension.

The construction of a filter is particularly simple in Fourier space. For a low-pass filter, the Fourier transform of the filter kernel, the so-called transfer function Ĝ(k), drops rapidly to zero for wavenumbers kk c = π G . Since

$${\langle \hat q\rangle _G}(k,t) = \hat q(k,t)\hat G(k),$$
(6)

only the Fourier modes \(\hat q(k,t)\) with kk c contribute significantly to the corresponding filtered field 〈q G (x) in physical space. The simplest case is the sharp cutoff filter, for which

$$\hat G_{sharp} (k) = \left\{ {\begin{array}{*{20}c} {1 if k \leqslant k_c ,} \\ {0 otherwise.} \\ \end{array} } \right.$$
(7)

The sharp cutoff filter, however, is not equivalent to the box filter, which has the Fourier representation

$${\hat G_{{\rm{box}}}}(k) = \prod\limits_{i = 1}^3 {{{\sin (k{\Delta _i}/2)} \over {k{\Delta _i}/2}}.} $$
(8)

A filter that is intermediate between these two cases is the Gaussian filter.

2.1 Decomposition of the compressible Navier-Stokes equations

The compressible Navier-Stokes equations for the mass density ρ, the momentum density ρu, and the energy density ρE of an electrically neutral fluid subject to gravitational and mechanical accelerations g and f, respectively, are

$${\partial \over {\partial t}}\rho + \nabla \cdot (u\rho) = 0,$$
(9)
$${\partial \over {\partial t}}\rho u + \nabla \cdot (\rho u \otimes u) = \rho (g + f) - \nabla P + \nabla \cdot \sigma,$$
(10)
$${\partial \over {\partial t}}\rho E + \nabla \cdot (\rho uE) = \rho u \cdot (g + f) - \nabla \cdot (uP) + \nabla \cdot (u \cdot \sigma).$$
(11)

Thermal conduction is neglected here. The energy per unit mass can be expressed as

$$E = e + {1 \over 2}{u^2},$$
(12)

where e is the internal or thermal gas energy. For a perfect gas, e is related to the gas pressure P and the temperature T via the ideal gas law:

$$e = {P \over {(\gamma - 1)\rho }} = {{{k_{\rm{B}}}T} \over {(\gamma - 1)\mu {m_{\rm{H}}}}},$$
(13)

where γ is the adiabatic exponent, kB the Boltzmann constant, μ the mean molecular weight, and mH the mass of the hydrogen atom. The viscous stress tensor σ is defined by

$${\sigma _{ij}} = 2\eta S_{ij}^{\ast} + \zeta d{\delta _{ij}} = 2\eta \left({{S_{ij}} - {1 \over 3}d{\delta _{ij}}} \right) + \zeta d{\delta _{ij}},$$
(14)

where the two coefficients η and ζ are the dynamic and bulk viscosities of the fluid,

$${S_{ij}} = {1 \over 2}\left({{{\partial {u_i}} \over {\partial {x_j}}} + {{\partial {u_j}} \over {\partial {x_i}}}} \right)$$
(15)

is the rate-of-strain tensor, and the trace S ii is equal to the divergence d = u. The gravitational acceleration is given by g = − ϕ, where the gravitational potential ϕ is determined by the Poisson equation

$${\nabla ^2}\phi = 4\pi G\rho,$$
(16)

where G is Newton’s constant.Footnote 4

Mean-field equations for compressible turbulence are derived in Canuto (1997). Much in the same way, a general low-pass filter 〈〉 G can be applied to the system of PDEs (9)(11). Alternative formulations can be found in Garnier et al. (2009), Section 2.4. For brevity, we omit the subscript G in the following. Since 〈 〉 commutes with differential operators, the smoothed mass density 〈p〉 obeys an equation of exactly the same form as the continuity equation,

$${\partial \over {\partial t}}\langle \rho \rangle + \nabla \cdot \langle \rho \rangle \tilde u = 0,$$
(17)

if we set 〈ρu〉 = 〈ρũ. This relation defines the Favre-filtered velocity

$$\tilde u = {{\langle \rho u\rangle } \over {\langle \rho \rangle }}.$$
(18)

Filtering the momentum equation results in

$${\partial \over {\partial t}}\langle \rho \rangle \tilde u + \nabla \cdot \langle \rho u \otimes u\rangle = \langle \rho (- \nabla \phi + f)\rangle - \nabla \langle P\rangle + \nabla \cdot \langle \sigma \rangle $$

Owing to the non-linearities, however, we are facing some difficulties here. To obtain a PDE with the same basic structure as the unfiltered momentum equation, the advection term on the left-hand side should read • [〈ρũũ]. The solution is to split the filtered non-linear terms:

$$\langle \rho u \otimes u\rangle = \langle \rho \rangle \tilde u \otimes \tilde u - \tau (\rho u,u)\quad {\rm{where}}\quad \tau: = - \langle \rho u \otimes u\rangle + {{\langle \rho u\rangle \otimes \langle \rho u\rangle } \over {\langle \rho \rangle }}.$$
(19)

Since the Poisson equation (16) is linear, the smoothed potential 〈ϕ〉 is solely determined by 〈ρ〉. The self-gravity term 〈ρϕ〉, however, has to be split by defining

$$\gamma: = - \langle \rho \nabla \phi \rangle + \langle \rho \rangle \nabla \langle \phi \rangle.$$
(20)

The specific external force f, on the other hand, usually varies only over the largest scales of the system (see Wagner et al., 2012). This can be any type of “stirring” force or large-scale gravitational forces acting on the system. If the filter length is small compared to these scales, 〈ρf 〉 ≃ 〈ρf is a good approximation. Thus, the filtered momentum equation can be casted into the following form (Garnier et al., 2009; Moin et al., 1991; Yoshizawa, 1991; Germano, 1992; Canuto, 1997; Schmidt et al., 2006b):

$${\partial \over {\partial t}}\langle \rho \rangle \tilde u + \nabla \cdot [\langle \rho \rangle \tilde u \otimes \tilde u] = \langle \rho \rangle (- \nabla \langle \phi \rangle + f) - \nabla \langle P\rangle + \nabla \cdot [\langle \sigma \rangle + \tau ] + \gamma.$$
(21)

Now, what is the physical interpretation of the terms τ and γ? Let us first consider the weakly compressible limit. By assuming that ρ varies only little over the filter length, density factors can be pulled out of brackets. In this case, ũ ≃ 〈u〉. By defining the fluctuation of the velocity as u = uũ, it follows that

$$\tau \simeq \rho [\langle u\rangle \otimes \langle u\rangle - \langle \langle u\rangle \otimes \langle u\rangle \rangle - 2\langle \langle u\rangle \otimes u^{\prime}\rangle - \langle u^{\prime} \otimes u^{\prime}\rangle ].$$

If we further assume that 〈 〉 is a Reynolds operator (see Section 3.3 in Sagaut, 2006), which is not generally true for filters but applies, for example, to global averages, filtered quantities can be pulled out of brackets and the above expression simplifies to

$$\tau \simeq - \rho \langle u^{\prime} \otimes u^{\prime}\rangle.$$

Although this simple relation holds only for a Reynolds operator in the weakly compressible limit, τ is generally interpreted as the stress tensor associated with the turbulent velocity fluctuations below the filter length. For this reason, τ is called the subgrid-scale turbulence stress tensor in the context of LES. The non-linear interactions of the filtered flow (the “large eddies”) with small-scale fluctuations below the grid scale ∆ are given by τ in Eq. (21). Likewise, the term γ defined by Eq. (20) accounts for the momentum transfer due self-gravitating fluctuations in the density. The trace of τ defines the fraction of kinetic energy on length scales smaller than the filter length:

$$\langle \rho \rangle K: = - {1 \over 2}{\tau _{ii}} = {1 \over 2}\langle \rho {u^2}\rangle - {1 \over 2}\langle \rho \rangle {\tilde u^2}.$$
(22)

If the filter length is the grid scale, ρK is called the subgrid-scale turbulence energy. The first term on the right-hand side of Eq. (22) is the total kinetic energy, the second term the kinetic energy on length scales greater than the filter length (i.e., the numerically resolved kinetic energy in LES).

In the limit of high Reynolds numbers, the viscous dissipation scale (also known as Kolmogorov scale) is typically much smaller than the filter length. In this case, scaling arguments for incompressible turbulence imply 〈σ〉 ≪ τ, i.e., the filtered viscous stresses are negligible compared to the stresses associated with the turbulent velocity fluctuations (Röpke and Schmidt, 2009). Since the scaling of compressible turbulence tends to be stiffer than for incompressible turbulence (Kritsuk et al., 2007; Schmidt et al., 2008, 2009), one can reasonably assume that this conclusion is generally applicable. The filtered momentum equation (21) thus can be written as

$${\partial \over {\partial t}}\langle \rho \rangle \tilde u + \nabla \cdot [\langle \rho \rangle \tilde u \otimes \tilde u] = \langle \rho \rangle (- \nabla \langle \phi \rangle + f) - \nabla \left({\langle P\rangle + {2 \over 3}\rho K} \right) + \nabla \cdot {\tau ^{\ast}} + \gamma.$$
(23)

where τ* is the trace-free part of τ:

$$\tau _{ij}^{\ast} = {\tau _{ij}} - {1 \over 3}{\tau _{kk}}{\delta _{ij}} = {\tau _{ij}} + {2 \over 3}\rho K{\delta _{ij}}.$$
(24)

As one can see from Eq. (23), the trace of τ is associated with the turbulent pressure 2/3ρK at the filter length scale.

In contrast to the filtered momentum density, which can be expressed as 〈ρu〉 = 〈ρũ, the energy density on length scales greater than the filter length is given by

$$\langle \rho \rangle \tilde E: = \langle \rho \rangle \left({\tilde e + {1 \over 2}{{\tilde u}^2}} \right) = \langle \rho E\rangle - \langle \rho \rangle K,$$
(25)

where the second equality follows from Eqs. (12) and (22). Consequently, 〈ρ ≠ 〈ρE〉.Footnote 5 A PDE for 〈ρ follows from the contraction of Eq. (23) with ũ plus the filtered internal energy equation. The subtraction of this PDE from the filtered equation for the total energy yields the PDE for ρK (see Section 3.3 in Sagaut, 2006 and Moin et al., 1991; Yoshizawa, 1991; Germano, 1992; Canuto, 1997; Schmidt et al., 2006b). In the limit of high Reynolds numbers, the resulting equations are:

$$\begin{array}{*{20}c} {\frac{\partial } {{\partial t}}\left\langle \rho \right\rangle \tilde E + \nabla \cdot \left\langle \rho \right\rangle \tilde u\tilde E = \left\langle \rho \right\rangle \tilde u \cdot \left( { - \nabla \left\langle \varphi \right\rangle + f} \right) + \nabla \cdot \left[ { - \tilde u\left( {\left\langle P \right\rangle + \frac{2} {3}\rho K} \right) + \tilde u \cdot \tau * + \mathfrak{F}(conv)} \right]} \\ { - \Sigma + \left\langle \rho \right\rangle (\varepsilon + \lambda ) + \tilde u \cdot \gamma ,} \\ \end{array} $$
(26)
$${\partial \over {\partial t}}\langle \rho \rangle K + \nabla \cdot \langle \rho \rangle \tilde uK = \Gamma + \sum - \langle \rho \rangle (\epsilon + \lambda) + \nabla \cdot \left[ {{\mathfrak{F}^{({\rm{kin}})}} + {\mathfrak{F}^{({\rm{press}})}}} \right].$$
(27)

The additional source and transport terms resulting form the scale separation of the energy are defined as follows.

  • Gravitational energy injection on subgrid scales:

    $$\Gamma = - \langle \rho \tilde u \cdot \nabla \phi \rangle + \tilde u \cdot \langle \rho \nabla \phi \rangle = - \langle \rho \tilde u \cdot \nabla \phi \rangle + \langle \rho \rangle \tilde u \cdot \nabla \langle \phi \rangle - \tilde u \cdot \gamma.$$
    (28)
  • Rate of subgrid-scale turbulence energy production:Footnote 6

    $$\Sigma = {\tau _{ij}}{\tilde S_{ij}},$$
    (29)

    where is defined by Eqs. (19) and \({\tilde S_{ij}}\) is the rate-of-strain tensor associated with the Favre-filtered velocity:Footnote 7

    $${\tilde S_{ij}}: = \frac{1}{2}\left( {\frac{{\partial \tilde u}}{{\partial {x_j}}}\frac{{\partial {{\tilde u}_j}}}{{\partial {x_i}}}} \right).$$
    (30)
  • Rate of viscous energy dissipation in the limit of high Reynolds numbers:Footnote 8

    $$\langle \rho \rangle \epsilon = \langle {\sigma _{ij}}{S_{ij}}\rangle - \langle {\sigma _{ij}}\rangle {\tilde S_{ij}} \simeq \langle {\sigma _{ij}}{S_{ij}}\rangle = \langle \eta \vert{S^{\ast}}{\vert^2} + \zeta {d^2}\rangle,$$
    (31)

    where is defined by Eq. (15), \(\vert{S^{\ast}}{\vert^2} = 2S_{ij}^{\ast}S_{ij}^{\ast}\) is the squared norm of the trace-free rate-of-strain tensor \(S_{ij}^{\ast} = {S_{ij}} - {1 \over 3}d{\delta _{ij}}\) and d = S ii . Although the viscous stresses can be neglected in the filtered momentum equation, viscous dissipation is crucial for the energy balance of turbulent flows.

  • Rate of subgrid-scale pressure dilatation:

    $$\langle \rho \rangle \lambda = - \langle dP\rangle + \tilde d\langle P\rangle,$$
    (32)

    where \(\tilde d = {\tilde S_{ii}} = \partial {\tilde u_i}/\partial {x_i}\).

  • Convective internal energy flux on sub-grid scales:Footnote 9

    $${\mathfrak{F}^{({\rm{conv}})}} = - \langle \rho ue\rangle + \tilde \rho \tilde u\tilde e.$$
    (33)
  • The flux associated with pressure fluctuations:

    $${\mathfrak{F}^{({\rm{press}})}} = - \langle uP\rangle + \tilde u\langle P\rangle.$$
    (34)

    For ideal gas with adiabatic exponent \(\gamma,\;{\mathfrak{F}^{({\rm{press}})}} = (\gamma - 1){\mathfrak{F}^{({\rm{conv}})}}\).

  • The diffusive flux of turbulent energy on sub-grid scales:

    $${\mathfrak{F}^{({\rm{kin}})}} = - {1 \over 2}\langle \rho {u^2}u\rangle + {1 \over 2}\langle \rho {u^2}\rangle \tilde u - \tilde u \cdot \tau.$$
    (35)

    For Reynolds operators in the weakly compressible limit, \({\mathfrak{F}^{({\rm{kin}})}}\) can be expressed as a third-order moment of the velocity fluctuation: \(2\mathfrak{F}_j^{({\rm{kin}})} \simeq - \rho \langle u_i^{\prime}u_i^{\prime}u_j^{\prime}\rangle \) (Germano, 1992).

  • There is also a viscous flux, which can be neglected relative to other flux terms if the Reynolds number is sufficiently high.

By adding the Eqs. (26) and (27), we obtain an equation for the filtered total energy

$$\begin{array}{*{20}c} {\frac{\partial } {{\partial t}}\left\langle \rho \right\rangle (\tilde E + K) + \nabla \cdot \left\langle \rho \right\rangle \tilde u(\tilde E + K) = \left\langle \rho \right\rangle \tilde u\left( { - \nabla \left\langle \varphi \right\rangle + f} \right) + \Gamma + \tilde u \cdot \gamma } \\ { + \nabla \cdot \left[ { - \tilde u\left( {\left\langle P \right\rangle + \frac{2} {3}\rho K} \right) + \tilde u \cdot \tau * + \mathfrak{F}^{(kin)} + \mathfrak{F}^{(conv)} + \mathfrak{F}^{(press)} } \right].} \\ \end{array} $$
(36)

Except for the source terms related to self-gravity and external forces, production and dissipation rates cancel out. The rightmost fluxes are related to turbulent transport processes below the filter length. In particular, the sum of \({\mathfrak{F}^{({\rm{conv}})}}\) and \({\mathfrak{F}^{({\rm{press}})}}\) can be expressed as convective enthalpy flux,

$$\mathfrak{F}^{(conv)} + \mathfrak{F}^{(conv)} = - \left\langle {\rho uh} \right\rangle + \rho \tilde u\tilde h,$$
(37)

where ρh = ρe + P, which corresponds to −ρuh′〉 in the weakly compressible limit. For a closed system of PDEs, it is necessary to compute all terms defined above in terms of known quantities. A rigorous calculation requires further PDEs, which involve higher-order moments and so on ad infinitum. This is known as the closure problem. A subgrid-scale model truncates the closure problem by approximating moments above a given order by lower-order moments.

2.2 Cosmological fluid dynamics

In cosmological simulations, the equations of fluid dynamics are solved in a comoving coordinate system. Coordinates of observers that are stationary relative to the Hubble expansion of the Universe are constant in this system. The expansion is characterized by the scale factor a(t), which is determined by the Friedmann equations for a homogeneous and isotropic cosmology (Peacock, 1999). If the proper coordinates, which include changes of position due to the expansion of the Universe, are denoted by xproper and tproper, the corresponding comoving coordinates are x = xproper/a and t = tproper. Derivative operators transform as

$${\left. {{\partial \over {\partial t}}} \right\vert_{{\rm{proper}}}} = {\partial \over {\partial t}} - {{\dot a} \over a}x \cdot \nabla \quad {\rm{and}}\quad {\nabla _{{\rm{proper}}}} = {1 \over a}\nabla.$$

Furthermore, the invariance of mass implies that the comoving baryonic density ρ is related to the proper density by ρ = a3ρproper. It can then be shown that the continuity equation for ρ in comoving coordinates assumes exactly the same form as Eq. (9):

$${{\partial \rho } \over {\partial t}} + \nabla \cdot (\rho u) = 0.$$

Here, u is the so-called peculiar velocity, which is defined as

$$u = \dot x = {1 \over a}{u_{{\rm{proper}}}} - Hx,$$
(38)

where uproper = proper is the proper velocity and H = ȧ/a the Hubble constant. This means that, in the comoving coordinate system, matter moves with velocity u relative to the Hubble flow Hx. With some algebra, also the momentum and energy equations can be transformed to comoving coordinates. The resulting equations do not have the same form as Eqs. (10) and (11), but include additional terms with prefactors H. However, a particularly simple representation of the momentum and energy equations is obtained if the proper peculiar velocity

$$U = au = {u_{{\rm{proper}}}} - \dot ax$$
(39)

is used in place of u.

Filtered dynamical equations for cosmological fluids were first derived in Maier (2008) and Maier et al. (2009), and presented in an alternative formulation in Schmidt et al. (2014). The applied filter kernel is static in comoving coordinates, i.e., the filter length increases proportional to the cosmological scale factor a. Consequently, commutation of the filter with time derivatives is unaffected by the cosmological expansion and equations for filtered dynamical variables follow completely analogous to Section 2.1. By neglecting gravitational terms associated with fluctuations below the filter length, the following equations for the filtered mass density 〈ρ〉, the filtered momentum density 〈ρU〉 = 〈ρŨ, and energy density 〈ρ, where = + ½Ũ2, are obtained:

$${{\partial \langle \rho \rangle } \over {\partial t}} + {1 \over a}\nabla \cdot [\langle \rho \rangle \tilde U] = 0,$$
(40)
$${\partial \over {\partial t}}a\langle \rho \rangle \tilde U + \nabla \cdot [\langle \rho \rangle \tilde U \otimes \tilde U] = - \langle \rho \rangle \nabla \langle \phi \rangle - \nabla \langle P\rangle + \nabla \cdot \tau + \gamma,$$
(41)
$$\begin{array}{*{20}c} {\frac{\partial } {{\partial t}}a^2 \left\langle \rho \right\rangle \tilde E + a\nabla \cdot \left[ {\left\langle \rho \right\rangle \tilde U\tilde E} \right] = - a\left\langle \rho \right\rangle \tilde U \cdot \nabla \left\langle \varphi \right\rangle + a\nabla \left[ { - \tilde U\left\langle P \right\rangle + \tilde U \cdot \tau + \mathfrak{F}^{(conv)} } \right]} \\ { - a\dot a\left[ {2 - 3\left( {\gamma - 1} \right)} \right]\left\langle \rho \right\rangle \tilde e - a\left[ {\Sigma + \left\langle \rho \right\rangle \left( {\varepsilon + \lambda } \right)} \right] + a\tilde U \cdot \lambda ,} \\ \end{array} $$
(42)

Here, the filtered internal energy density is 〈ρ = 〈ρe〉 = 〈P〉/ (γ − 1), where P = a3Pproper, and the gravitational potential 〈ϕ〉 of baryonic and dark matter density fluctuations is given by the cosmological Poisson equation

$${\nabla ^2}\langle \phi \rangle = {{3H_0^2{\Omega _{\rm{m}}}({t_0})} \over {2a}}\langle {\delta _{\rm{m}}}\rangle,$$
(43)

where H0 = ȧ(t0) is the Hubble constant and Ωm(t0) the density parameter of matter at redshift zero. Since the mean matter density ρm,0 is constant in comoving coordinates, the source term of the Poisson equation can expressed in terms of the density fluctuation δm = (〈ρdm + ρ〉 − ρm,0)/ρm,0 for the local dark matter density ρdm and baryonic mass density ρ. The density parameter is defined by Ωm(t0) = ρm,0/ρcrit,0, where \({\rho _{{\rm{crit}},0}} = 3H_0^2/(8\pi G)\) is the critical density at time t = t0.

The kinetic energy associated with peculiar velocity fluctuations below the filter length,

$$\langle \rho \rangle K = {1 \over 2}\langle \rho {U^2}\rangle - {1 \over 2}\langle \rho \rangle {\tilde U^2},$$

is given by the dynamical equation

$${\partial \over {\partial t}}{a^2}\langle \rho \rangle K + a\nabla \cdot [\langle \rho \rangle \tilde UK] = a[\Gamma + \Sigma - \langle \rho \rangle (\epsilon + \lambda)] + a\nabla \cdot \left[ {{\mathfrak{F}^{({\rm{kin}})}} + {\mathfrak{F}^{({\rm{press}})}}} \right].$$
(44)

Since the prefactors of all terms except for the time derivatives in the momentum and energy equations are unity and a, respectively, the definitions of all source and transport terms in Eqs. (41), (42) and (44) are analogous to the definitions given in Section 2.1, with u i being replaced by U i . This suggests that closures for turbulence in a static space are applicable to cosmological fluids as well. Although cosmological expansion, in principle, causes a dampening of the kinetic energy (Schmidt et al., 2014), this effect is subdominant for turbulent eddies even on the largest scales in galaxy clusters because turbulence is driven in gravitationally bound gas on time scales shorter than the current Hubble time 1/H0.

3 Subgrid-Scale Models

There is a beautiful correspondence between finite-volume discretization and filtering. Finite-volume methods solve an equation for the cell averages of some dynamical variable q(x):

$${Q_{ijk}} = \int_{{z_k} - \Delta/2}^{{z_k} + \Delta/2} {\int_{{y_j} - \Delta/2}^{{y_j} + \Delta/2} {\int_{{x_i} - \Delta/2}^{{x_i} + \Delta/2} {q(x,y,z){\rm{d}}x\,{\rm{d}}y\,{\rm{d}}z} } } $$

Here, (x i , y j , z k ) are the cell-centered coordinates and Δ is the linear size of a grid cell. It is not difficult to see that Q ijk equals the box-filtered variable 〈q G for a box filter G with filter length Δ at the discrete points (x i , y j , z k ).Footnote 10 Also finite differences correspond to low-pass filters. Thus, numerical discretization can be interpreted as implicit filtering. The numerical errors in the approximations to Q ijk can be characterized by the truncation error of the finite-volume method. For stable schemes with flux limiters, these errors are usually associated with terms of the diffusion type, i.e., proportional to ∇2q.Footnote 11 In ILES, this is what causes the dissipation of kinetic energy into heat (a detailed account of ILES is given in Garnier et al., 2009). A similar expression follows if the Boussinesq expression for the turbulent stresses is used as an explicit SGS model for the interaction between numerically resolved and unresolved turbulent eddies. An important difference, however, is that the turbulent viscosity in the resulting diffusion terms in the momentum and energy equations is controlled by the dynamical variable K, which is the kinetic energy associated with turbulent velocity fluctuations below the grid scale. In this section, we mainly discuss the computation of K in LES. The approach we follow here is known as functional modeling. The aim is to model only statistical effects of SGS turbulence on the dynamics of the filtered fields, which are identified with the numerical solution. An alternative strategy is structural modeling (see Chapter 5 in Garnier et al. 2009), which is not covered here.

3.1 Closures for the turbulence stress tensor

The SGS turbulence stress tensor, which is associated with the non-linear energy transfer between large and small scales, is the central quantity that has to be modeled in LES. The most commonly used closure is the eddy-viscosity closure. The underlying assumption is that the form of the trace-free part τ* is analogous to the anisotropic viscous stress tensor σ*, with the correspondenceFootnote 12

$${S_{ij}} \leftrightarrow {\tilde S_{ij}},$$
(45)
$$\eta \leftrightarrow \langle \rho \rangle {\nu _{{\rm{sgs}}}},$$
(46)

where S ij and \({\tilde S_{ij}}\) are defined by Eqs. (15) and (30), respectively. The turbulent viscosity νsgs is assumed to depend on the grid scale and the unresolved turbulent velocity fluctuation (see, for example, Section 4.3 in Sagaut, 2006):

$${\nu _{{\rm{sgs}}}} = {C_\nu }\Delta \sqrt K.$$
(47)

Hence,

$$\tau _{ij}^{({\rm{eddy}})} = 2\langle \rho \rangle \left({{\nu _{{\rm{sgs}}}}\tilde S_{ij}^{\ast} - {1 \over 3}K{\delta _{ij}}} \right).$$
(48)

For brevity, we drop brackets and tildes indicating filtered and Favre-filtered quantities from now onwards, so that the turbulent stresses can be written as

$$\tau _{ij}^{({\rm{eddy}})} = 2\rho \left({{\nu _{{\rm{sgs}}}}S_{ij}^{\ast} - {1 \over 3}K{\delta _{ij}}} \right).$$
(49)

In the following, it is understood that all quantities are either numerically resolved variables or modeled in terms of these variables. The production rate (energy flux) corresponding to the eddy-viscosity closure is

$${\Sigma ^{({\rm{eddy}})}} = {C_\nu }\rho \Delta {K^{1/2}}\vert{S^{\ast}}{\vert^2} - {2 \over 3}\rho Kd,$$
(50)

where d = S ii is the divergence. The eddy-viscosity coefficient C ν is typically in the range from 0.05 and 0.1 (Sagaut, 2006; Schmidt et al., 2006b; Schmidt and Federrath, 2011).

In the incompressible case (d = 0), the eddy-viscosity closure admits only positive energy flux. However, direct numerical simulation data show that there is a certain amount of backscattering from smaller to larger scales, corresponding to a negative energy flux (Schmidt et al., 2006b; Schmidt and Federrath, 2011). This motivated a closure for the turbulent viscosity that is constructed from the determinant of the trace-free rate-of-strain tensor (Woodward et al., 2001):Footnote 13

$${\nu _{{\rm{sgs}}}} = - {{{C_1}{\Delta ^2}\det {{\bf{S}}^{\ast}}} \over {\vert{S^{\ast}}{\vert^2}}}.$$
(51)

By substituting the above expression into Eq. (49) for the turbulent stresses, it follows that the production rate is given by

$${\Sigma ^{(\det)}} = - {C_1}\rho {\Delta ^2}\det {{\bf{S}}^{\ast}} - {2 \over 3}\rho Kd,$$
(52)

Since the determinant can be positive under certain flow conditions, in principle, this closure accounts for backscattering (also known as inverse cascade). This phenomenon can be explained by the so-called the “tornado” topology, i.e., the alignment of vortices along a single stretching direction (Woodward et al., 2001). Then the flow is contracting in one dimension and expanding in the other two, which results in a positive determinant. A negative determinant, on the other hand, corresponds to the standard situation of a forward cascade transporting energy from larger to smaller eddies.

While the determinant closure modifies only the turbulent viscosity, a different expression for the SGS turbulence stress tensor is proposed for compressible turbulence in Woodward et al. (2006). Based on Taylor series expansions of the velocity around grid cell centers, an appropriate normalization leads to the non-linear closure

$$\tau _{ij}^{({\rm{nonlin}})} = 4\rho K{{{u_{i,k}}{u_{j,k}}} \over {{{\vert {\nabla \otimes u} \vert}^2}}},$$
(53)

where ∣u∣ = (2u i,k u i,k )1/2 is the norm of the velocity derivative u.Footnote 14 The above expression satisfies the identity τ ii = −2ρK. However, it is generally not adequate as a model for the turbulence stress tensor in LES (Schmidt and Federrath, 2011). In contrast to the eddy-viscosity closure, rotation invariance is violated because of the antisymmetric part of u. This would cause spurious production of K in a uniformly rotating fluid. A further problem is that K = 0 would be a fixed point of Eq. (27) if all other sources of turbulence energy are zero. This results in unphysical behavior. With the eddy-viscosity closure, on the other hand, K can grow sufficiently fast from arbitrarily small initial values because νsgs is proportional to √K rather than K. For this reason, a linear combination of \(\tau _{ij}^{({\rm{nonlin}})}\) and \(2{\nu _{{\rm{sgs}}}}S_{ij}^{\ast}\) is used in Woodward et al. (2001), where νsgs is given by Eq. (51). The additional determinant term, with a small coefficient C1, has the function of a seed term that triggers the production of turbulence energy, while the production rate vanishes for a uniformly rotating fluid.

With the standard turbulent viscosity defined by Eq. (47), the same idea leads to the following generalized two-coefficient closure (Schmidt and Federrath, 2011):

$${\tau _{ij}} = 2\rho \left[ {{C_1}\Delta {{(2K)}^{1/2}}S_{ij}^{\ast} - 2{C_2}\rho K{{{u_{i,k}}{u_{j,k}}} \over {{{\vert {\nabla \otimes u} \vert}^2}}} - {1 \over 3}(1 - {C_2})K{\delta _{ij}}} \right].$$
(54)

The coefficient C2 determines the relative contributions from the non-linear and divergence terms to the trace τ ii . The purely non-linear closure corresponds to C1 =0 and C2 = 1. Equation (49), on the other hand, is obtained if C1 = C ν /√2 and C2 = 0. For the application in LES, it is necessary to calibrate the closure coefficients C1 and C2. For supersonic turbulence, C1 = 0.02 and C2 = 0.7 appear to be robust values (see Section 4). The rate of production following from the generalized closure is

$$\Sigma = {C_1}\rho \Delta {(2K)^{1/2}}\vert{S^{\ast}}{\vert^2} - 4{C_2}\rho K{{{u_{i,k}}{u_{j,k}}S_{ij}^{\ast}} \over {{{\left\vert {\nabla \otimes u} \right\vert}^2}}} - {2 \over 3}\rho Kd.$$
(55)

The first term dominates if K1/2 is small compared to Δ∣S*∣. For strong turbulence intensity, i. e., K1/2 ≳ Δ∣u∣, the second term contributes significantly. The transition is further influenced by the ratio C2/C1.

3.2 The Sarkar-Smagorinsky model for weakly compressible turbulence

In the case of isotropic incompressible turbulence, the mean SGS turbulence energy \(\overline K \) for a sharp cutoff at the length scale Δ is obtained by integrating the Kolmogorov spectrum E(k) over wavenumbers kπ/Δ:

$$\overline K = \int_{\pi/\Delta } {E(k){\rm{d}}k = {3 \over 2}C{{\bar \epsilon }^{2/3}}{{\left({{\pi \over \Delta }} \right)}^{ - 2/3}}.} $$
(56)

The mean dissipation rate is therefore given by

$$\bar \epsilon = \pi {\left({{{3C} \over 2}} \right)^{ - 3/2}}{{{{\overline K }^{3/2}}} \over \Delta } \approx 0.81{{{{\overline K }^{3/2}}} \over \Delta }$$
(57)

for the Kolmogorov constant C ≈ 1.65 (Pope, 2000). It is commonly assumed that an expression of this form also holds for the local dissipation rate in LES (see, for example, Sagaut, 2006):

$$\epsilon = {C_\epsilon }{{{K^{3/2}}} \over \Delta },$$
(58)

with C ε ∼ 1. The above dimensional closure for the dissipation rate basically means that the time scale of energy dissipation is given by τ ε ∼ Δ/√K.

For subsonic compressible turbulence, closures for the dissipation rate and pressure dilatation are obtained by separating the pressure fluctuations into a rapid oscillatory and a slow component (Sarkar, 1992). The resulting combined expression for ε + λ reads

$$\rho (\epsilon + \lambda) = {C_\epsilon }[1 + ({\alpha _1} - {\alpha _3})\mathcal{M}_{{\rm{sgs}}}^2]{{\rho {K^{3/2}}} \over \Delta } + {\alpha _2}{\mathcal{M}_{{\rm{sgs}}}}\tau _{ij}^{\ast}S_{ij}^{\ast} - {{16} \over 3}{\alpha _4}\mathcal{M}_{{\rm{sgs}}}^2\rho Kd,$$
(59)

where

$${\mathcal{M}_{{\rm{sgs}}}} = {{\sqrt {2K} } \over {{c_{\rm{s}}}}}$$
(60)

is the turbulent Mach number associated with the SGS velocity fluctuation √2K.

A particularly simple SGS model can be formulated by neglecting all gravitational and transport terms associated with subgrid-scale effects in Eq. (27). If furthermore a balance between production and dissipation is assumed, then

$$\Sigma \simeq \rho (\epsilon + \lambda)$$

implies

$${C_\nu }(1 - {\alpha _2}{\mathcal{M}_{{\rm{sgs}}}})\Delta {K^{1/2}}\vert {S^{\ast}}{\vert ^2} - {2 \over 3}(1 - 8{\alpha _4}\mathcal{M}_{{\rm{sgs}}}^2)Kd \simeq {C_\epsilon }[1 + ({\alpha _1} - {\alpha _3})\mathcal{M}_{{\rm{sgs}}}^2]{{{K^{3/2}}} \over \Delta }.$$
(61)

Here, the eddy-viscosity closure (49) is substituted for \(\tau _{ij}^{\ast}\). If the above algebraic equation is solved for K, the PDEs (17), (23), and (26) form a closed system. The effect of the compressibility corrections is a reduction of the production due to anisotropic shear, νsgsS∣*, by the factor \((1 - {\alpha _2}{\mathcal{M}_{{\rm{sgs}}}})\) and an enhancement of the solenoidal dissipation rate C ν K3/2/Δ by a factor that increases with the square of \({\mathcal{M}_{{\rm{sgs}}}}\) (α1 tends to be greater than α3). An extension to a non-equilibrium model based on the dynamical equation (27) for K was exploited in Maier (2008) and Maier et al. (2009) for cosmological LES of the gas in galaxy clusters. However, this model applies only if \({\mathcal{M}_{{\rm{sgs}}}}\) is small compared to unity, which is the case for turbulence in the intracluster medium. On the other hand, the correction factors are close to unity for small \({\mathcal{M}_{{\rm{sgs}}}}\) and, given the many approximations involved, it is not clear whether they have any significant effect. Apart from that, the model definitely breaks down in the vicinity of accretion shocks and in the cooler regions of the intergalactic medium, where \({\mathcal{M}_{{\rm{sgs}}}}\) can become large compared to unity.

In the limit \({\mathcal{M}_{{\rm{sgs}}}} \rightarrow 0\) and d → 0, the classical Smagorinsky model for incompressible turbulence (Smagorinsky, 1963) follows from Eq. (61). In this case,

$$K \simeq {{{C_\nu }} \over {{C_\epsilon }}}{\Delta ^2}\vert S{\vert ^2}\quad {\rm{and}}\quad {\nu _{{\rm{sgs}}}} = {({C_s}\Delta)^2}\vert S\vert,$$
(62)

where \({C_{\rm{s}}} = {(C_\nu ^3/{C_\epsilon })^{1/4}}\). The corresponding equilibrium dissipation rate is

$$\epsilon \simeq {({C_{\rm{s}}}\Delta)^2}\vert S{\vert ^3}.$$
(63)

This expression has an important implication. One could calculate from ILES data analogous to the viscous dissipation rate following from the Navier-Stokes equations, i.e.,

$$\epsilon \sim {\nu _{{\rm{eff}}}}\vert S{\vert ^2}.$$

Here, νeff = VL/Reeff is assumed to be the constant numerical viscosity, which is given by the effective Reynolds number Reeff of the simulation (Pan et al., 2009). However, the above estimation of the dissipation rate is clearly at odds with Eq. (63). This is not only an inaccuracy, but a genuine inconsistency. If Reeff ∼ Re, where Re is the physical Reynolds number defined by Eq. (1), then νS2 is the physical dissipation rate in a direct numerical simulation of a fluid with microscopic viscosity ν. In an LES with Reeff ≪ Re, on the other hand, the Smagorinsky model implies a turbulent viscosity of the order Δ2S∣ for steady-state turbulence, which is not a constant. In this case, the dissipation rate is approximately given by Eq. (63). This is a consequence of Eq. (31), which implies that ρε cannot be expressed as the contraction of the filtered viscous stress tensor with the filtered rate-of-strain tensor. The dissipation rate is instead given by the filtered contraction of the two tensors. As shown in Schmidt and Federrath (2011), the argument remains valid even if the dissipation rate is calculated for LES of supersonic turbulence with the advanced SGS model presented in the following section.

3.3 The compressible subgrid-scale turbulence energy model

To determine K, one can either invoke the equilibrium condition, such as in the Smagorinsky model, or numerically solve the PDE (27). The latter is called the SGS turbulence energy model (Sagaut, 2006; Schumann, 1975; Moin et al., 1991; Yoshizawa, 1991; Germano, 1992; Schmidt et al., 2006b; Schmidt and Federrath, 2011). If gravitational terms are negligible, the turbulence energy equation can be explicitly written as

$$\begin{array}{*{20}c} {\frac{\partial } {{\partial t}}\rho K + \nabla \cdot (\rho uK) = C_1 \rho \Delta (2K)^{1/2} \left| {S*} \right|^2 - 4C_2 \rho K\frac{{u_{i,k} u_{j,k} S_{ij}^* }} {{\left| {\nabla \otimes u} \right|^2 }} - \frac{2} {3}\rho Kd} \\ { - C_\varepsilon \frac{{K^{3/2} }} {\Delta } + \nabla \cdot \left[ {C_\kappa \rho \Delta K^{1/2} \nabla K} \right],} \\ \end{array} $$
(64)

where the closure (55) for the production rate Σ, the dissipation rate ε defined by Eq. (58), and the gradient-diffusion closure for \({\mathfrak{F}^{({\rm{kin}})}} + {\mathfrak{F}^{({\rm{press}})}}\) were substituted into Eq. (27). The gradient-diffusion hypothesis, which is also known as Kolmogorov-Prandtl relation, is based on the assumption that the turbulent transport of K is a diffusion process satisfying Fick’s law (see, for example, Pope, 2000; Schmidt et al., 2006b):

$${\Im ^{({\rm{kin}})}} + {\Im ^{({\rm{press}})}} = \rho {\kappa _{{\rm{sgs}}}}\nabla K,$$
(65)

with a turbulent diffusivity

$${\kappa _{{\rm{sgs}}}} = {C_\kappa }\Delta \sqrt K = {{{C_\kappa }} \over {{C_\nu }}}{\nu _{{\rm{sgs}}}}.$$
(66)

The Prandtl number of turbulent transport, C κ /C ν , is often assumed to be of the order unity. For a calibration of C κ , see Section 4. The pressure-dilatation λ is assumed to be negligible in Eq. (64) because no satisfactory closure is known for the highly compressible regime (Woodward et al., 2006; Schmidt and Federrath, 2011). For weakly compressible turbulence, Eq. (59) could be used, but the applicability of this closure and an appropriate calibration of the coefficients α1,…,α4 requires further investigation. As pointed out in Section 3.2, the contribution from λ is too small to significantly influence K in the weakly compressible regime. In this case, the equation with the eddy-viscosity closure, i. e., C1 = C ν /√2 and C2 = 0, can be regarded as a sufficient model for most applications. In particular, this variant of the SGS turbulence energy model was used for simulations of thermonuclear combustion in white dwarfs (Section 6.1) and cosmological simulations (Section 6.3).

Negligible pressure-dilatation is also a reasonable assumption at high Mach numbers because the kinetic energy is large compared to the internal energy and non-linear interactions between turbulent velocity fluctuations should be the dominant mode of energy transfer. This is manifest in the closure (50) for Σ, which is solely constructed from the gradient of the resolved velocity field, but does not depend on density or pressure gradients. There are both theoretical and numerical studies in support of this conjecture. In Aluie (2011, 2013), it is argued that a range of length scales exists, in which the kinetic and internal energies decouple and the flux through the kinetic energy cascade becomes asymptotically constant, while ρλ is subdominant. The computation of the different contributions to the total energy flux at varying length scales from supersonic turbulence data in Kritsuk et al. (2013) confirms this conclusion. Based on an analytical theory for the two-point correlations of compressible turbulence (Galtier and Banerjee, 2011), it is shown that the main contribution to the energy flux is

$${F_{\Vert}}(r) = \langle \delta (\rho u) \cdot \delta u\;\delta {u_{\Vert}}\rangle,$$
(67)

where δu is the velocity difference between two points separated by a distance r, δu is the longitudinal component of the velocity difference in the direction of r, and the brackets denote the ensemble average. The closure (55) for the Σ has a similar structure, with factors of ρK1/2 and derivatives of u corresponding to fluctuations on the grid scale.

However, as pointed out in Section 2.5.2 of Garnier et al. (2009), a subtlety arises in the presence of shocks because the Rankine-Hugeniot conditions for jumps across shock fronts should be filtered in place of the PDEs. This entails SGS terms that are different from the terms in the filtered PDEs. However, it is questionable whether any attempt to model these terms would be useful. The assumption that the unmodified jump conditions apply to the numerical solution amounts to a fallback from LES with an explicit SGS model to ILES. Since shock-capturing schemes, such as PPM, fall back to stronger diffusion in the vicinity of shocks, this is probably the most reasonable thing one can do. Nevertheless, the closure (55) captures the non-linear interscale transfer of energy due to supersonic turbulent velocity fluctuations. The SGS model outlined above accounts for the statistical effect of shocks as well as vortices interacting with each other across the grid scale (Schmidt et al., 2008, 2009), while any SGS terms in the jump conditions would mainly correct geometric differences between smoothed shock fronts and the corresponding unfiltered fronts with substructure on smaller length scales (just like the turbulent flame fronts discussed in Section 6.1).

Figure 1:
figure 1

Visualization of the SGS turbulence energy density ρK in a 5123 LES with solenoidal forcing. Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

Schmidt and Federrath (2011) demonstrated that Eq. (64) for ρK works very well in the highly compressible regime. As an example, Figure 1 shows a visualization of ρK from an LES of isotropic supersonic turbulence, where solenoidal stochastic forcing maintains a root mean square Mach number of about 5 in the statistically stationary regime.Footnote 15 The numerical resolution is 5123. In the reddish regions of the plot, Ksgs is higher than the spatial average, while it is lower in the bluish regions. The structure of the numerically resolved turbulent flow is illustrated by the so-called denstrophy,

$${\Omega _{1/2}} = {1 \over 2}{\left\vert{\nabla \times \left({{\rho ^{1/2}}u} \right)} \right\vert^2},$$

in Figure 2. Since Ω1/2 combines density fluctuations and the rotation of the velocity, ∇ × u, it indicates both small-scale compression and eddy-like motion (Kritsuk et al., 2007). There is clearly a correlation between Ω1/2 and ρK, which reflects the local interaction between resolved small-scale modes and subgrid-scale turbulence, as expressed by the production terms in Eq. (64). This correlation is akin to the equilibrium condition (62) following from the Smagorinsky model for incompressible turbulence. Owing to the non-local effects in the PDE (64), however, the SGS turbulence energy cannot be reliably estimated from local quantities such as Ω1/2 (Schmidt and Federrath, 2011). In particular, turbulent diffusion smears out ρK in comparison to Ω1/2.

Figure 2:
figure 2

Visualization of the the denstrophy Ω1/2 for the same LES as in Figure 1. Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

A critical property is the scaling behavior of the SGS turbulence energy. For statistically stationary homogeneous turbulence, the mean value of ρK should scale as a power of the grid resolution because the fraction of unresolved kinetic energy changes as the the cutoff of the energy spectrum is shifted (see Eq. 56). This was verified in Schmidt and Federrath (2011) by running LES with different grid scales Δ and fixed forcing length L. The global spatial averages 〈ρK〉 in these simulations are plotted in Figure 3 (left panel) for Δ ranging from L/256 to L/32, where the case Δ = L/256 corresponds to the 5123 simulation depicted in Figures 1 and 2. Although there are substantial fluctuations, one can qualitatively see that 〈ρK〉 decreases with Δ. Time-averaging over the statistically stationary regime yields mean values that are close to the power law

$$\langle \rho K\rangle \propto {\Delta ^\alpha },$$
(68)

with α ≈ 0.799 ± 0.009 (power-law fits are shown in Figure 4). The scaling exponent is in between the Kolmogorov and Burgers exponents and roughly comparable to the slope of the second-order structure functions with fractional mass-weighing reported in Schmidt et al. (2008).

In the filtered momentum and energy equations (21) and (26), respectively, the trace of the SGS turbulence stress tensor acts as an additional turbulent pressure. The sum of the thermal and turbulent pressures is sometimes called the effective pressure:

$${P_{{\rm{eff}}}} = P + {2 \over 3}\rho K.$$
(69)

It is important to keep in mind that Peff depends on the numerical resolution and PeffP in the limit Δ → 0 (DNS). Figure 5 shows a phase plot of the effective pressure vs. the mass density for the highest-resolution case. One can see that the average of the effective pressure for a given mass density closely follows the isothermal relation Pρ. Although the mean turbulent pressure 2/3ρK is small compared to the thermal pressure for the resolution Δ = L/256, the intermittency of turbulence can locally produce an effective pressure that exceeds the thermal pressure by one order of magnitude. Consequently, the turbulent pressure can become important for compressible turbulence, particularly if there are other sources than the turbulent cascade. As an example, turbulent feedback in galaxy simulations is discussed in Section 6.2. In addition to the turbulent pressure, the non-diagonal turbulent stresses \(\tau _{ij}^{\ast}\) act on the resolved flow. In the case of the eddy-viscosity closure, \(\tau _{ij}^{\ast}\) causes a diffusion effect on top of the numerical diffusion, which occurs regardless of the compressibility of the flow.Footnote 16 For strongly diffusive numerical schemes, this effect is marginal. For high-resolution schemes, on the other hand, the explicit turbulent stresses in LES can become significant. Moreover, the non-linear term in Eq. (54) modifies the diffusion-like tensor in the case of supersonic turbulence. Of course, adding the turbulent stresses in LES does not merely degrade a high-resolution scheme to a more diffusive scheme because the diffusion is linked to the non-linear turbulent interactions across the grid scale. For non-turbulent flows, the turbulent stresses should vanish if the SGS model is consistent.

Figure 3:
figure 3

Temporal evolution of the spatially averaged SGS turbulence energy (left) and the dissipation rate (right) for forced supersonic turbulence. The grid scale Δ decreases from L/32 (light colour) to L/256 (full colour). Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

Figure 4:
figure 4

Time-averaged mean values of the SGS turbulence energy in LES with different resolutions (dots) and power-law fits (dashed lines) for solenoidal and compressive forcing. Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

Figure 5:
figure 5

Phase diagram of the effective pressure defined by Eq. (69) vs. the mass density, showing a two-dimensional probability density function. The spacing of the contour lines is logarithmic. The variables are normalized by their mean values, P0 and ρ0. The black dashed line corresponds to the isothermal relation Pρ, and the red circles indicate the mean effective pressure for given mass density. Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

3.4 Two equation models and gravity

In the framework of the Reynolds-averaged Navier-Stokes equations (RANS), which are equivalent to the filtered Navier-Stokes equation in the limit of a filter scale comparable the integral scale of the flow, the K-ε turbulence model can be used to calculate both the turbulence energy K and the dissipation rate ε. Two inhomogeneous PDEs of the advection-diffusion type determine K and ε (Pope, 2000). In contrast to the simple dimensional closure (58) with a single coefficient of order unity, the diffusion and source terms in the equation for ε come with several additional closure coefficients. This type of model is commonly used for industrial and environmental flows.

A two equation model of similar structure is proposed for buoyancy-driven flows in Dimonte and Tipton (2006). The model predicts the energy K and characteristic size L of the dominant eddies produced by Rayleigh-Taylor and Richtmeyer-Meshkov instabilities. The evolution of these variables is given by the following two equations (in notation adapted to this review):

$${\partial \over {\partial t}}\rho L + \nabla \cdot (\rho uL) = \rho \sqrt {2K} + {C_L}\rho Ld + \nabla \cdot (\rho {k_L}\nabla L),$$
(70)
$${\partial \over {\partial t}}\rho L + \nabla \cdot (\rho uK) = \Gamma + \tau _{ij}^{\ast}{S_{ij}} - {2 \over 3}\rho Kd - {C_\epsilon }{{{K^{3/2}}} \over L} + \nabla \cdot (\rho {k_K}\nabla K).$$
(71)

The production rate due to the Rayleigh-Taylor instability is basically given by Γ ∝ ρ √2 K geff. For each grid cell, the buoyant acceleration geff is assumed to depend on the density contrast across the cell faces, the length scale L, and the components of gravity along the three coordinate axes. Apart from differences in the determination of geff, the above expression for Γ is the same as in Eq. (128) for the enhanced production of turbulence at the interface between low-density and high-density material, which is applied to the propagation of turbulent flame fronts in thermonuclear supernovae (see Section 6.1).

The K-L model outlined above was adopted as an SGS model for simulations of turbulence driven by active galactic nuclei (AGNs) in galaxy clusters (Scannapieco and Brüggen, 2008). In this case, turbulence is thought to be stirred by hot bubbles rising due to their buoyancy in the intracluster medium (ICM). These bubbles originate from the AGNs in the cluster. For this reason, production through the turbulent cascade is set to zero, i.e., \(\tau _{ij}^{\ast} = 0\).Footnote 17 In contrast to LES based on the consistent decomposition derived in Section 2.1, it follows from the very concept of the K-L model that K cannot be interpreted as the kinetic energy associated with velocity fluctuations below the grid scale. Since K is the kinetic energy associated with the dominant eddies driven by the RT instability on a length scale L, where L is a dynamical variable that can become large compared to the grid scale, K generally encompasses some fraction of the numerically resolved turbulence energy on top of the SGS turbulence energy (if L falls below the grid resolution, on the other hand, K will represent only a fraction of the SGS turbulence energy). Consequently, the K-L model works as a hybrid model in these simulations, which stands in between RANS and LES. As a result, turbulent fluctuations in the RT-unstable bubbles are smeared out over length scales smaller than L and only the coherent structures on larger scales are captured in Scannapieco and Brüggen (2008). LES, on the other hand, always resolve the small-scale structure of turbulent flows down to the grid scale. LES in the sense discussed here, however, would break down if there is no physical turbulence on length scales of the order of the grid scale and below (i.e., if the scale of viscous dissipation is within the range of numerically resolved scales).

For the general case of self-gravitating turbulent gas, no satisfactory closure for Γ has been found yet. A conceptual difficulty is that the acceleration caused by gravity is genuinely anisotropic, while SGS models such as the turbulence energy model are based on local isotropy. The usual solution to this problem is to resolve the flow down to length scales that are not strongly affected by gravity. In AMR simulations, this is achieved by imposing a Truelove-like resolution criterion such that a sufficiently large ratio between the local Jeans length and the grid scale is maintained (Truelove et al., 1997; Federrath et al., 2011). Since the density of gravitationally unstable gas would increase indefinitely, excess mass is usually dumped into sink particles at the highest refinement level (Krumholz et al., 2004; Federrath et al., 2010a; Wang et al., 2010). Thereby, collapsing gas is decoupled from the numerically computed gas dynamics. In a certain sense, a sink particle is nothing but an SGS model for a self-gravitating overdense cloud that collapses down to scales below the minimal grid scale. Despite being a crude model, sink particles are a reasonable approximation to collapsed clouds because they mainly interact through accretion (i.e., mass accumulation) with the numerically resolved gas dynamics. A more complex situation is encountered if the objects represented by the sink particles produce feedback onto the gas. An example is stellar feedback in galaxy simulations, which can be treated with the approach discussed in Section 6.2.

4 Determination of Closure Coefficients

One of the basic assumptions of the Kolmogorov theory is that turbulence is statistically self-similar in the inertial subrange (see, for example, Frisch, 1995). With regard to subgrid-scale closures, the self-similarity of turbulence implies that dimensionless coefficients such as C ν in Eq. (47) should be independent of the chosen filter scale. This is not only a necessary condition for the feasability of LES, but it also allows for the calibration of closure coefficients by explicitly filtering turbulence data. Since closures do not exactly match SGS terms, an improved approximation can be achieved by so-called dynamic procedures, which estimate coefficients from properties of the numerically resolved flow under the assumption of local self-similarity.

4.1 Hierarchical filtering

As a formal framework, let us consider an infinite series of isotropic and time-independent filter operators 〈 〉 n . Each filter is defined by a kernel G n (r) with filter length Δ n (see Section 2). We shall assume that Δ0L, where L is the integral length scale of the flow, and

$$\forall n \in {\mathbb{N}_0}:{G_n}(x) = {\lambda ^3}{G_{n - 1}}(\lambda x),\quad {\rm{where}}\;\lambda {\rm{ > 1,}}$$
(72)

i.e., 〈 〉 n for n = 0, 1, 2, … is a self-similar hierarchy of filters. Typical examples are the box filter defined by Eq. (5) or the Gaussian filter, which has the kernel (Sagaut, 2006; Pope, 2000)

$${G_n}(r) = {\left({{6 \over {\pi \Delta _n^2}}} \right)^{3/2}}\exp \left({{{6{r^2}} \over {\Delta _n^2}}} \right)$$
(73)

for isotropic filter lengths Δ n = Δ0/λn. Since G n (r) → δ(r) in the limit n → ∞, 〈 〉 is the identity operator. Because filtering in physical space corresponds to a multiplication with the transfer function of the filter in Fourier space, it follows that

$${\langle {\langle q\rangle _m}\rangle _n} \simeq {\langle q\rangle _n}\quad {\rm{if}}\;{\Delta _n} \gg {\Delta _m}.$$
(74)

For Gaussian filters, the validity of this approximation becomes immediately clear by calculating the product of the transfer functions:

$${\hat G_m}(k){\hat G_n}(k) = \exp \left[ { - {{{k^2}(\Delta _m^2 + \Delta _n^2)} \over {24}}} \right] \simeq \exp \left[ { - {{{k^2}\Delta _n^2} \over {24}}} \right] = {\hat G_n}(k).$$
(75)

We can now apply the scale separation of the Navier-Stokes equations introduced in Section 2 at different levels of the filter hierarchy. In particular, the filtered density field at the n-th level is 〈ρ n , and the Favre-filtered velocity is given by

$${\tilde u^{[n]}} = {{{{\langle \rho u\rangle }_n}} \over {{{\langle \rho \rangle }_n}}}.$$
(76)

By filtering twice at levels m and n, we obtain

$${\tilde u^{[m][n]}}{\langle {\langle \rho \rangle _m}\rangle _n} = {\langle {\langle \rho \rangle _m}{u^{[m]}}\rangle _n} = {\langle {\langle \rho \rangle _m}\rangle _n}.$$
(77)

If the n-th level is much coarser than the m-th level, the asymptotic relation (74) for Δ n ≫ Δ m implies

$${\tilde u^{[m][n]}}{\langle {\langle \rho \rangle _m}\rangle _n} = {\langle {\langle \rho u\rangle _m}\rangle _n} \simeq {\langle \rho u\rangle _n} = {\tilde u^{[n]}}{\langle \rho \rangle _n}.$$
(78)

The turbulence stress tensor on the length scale Δ n of the n-th filter is defined by

$${\tau ^{[n]}} = - {\langle \rho u \otimes u\rangle _n} + {\langle \rho \rangle _n}{\tilde u^{[n]}}{\tilde u^{[n]}}.$$
(79)

The stress tensors for two filter levels m and n, where Δ m < Δ n , are related by the Germano identity (see Section 3.3.3 in Sagaut, 2006 and Germano, 1992; Schmidt, 2004):

$${\tau ^{[m][n]}} = {\langle {\tau ^{[m]}}\rangle _n} + {\tau ^{[m,n]}}.$$
(80)

The stress tensor associated with the double-filtered variables is defined by

$$\begin{array}{*{20}c} {\tau ^{[m][n]} = - \left\langle {\left\langle {\rho u \otimes u} \right\rangle _m } \right\rangle _n + \left\langle {\left\langle \rho \right\rangle _m } \right\rangle _n \tilde u^{[m][n]} \tilde u^{[m][n]} = } \\ { - \left\langle {\left\langle {\rho u \otimes u} \right\rangle _m } \right\rangle _n + \frac{{\left\langle {\left\langle {\rho u} \right\rangle _m } \right\rangle _n \otimes \left\langle {\left\langle {\rho u} \right\rangle _m } \right\rangle _n }} {{\left\langle {\left\langle \rho \right\rangle _m } \right\rangle _n }}} \\ \end{array} $$
(81)

and

$${\tau ^{[m,n]}} = - {\langle {\langle \rho \rangle _m}{\tilde u^{[m]}} \otimes {\tilde u^{[m]}}\rangle _n} + {{{{\langle {{\langle \rho \rangle }_m}{{\tilde u}^{[m]}}\rangle }_n} \otimes {{\langle {{\langle \rho \rangle }_m}{{\tilde u}^{[m]}}\rangle }_n}} \over {{{\langle {{\langle \rho \rangle }_m}\rangle }_n}}}$$
(82)

is the Leonard stress tensor, which is associated with velocity fluctuations in the intermediate range of length scales Δ m ≤ ℓ ≤ Δ n . The Germano identity also holds for two arbitrary filters in the hierarchy. In the limit Δ n ≫ Δ m , the contribution from 〈τ[m] n becomes negligible and

$${\tau ^{[m][n]}} \simeq {\tau ^{[m,n]}} \simeq {\tau ^{[n]}},$$
(83)

where the second relation follows from Eq. (74). As a consequence, the turbulent stresses associated with the scale Δ n are not sensitive to the flow structure on much smaller scales. In particular, it follows that K[n]K[m, n] if Δ m ≪ Δ n .

In Schmidt (2004), Schmidt et al. (2006b), and Schmidt and Federrath (2011), Gaussian filters (see Eq. 73) are applied to data from ILES of forced compressible turbulence for the verification of closures. The following line of reasoning is of central importance for estimating closure coefficients from finite-resolution data. To begin with, let us assume that ρ and u are the physical density and velocity fields. Let us further assume that the implicit filter of the ILES correspond to the filter level m = I, i.e., 〈ρI and ũ[I] represent the numerically computed density and velocity fields. Now, if the numerical data are coarse-grained by an explicit filter 〈 〉 n in the inertial subrange, the turbulence stress tensor τ[I, n] defined by Eq. (82) can be calculated. But a closure for the turbulent stresses on the the length scale Δ n applies to τ[n], which is defined by Eq. (79). If Δ n is reasonably large compared to Δ I , however, one can make use of the approximation τ[n]τ[I,n] (see relation 83 for m = I). Owing to Eq. (78), the distinction between, on the one hand, the physical densities and velocities (or DNS data) and, on the other hand, the ILES data becomes immaterial. It is thus possible to calculate coarse-grained eddy-viscosity coefficients as:

$${C_\nu } \simeq {{\tau _{ij}^{[n]{\ast}}S_{ij}^{[n]}} \over {{{\langle \rho \rangle }_n}{\Delta _n}{{({K^{[n]}})}^{1/2}}\vert {S^{[n]{\ast}}}{\vert ^2}}}.$$

Of course, since the eddy-viscosity closure is not exact, the value of C ν varies. But the mean value turns out to be roughly 0.05 for different filter lengths and simulation parameters (Schmidt, 2004), in good agreement with other estimates in the literature. The same method was used to determine the coefficient C κ for the gradient-diffusion closure (65). The result C ν ≈ 0.4 implies a Prandtl number C κ /C ν around 10, contrary to the common assumption that the kinetic Prandtl number is of the order unity (Sagaut, 2006).

4.2 Dynamic procedures

Subgrid-scale models in their standard form apply to statistically stationary and isotropic turbulence. But turbulent flows in nature often deviate from this idealization: In terrestrial applications, flow inhomogeneities are inevitably caused by boundary conditions (“walls”). In astrophysics, one of the major energy sources is gravity. It causes matter to clump (galaxies and clusters) or to move under the action of central gravitational fields (stars), which produces inherently inhomogeneous and anisotropic flows. For example, turbulent convection in stars introduces a vertical anisotropy of the flow. Turbulence driven by violent energy release (supernovae) can also be highly inhomogeneous.

One of the solutions to this problem is to localize closures, i.e., to calculate local closure co-efficients. This requires local estimators that take properties of the flow in some small region as input. Obviously, this works only if the size of this region is not significantly affected by the flow inhomogeneity on larger scales. In other words, the flow must be asymptotically homogeneous and isotropic at least on length scales of the order of the grid scale. In this case, a so-called test filter 〈 〉T can be applied in LES, with a filter length ΔT that is a small multiple of the grid scale Δ. Test filters are usually implemented as discrete filters over several grid cells (see Section 2.3.2 in Garnier et al., 2009). A multi-dimensional test filter can be composed as a succession of one-dimensional filters.Footnote 18 The test filter length ΔT can be adjusted by varying the weights of the cells. An optimal ratio γT = ΔT/Δ is given by the closest match between the filter transfer functions of the discrete and analytical box filters with filter length ΔT (Vasilyev et al., 1998). For instance, a test filter with γT = 2.771 is optimal if a five-point stencil is used in each spatial dimension (Schmidt, 2004).

By identifying Δ m with Δ and Δ n with ΔT, the Germano identity (80) allows us to express the turbulence stress tensor on the length scale of the test filter as the sum of the test-filtered SGS turbulence stress tensor and the Leonard tensor for the intermediate velocity fluctuations (see also Section 4.3 in Sagaut (2006)):

$${\bf{T}} = {\langle \tau \rangle _{\rm{T}}} + {\bf{L}}.$$
(84)

Here, the Leonard tensor L associated with the test filter is defined by

$${\bf{L}} = - {\langle \rho u \otimes u\rangle _{\rm{T}}} + {{{{\langle \rho u\rangle }_{\rm{T}}} \otimes {{\langle \rho u\rangle }_{\rm{T}}}} \over {{{\langle \rho \rangle }_{\rm{T}}}}},$$
(85)

where we use the simplified notation ρ and u for the density and velocity on the grid scale, as in Section 3.1. Because of the scale-invariance of turbulence, Germano et al. (1991) proposed that the eddy-viscosity closure (49) holds for both T and τ. In the case of the Smagorinsky model (see Eq. 62 for νsgs), the corresponding tensors are:

$$\tau _{ij}^{\ast} = 2\rho {({C_{\rm{S}}}\Delta)^2}\vert S\vert {S_{ij}} =:C_{\rm{S}}^2{\beta _{ij}},$$
(86)
$$T_{ij}^{\ast} \simeq 2{\rho _{\rm{T}}}{({C_{\rm{S}}}{\Delta _{\rm{T}}})^2}\vert {S_{\rm{T}}}\vert {({S_{\rm{T}}})_{ij}} =:C_{\rm{S}}^2{\alpha _{ij}}.$$
(87)

The rate-of-strain tensor (ST) ij at the test filter level is given by the symmetrized derivative of the test-filtered numerically resolved velocity field, ∂ i u j T, analogous to Eq. (30). The variable CS(x, t) needs to be determined. This can be achieved by substituting the above expressions for \(\tau _{ij}^{\ast}\) and \(T_{ij}^{\ast}\) into the trace-free part of the Germano identity (84), which implies

$$L_{ij}^{\ast} \simeq C_{\rm{S}}^2{\alpha _{ij}} - {\langle C_{\rm{S}}^2\beta \rangle _{\rm{T}}}.$$
(88)

Under the assumption that CS varies only little over the smoothing length of the test filter, one can set \({\langle C_{\rm{s}}^2{\beta _{ij}}\rangle _{\rm{T}}} \simeq C_{\rm{s}}^2{\langle {\beta _{ij}}\rangle _{\rm{T}}}\). Since L* ij can be evaluated from Eq. (85), minimalization of the residual error between \(L_{ij}^{\ast}\) and the expression on the right-hand side of Eq. (88) yields

$$C_{\rm{S}}^2 = {{{m_{ij}}L_{ij}^{\ast}} \over {{m_{ij}}{m_{ij}}}},$$
(89)

where m ij = αij − 〈β ij T. This is the Germano-Lilly dynamic procedure, which was applied, for example, in LES of turbulent channel flows (Piomelli, 1993). In principle, this procedure could also be applied to the non-equilibrium model with the turbulent viscosity defined by Eq. (47). In this case, the turbulence energy associated with T is given by the contracted Germano identity,

$$ - {1 \over 2}{T_{ii}} = {\langle \rho K\rangle _{\rm{T}}} + {\rho _{\rm{T}}}{K_{\rm{T}}},$$

where ρK = − τ ii /2 and ρTKT = − L ii /2.

However, the dynamic procedure as outlined above has several caveats. In particular, the assumption of negligible variation of CS over the the test filter length is found to be violated significantly. Moreover, CS diverges if m ij vanishes. Consequently, several attempts were made to improve the dynamic procedure (Liu et al., 1994; Piomelli and Liu, 1995; Ghosal et al., 1995). A particularly simple modification was found by analyzing experimental measurements of turbulent velocity fluctuations in consecutive wave number bands [k n −1, k n ], corresponding to a hierarchy of filters. By explicitly evaluating the turbulent stresses τ[n] associated with the wave numbers k n = π n , the correlations with localized closures were verified. Although some correlation between the turbulent stresses at different filter levels was found, the correlation of τ[n][n−1] with the Leonard stresses L[n, n−1] turned out to be significantly better. This observation can be understood as a consequence of the locality of the energy transfer (Kraichnan, 1976; Sagaut, 2006), i.e., the energy transfer across a certain wave number k is mainly caused by interactions in the narrow spectral band [½k, 2k]. With regard to test filtering in LES, this implies that the eddy-viscosity closure should be applied to L in place of T. The localized coefficient C ν (x, t) of the turbulent viscosity is then given by (Kim and Menon, 1999; Schmidt, 2004; Schmidt et al., 2006b; Röpke and Schmidt, 2009)

$${C_\nu } \simeq {{L_{ij}^{\ast}{{({S_{\rm{T}}})}_{ij}}} \over {{\rho _{\rm{T}}}{\Delta _{\rm{T}}}K_{\rm{T}}^{1/2}\vert S_{\rm{T}}^{\ast}{\vert ^2}}},$$
(90)

where KT = −L ii /(2ρT) is the resolved kinetic energy on length scales Δ ≤ ℓ ≤ Δ T . Substitution of the above expression for C ν into Eq. (50) for the localized rate of production yields

$$\Sigma = {\tau _{ij}}{S_{ij}} \simeq {{\rho \Delta } \over {{\rho _{\rm{T}}}{\Delta _{\rm{T}}}}}{\left({{K \over {{K_{\rm{T}}}}}} \right)^{1/2}}{\left({{{\vert {S^{\ast}}\vert } \over {\vert S_{\rm{T}}^{\ast}\vert }}} \right)^2}L_{ij}^{\ast}{({S_{\rm{T}}})_{ij}} - {2 \over 3}\rho Kd.$$
(91)

The above formula was used for simulations of turbulent thermonuclear combustion in white dwarfs (see Section 6.1). A generalization of the dynamic procedure to localize both 1 and 2 in the closure (54) would be straightforward, but has not been applied so far. In this case, a linear system in the coefficients C1 and C2 has to be solved to minimalize the residual.

For ILES of subsonic and transonic turbulence produced by stochastic forcing in periodic boxes (Schmidt et al., 2006a), the enhanced fidelity of the localized closure (91) can be verified by coarse-graining the data with hierarchical Gaussian filters 〈 〉 n as explained in Section (4.1) (Schmidt et al., 2006b). Let us first consider the turbulence energy flux produced by anisotropic shear on the length scale Δ n , assuming a turbulent viscosity with a constant coefficient 〈C ν 〉, which is obtained by averaging Eq. (90) over the domain of the flow. If the closure with this coefficient were exact, we would have

$${\Sigma ^{[n]({\rm{eddy)}}}} + {2 \over 3}{\langle \rho \rangle _n}{K^{[n]}}{d^{[n]}} = \langle {C_\nu }\rangle {\langle \rho \rangle _n}{\Delta _n}\sqrt {{K^{[n]}}} \vert {S^{[n]\ast}}{\vert ^2}.$$
(92)

Here, the divergence term is added to Σ[n](eddy) to express the flux associated with the trace-free rate-of-strain tensor S[n]*. However, coarse-grained data show that Eq. (92) is not very well-satisfied. The probability density functions of the expression on the right-hand side and the

Figure 6:
figure 6

Left: comparison of the probability density functions of the coarse-grained turbulence energy flux due to anisotropic shear with different closures (left). Right plot: probability density functions of the localized closure coefficient C ν obtained by test-filtering different coarse-grained ILES. The ratio of the test filter length to the coarse-graining length is γT. The static closure refers to the case with a constant coefficient. The parameters of the random forcing are the characteristic Mach number V/c0 and the weight ζ of the Helmholtz decomposition into solenoidal and compressive modes (Schmidt, 2004; Schmidt et al., 2006a, b).

explicitly calculated energy flux on the left-hand side are compared in the left plot in Figure 6. One can see that the latter is negative in about 20% of the domain (purple line), corresponding to backscattering from smaller to larger scales. This is excluded by the eddy-viscosity closure with a fixed coefficient (light blue line). In this case, negative values of the total energy flux Σ[n] (eddy) are solely due to the divergence term. The bias toward positive fluxes can be avoided by localizing the eddy-viscosity closure (Schmidt et al., 2006b). For a test filter 〈 〉 n −1 with filter length γT = Δ n −1 n > 1, the energy flux is given by the following analogue of Eq. (91):

$${\Sigma ^{[n]({\rm{locl)}}}} + {2 \over 3}{\langle \rho \rangle _n}{K^{[n]}}{d^{[n]}} = {{{{\langle \rho \rangle }_n}{\Delta _n}} \over {{{\langle \rho \rangle }_{n - 1}}{\Delta ^{[n - 1]}}}}{\left({{{{K^{[n]}}} \over {{K^{[n - 1]}}}}} \right)^{1/2}}{\left({{{\vert {S^{[n]\ast}}\vert } \over {\vert {S^{[n - 1]\ast}}\vert }}} \right)^2}L_{ij}^{[n - 1]\ast}S_{ij}^{[n - 1]}.$$
(93)

The probability density functions that are plotted for different test filtering ratios γT in Figure 6 (left plot) indicate a substantially improved match between the localized closure and the explicitly calculated energy flux. Indeed, distributions of the localized closure coefficients show that C ν has a negative branch (see right plot in Figure 6). This result suggests an improvement due to the dynamic procedure even in the case of homogeneous turbulence. On the other hand, the mean values of C ν appear to be fairly robust for different forcing parameters.

4.3 Global least squares method

Closures can also be tested by analyzing correlations. This allows for the calibration of the closure coefficients by least squares minimalization of the integrated residual (Schmidt and Federrath, 2011). For example, let us consider a generic closure with a single coefficient C1 at the n-th filter level:

$${C_1}{f^{[n]({\rm{cls}})}} = {\Sigma ^{[n]({\rm{cls}})}} + {2 \over 3}{\langle \rho \rangle _n}{K^{[n]}}{d^{[n]}},$$
(94)

The global residual of the explicitly computed turbulence energy flux \({\sum ^{[n]}} = \tau _{ij}^{[n]}S_{ij}^{[n]}\), where \(\tau _{ij}^{[n]}\) is defined by Eq. (79), can be quantified by the squared error integrated over the whole domain \(\mathcal{V}\) of the turbulent flow:

$${\rm{er}}{{\rm{r}}^2}({C_1}) = \int_\mathcal{V} {{{\left\vert {\Sigma [n] + {2 \over 3}{{\langle \rho \rangle }_n}{K^{[n]}}{d^{[n]}} - {C_1}{f^{[n]({\rm{cls}})}}} \right\vert }^2}{{\rm{d}}^3}x.} $$
(95)

The minimum of err2(C1) is obtained by setting the derivative with respect to C1 equal to zero:

$${C_1} = {{\int_\mathcal{V} {{f^{[n]({\rm{cls}})}}[{\Sigma ^{[n]({\rm{cls}})}} + {\textstyle{2 \over 3}}{{\langle \rho \rangle }_n}{K^{[n]}}{d^{[n]}}]} \;{{\rm{d}}^3}x} \over {\int_\mathcal{V} {\vert {f^{[n]({\rm{cls}})}}{\vert ^2}{{\rm{d}}^3}x} }}.$$
(96)

In contrast to the dynamic procedure, the resulting closure coefficients are constants.

The method of least squares described above is applied in Schmidt and Federrath (2011) to various ILES of supersonic isothermal turbulence produced by stochastic forcing (Schmidt et al., 2009; Federrath et al., 2010b). To coarse-grain the data, a Gaussian filter with a smoothing length Δ4 = L/16 = 32Δ I is used, where Δ I = Δ9 ≡ Δ is the grid resolution of the ILEs.Footnote 19 Since the filter length is large compared to the grid resolution in this case, it is advantageous to apply the filter operation in Fourier space. For the eddy-viscosity closure, we have

$${f^{[n]({\rm{cls}})}} = {\Delta _n}{\langle \rho \rangle _n}\sqrt {2{K^{[n]}}} \vert {S^{[n]\ast}}{\vert ^2}.$$

The closure coefficients following from Eq. (96) are, for instance, C1 ≈ 0.102 for the 10243 ILES with purely solenoidal (divergence-free) forcing, and C1 ≈ 0.092 in the case of compressive (rotation-free) forcing. The corresponding value of C ν = √2C1 is about 0.14. When comparing this value to the results in Schmidt (2004); Schmidt et al. (2006b) (see Section 4.1), one has to bear in mind not only that a different method is applied to determine the coefficients, but also that the turbulence properties differ substantially.

The correlation diagram for Σ[4](cls), with the least-squares coefficient C1, versus Σ[4] in the case of solenoidal forcing is shown in Figure 7 (left plot). The overall correlation is actually quite good. A quantitative measure is the correlation coefficient

$$\begin{array}{*{20}c} {corr\left[ {\Sigma ^{[n]} ,\Sigma ^{[n](cls)} } \right] = } \\ {\frac{{\int_\nu {\left( {\Sigma ^{[n]} - \left\langle {\Sigma ^{[n]} } \right\rangle } \right)\left( {\Sigma ^{[n](cls)} - \left\langle {\Sigma ^{[n][cls]} } \right\rangle } \right)} d^3 x}} {{std\left( {\Sigma ^{[n]} } \right)std\left( {\Sigma ^{[n](cls)} } \right)}},} \\ \end{array} $$
(97)

where std denotes the standard deviation and the angle brackets indicate an average over the whole domain. The correlation coefficients of the eddy-viscosity closure are found to be 0.95 and 0.93 for solenoidal and compressive forcing, respectively (Schmidt and Federrath, 2011). However, it becomes apparent that the closure breaks down for negative fluxes. This corresponds to the bias of the probability density function for the static closure with an averaged coefficient in Figure 6 (left plot). Rather than applying the dynamic procedure, it is shown in Schmidt and Federrath (2011) that the determinant closure (52) results in a largely improved approximation of negative turbulence energy flux. This is a consequence of the varying sign of the determinant, det S*, while ΔK1/2S*∣2 is positive. However, the scatter of the determinant closure is high, particularly for large positive flux. This is remedied by the non-linear closure (53) for the turbulence stress tensor, which produces an excellent correlation between Σ[n](cls) and Σ[n], as demonstrated by the right plot in Figure 7. Since the formulae are analogous to the eddy-viscosity closure, we refer to Schmidt and Federrath (2011) for details. The correlation coefficients are 0.991 for both solenoidal and compressive forcing.

Figure 7:
figure 7

Correlations of the coarse-grained turbulence energy flux with the eddy-viscosity (left) and non-linear (right) closures for supersonic isothermal turbulence produced by solenoidal forcing. The applied filter length is 32Δ, where Δ is the grid resolution. The average prediction of the closure for each bin is indicated by the blue dots. Image reproduced with permission from Schmidt and Federrath (2011), copyright by ESO.

However, as explained in Section 3.1, the purely non-linear closure is not suitable for an SGS model. This is why the least squares method was applied to the generalized closure with two coefficient, C1 and C2. For

$${C_1}{f^{[n]({\rm{cls}})}} + {C_2}{g^{([n]{\rm{cls}})}} = {\Sigma ^{[n]({\rm{cls}})}} + {2 \over 3}{\langle \rho \rangle _n}{K^{[n]}}{d^{[n]}},$$
(98)

the closure coefficients are given by the linear system of equations

$$\begin{array}{*{20}c} {\left( {\int_\mathcal{V} {\left| {f^{[n](cls)} } \right|} ^2 d^3 x} \right)C_1 + \left( {\int_\mathcal{V} {f^{[n](cls)} g^{[n](cls)} } d^3 x} \right)C_2 = } \\ {\int_\mathcal{V} {f^{[n](cls)} } \left( {\Sigma ^{[n](cls)} + \frac{2} {3}\left\langle \rho \right\rangle _n K^{[n]} d^{[n]} } \right)d^3 x,} \\ \end{array} $$
(99)
$$\begin{array}{*{20}c} {\left( {\int_\nu {f^{[n](cls)} g^{[n](cls)} d^3 x} } \right)C_1 + \left( {\int_\nu {\left| {g^{[n](cls)} } \right|^2 d^3 x} } \right)C_2 = } \\ {\int_\nu {g^{[n](cls)} \left( {\Sigma ^{[n](cls)} + \frac{2} {3}\left\langle \rho \right\rangle _n K^{[n]} d^{[n]} } \right)d^3 x,} } \\ \end{array} $$
(100)

where

$${f^{[n]({\rm{cls}})}} = {\Delta _n}{\langle \rho \rangle _n}\sqrt {2{K^{[n]}}} \vert {S^{[n]\ast}}{\vert ^2}\quad {\rm{and}}\quad {g^{[n]({\rm{cls}})}} = - 4{\langle \rho \rangle _n}{K^{[n]}}{{u_{i,k}^{[n]}u_{j,k}^{[n]}S_{ij}^{[n]\ast}} \over {\vert \nabla \otimes {u^{[n]}}{\vert ^2}}}.$$
(101)

The solution is C1 ≈ 0.02 and C2 ≈ 0.7, with a correlation coefficient 0.990 (Schmidt and Federrath, 2011). These coefficients also yield good approximations to the turbulence energy flux for isothermal and adiabatic turbulence simulations at lower Mach numbers (Schmidt et al., 2009; Schmidt and Federrath, 2011). Moreover, the coefficients appear to vary only little with the filter length, at least in the range from Δ3 = 64Δ to Δ5 = 16Δ. Choosing longer or shorter filter lengths is not sensible because of the influence of the forcing (Δn must be small compared to Δ0 = L) and numerical dissipation (Δn ≫ Δ).

5 Adaptive Methods

The most powerful technique for finite-volume codes to resolve localized and anisotropic structures in a flow is adaptive mesh refinement (AMR) (Berger and Oliger, 1984; Berger and Colella, 1989). Even with AMR, however, it is generally not possible to fully resolve turbulence. This entails the problem that the numerically resolved and unresolved turbulence energy fractions vary as regions are refined or de-refined. In the following, it shown how to address this problem in adaptively refined LES. In principle, global energy and momentum conservation can be achieved, while reducing the need for artificial changes in the internal energy, which is the standard method to restore energy consistency between different refinement levels in AMR simulations. Apart from that, localized and anisotropic flow structures pose the problem that SGS models with constant coefficients introduce systematic errors because they are usually calibrated for statistically stationary and isotropic turbulence. Shear-improved SGS models can alleviate this problem by adjusting the non-linear energy transfer across the grid scale to local flow conditions. For example, this is possible by applying an adaptive temporal filter, the so-called Kalman filter.

5.1 Energy and momentum conservation in AMR simulations

In AMR simulations, data have to be transferred between different refinement levels by conservative interpolation or averaging. For example, if a region is refined, data from coarser grids are interpolated to finer grids. The same operation is used for filling ghost cells at the boundaries between a finer and a coarser level, which is required to compute fluxes through the faces of adjacent finer and coarser cells. Moreover, block-structured AMR codes usually average down the data from the highest-level grid to the all coarser levels. The mass density, momentum, and energy variables at two levels, say, l and l + 1, are in the simplest case related by

$${\rho _{{\rm{crs}}}}: = \overline \rho = {1 \over N}\sum\limits_n {{\rho _n},}$$
(102)
$${(\rho U)_{{\rm{crs}}}}: = \overline {\rho U} = {1 \over N}\sum\limits_n {{{(\rho U)}_n}},$$
(103)
$${E_{{\rm{crs}}}}: = \overline {\rho E} = {1 \over N}\sum\limits_n {{{(\rho E)}_n} = {1 \over N}\sum\limits_n {\left[ {{{(\rho e)}_n} + {1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}}} \right]},} $$
(104)

where N grid cells at level l + 1 are summed up to a single value in a coarse cell at level l. Obviously, these relations guarantee mass, momentum, and energy conservation. For more sophisticated interpolation schemes, the fine-grid values have different weights w n in the above sums. Without loss of generality, we assume w n = 1 in the following. Now, if we work out the internal energy at the coarser level, we obtain

$${e_{{\rm{crs}}}} = {E_{{\rm{crs}}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}} = \overline {\rho e} + \Delta (\overline {\rho K}),$$

where the energy difference

$$\Delta (\overline {\rho K}): = {1 \over N}\sum\limits_n {{1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}}}$$
(105)

is generally non-zero because of Eq. (103). This implies that

$${e_{{\rm{crs}}}} \neq \overline {\rho e} = {1 \over N}\sum\limits_n {{{(\rho e)}_n}.} $$

To put it another way, the kinetic energy differences between refinement levels have to be compensated by numerical cooling or heating in order to maintain energy conservation.

This can be alleviated by incrementing the SGS turbulence energy by \(\Delta (\overline {\rho K})\) when the cutoff scale is shifted from Δ l +1 to Δ l = rΔ l +1, where r > 1 is the refinement ratio:

$${(\rho K)_{{\rm{crs}}}} = \overline {\rho K} + \nabla (\overline {\rho K}),\quad {\rm{where}}\quad \overline {\rho K} = {1 \over N}\sum\limits_n {{{(\rho K)}_n}.}$$
(106)

We then have conservation of the total (resolved plus unresolved) kinetic energies,

$${1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}} + {(\rho K)_{{\rm{crs}}}} = {1 \over N}\sum\limits_n {\left[ {{1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}} + {{(\rho K)}_n}} \right]},$$
(107)

and \({e_{{\rm{crs}}}} = \overline {\rho e} \).

As demonstrated in Schmidt et al. (2014), Eq. (107) works very well for fully developed homogeneous turbulence, but leads to erroneous projections of the SGS turbulence energy to coarser grids if the flow structure is strongly inhomogeneous. A particular problem is posed by non-turbulent bulk flows such as gas accretion into gravitational wells. In this case, applying the increment \(\Delta (\overline {\rho K})\) defined by Eq. (105) can greatly overestimate the energy difference associated with turbulent velocity fluctuations. A tentative solution to this problem follows along similar lines as the re-scaling procedure used in Maier et al. (2009). By extrapolating the turbulent energy on the length scale Δ l +1 via a power law to Δ l , we have

$${(\rho K)_{{\rm{crs}}}} = \overline {\rho K} {\left({{{{\Delta _l}} \over {{\Delta _{l + 1}}}}} \right)^{2\eta }} = \overline {\rho K} \;{r^{2\eta }}.$$
(108)

Substitution into Eq. (106) yields

$$\Delta (\overline {\rho K}) = ({r^{2\eta }} - 1)\overline {\rho K},$$

where η = 1/3 in the case of Kolmogorov scaling. A shortcoming of this estimate is, of course, that the turbulent velocity fluctuations follow power-law scaling only in a statistical sense. To avoid an overshoot of \(\Delta (\overline {\rho K})\) if turbulence dominates the energy difference between levels, it is necessary to set

$$\Delta (\overline {\rho K}) = \min \left[ {({r^{2\eta }} - 1)\overline {\rho K},{1 \over N}\sum\limits_n {{1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}}} } \right].$$
(109)

Now, energy conservation can only be satisfied if \(\Delta (\overline {\rho K})\) is complemented by a correction of the internal energy,

$$\Delta (\overline {\rho e}) = {1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}} - \Delta (\overline {\rho K}).$$
(110)

As a result, we have

$${(\rho K)_{{\rm{crs}}}} = \overline {\rho K} + \Delta (\overline {\rho K}),$$
(111)
$${(\rho e)_{{\rm{crs}}}} = \overline {\rho e} + \Delta (\overline {\rho e}),$$
(112)
$${(\rho E)_{{\rm{crs}}}} = \overline {\rho E} - \Delta (\overline {\rho K}),$$
(113)

and the total energy is conserved because

$${(\rho E)_{{\rm{crs}}}} + {(\rho K)_{{\rm{crs}}}} = \overline {\rho E} + \overline {\rho K}.$$

For

$$({r^{2\eta }} - 1)\overline {\rho K} \ll {1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}}$$

most of the resolved kinetic energy difference is compensated by internal energy, similar to the standard energy correction employed in AMR simulations without SGS model. Fully turbulent energy compensation, on the other hand, follows as a limiting case if \(\Delta (\overline {\rho K})\) is given by Eq. (105),

For grid refinement, the data from the parent grid are interpolated to fine-grid values \(\rho _n^{\ast},\;{(\rho {U_n})^{\ast}}\), (ρU n )*, etc. such that

$$\begin{array}{*{20}c} {\frac{1} {N}\sum\limits_n {\rho _n^* = \rho _{crs} ,} } & {\frac{1} {N}\sum\limits_n {(\rho U)_n^* = (\rho U)crs,} } & {} \\ {\frac{1} {N}\sum\limits_n {(\rho E)_n^* = E_{crs} ,} } & {\frac{1} {N}\sum\limits_n {(\rho e)_n^* = e_{crs} ,} } & {\frac{1} {N}\sum\limits_n {(\rho K)_n^* = (\rho K)_{crs} .} } \\ \end{array} $$

One can set \({\rho _n} = \rho _n^{\ast},\;{(\rho U)_n} = (\rho U)_n^{\ast}\), and \({(\rho E)_n} = (\rho E)_n^{\ast}\), but then the kinetic energy difference resulting from the conservative interpolation of the momenta has to be compensated. By defining the energy corrections

$$\Delta (\overline {\rho K}) = (1 - {r^{ - 2\eta }}){(\rho K)_{{\rm{crs}}}},$$
(114)
$$\Delta (\overline {\rho e}) = {1 \over N}\sum\limits_n {{1 \over 2}{{{{[(\rho U)_n^{\ast}]}^2}} \over {\rho _n^{\ast}}} - {1 \over 2}} {{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}} - \Delta (\overline {\rho K}),$$
(115)

the interpolated values of the SGS turbulence and internal energies can be adjusted as follows:

$${(\rho K)_n} = (\rho K)_n^{\ast} - {{(\rho K)_n^{\ast}} \over {{{(\rho K)}_{{\rm{crs}}}}}}\Delta (\overline {\rho K}) = {r^{ - 2\eta }}(\rho K)_n^{\ast},$$
(116)
$${(\rho e)_n} = (\rho e)_n^{\ast} - {{(\rho e)_n^{\ast}} \over {{{(\rho e)}_{{\rm{crs}}}}}}\Delta (\overline {\rho e}),$$
(117)

The conservation of total energy follows from

$$\begin{array}{*{20}c} {\frac{1} {N}\sum\limits_n {\left[ {(\rho e)_n + \frac{1} {2}\frac{{(\rho U)_n^2 }} {{\rho _n }} + (\rho K)_n } \right] = } \frac{1} {N}\sum\limits_n {\left[ {(\rho e)_n^* + \frac{1} {2}\frac{{\left[ {(\rho U)_n^* } \right]^2 }} {{\rho _n^* }} + (\rho K)_n^* } \right] - \Delta \left( {\overline {\rho e} } \right) - \Delta \left( {\overline {\rho K} } \right)} } \\ { = e_{crs} + \frac{1} {2}\frac{{(\rho U)_{crs}^2 }} {{\rho _{crs} }} + (\rho K)_{crs} ,} \\ \end{array} $$

where the last equality follows by substituting Eqs. (114) and (115) for \(\Delta (\overline {\rho e})\) and \(\Delta (\overline {\rho K})\), respectively. Equation (116) is formally the same as the correction rule for the SGS turbulence energy in Maier et al. (2009). In contrast to the method outlined above, however, momentum conservation is not fulfilled by Maier et al. (2009) because the full kinetic energy difference between the finer and coarser levels is compensated by internal energy and then the velocities are re-scaled to compensate the power-law correction of the SGS turbulence energy (see also Section 6.3).

5.2 Shear-improved model

For inhomogeneous and non-stationary flows, dynamic procedures can be applied to adjust SGS model coefficients to local flow conditions (see Section 4.2). A completely different idea was put forward by Lévêque to calculate the turbulent stresses in LES of wall-bounded turbulence (Lévêque et al., 2007). Rather than adjusting the eddy-viscosity coefficient C ν , the numerically resolved velocity field is decomposed into a mean flow 〈u〉 and turbulent fluctuations u′. It is then argued that the turbulence energy flux for the Smagorinsky model should linearly depend on the shear associated with the fluctuating component, i.e.,

$$\epsilon \simeq \Sigma = {({C_{\rm{s}}}\Delta \vert S\vert)^2}(\vert S\vert - \vert \langle S\rangle \vert).$$

where ∣〈S〉∣ is the rate of strain of the mean flow. If the flow is laminar, 〈u〉 ≃ u and turbulence production vanishes because ∣S∣ − ∣〈S〉∣ ≃ 0. For developed isotropic turbulence, on the other hand, u′ ≫ u. In this case, the standard Smagorinsky model with Σ ∝ ∣S3 applies. For intermediate cases, the model corrects the energy flux Σ by taking into account interactions of the grid-scale fluctuations with the mean shear. This is why it is called the shear-improved model.

As proposed in Schmidt et al. (2014), this idea can be carried over to the SGS turbulence energy model (see Section 3.3) by defining the shear-improved eddy-viscosity closure as

$${\tau _{ij}} = 2\rho \left({{\nu _{{\rm{sgs}}}}S_{ij}^{\prime \ast} - {1 \over 3}K{\delta _{ij}}} \right),$$
(118)

where \(S_{ij}^{\prime\ast}\) is the trace-free part of the rate-of-strain tensor associated with the fluctuating component of the flow:

$$S_{ij}^{\prime} = {1 \over 2}\left({{{\partial u_i^{\prime}} \over {\partial {x_j}}} + {{\partial u_j^{\prime}} \over {\partial {x_i}}}} \right) = {S_{ij}} - {\langle S\rangle _{ij}}.$$
(119)

The turbulence energy flux is then given by

$$\Sigma = {\tau _{ij}}{S_{ij}} = {C_\nu }\rho \Delta {K^{1/2}}(\vert {S^{\ast}}{\vert ^2} - 2\langle S\rangle _{ij}^{\ast}{S_{ij}}) - {2 \over 3}\rho Kd.$$
(120)

Apart from the divergence term, the expression on the right-hand side is motivated by the generalized Karman-Howarth equation in Lévêque et al. (2007).

A difficulty in implementing the shear-improved model is the computation of the mean flow 〈u〉, which is an ensemble average over the flow velocity u. A practical solution is to find an approximation to the mean flow by smoothing u with a temporal low-pass filter. For statistically stationary turbulence, it is possible to use an exponentially weighted recursive average. In component notation, an estimate for the mean velocity at time t n +1 is calculated as weighted sum of the estimate at the previous time step and the local velocity \(u_i^{(n + 1)}\) for each grid cell (cell indices are omitted):

$${[{u_i}]^{(n + 1)}} = \left({1 - \alpha _i^{(n + 1)}} \right)\;{[{u_i}]^{(n)}} + \alpha _i^{(n + 1)}u_i^{(n + 1)},$$
(121)

The weighing coefficients are defined by

$$\alpha _i^{(n + 1)} = {{2\pi ({t_{n + 1}} - {t_n})} \over {\sqrt 3 {T_c}}}.$$

where Tc is the constant integral time scale of the flow (Cahuzac et al., 2010). Changes occurring on time scales smaller than the smoothing scale Tc are suppressed in [u i ]. Compared to dynamic procedures, this algorithm is very easy to implement and computationally much cheaper. In Lévêque et al. (2007), it is demonstrated that the shear-improved Smagorinsky model with exponential smoothing performs well in LES of plane-channel flow and reproduces data from direct numerical simulations.

However, the simple exponential smoothing algorithm produces an estimate [u i ] that lags behind the ensemble average 〈u i 〉 if the mean flow evolves. To address this problem, the so-called Kalman filtering technique is introduced in Cahuzac et al. (2010). The Kalman filter adapts itself to an unsteady mean flow by dynamically adjusting the weights of the recursive relation (121), depending on the variances of the mean flow evolution, [u i ](n) − [u i ](n−1), and the detected deviation from the mean, \(u_i^{(n)} - {[{u_i}]^{(n)}}\). This is achieved by setting \(\alpha _i^{(n + 1)}\) equal to the so-called Kalman gain \(K_i^{(n + 1)}\), which is defined by the ratio between the error variance of the smoothed component and the total error variance, including the fluctuating component. Since the error variances have to be evaluated at time t(n+1), a predictor-corrector scheme is used:

  1. 1.

    Given the error variance \(P_i^{(n)}\) at time t(n), the prediction for t(n+1) is

    $$P_i^{(n + 1)\ast} = P_i^{(n)} + \sigma _{\delta [{u_i}]}^{2(n)},$$
    (122)

    where

    $$\sigma _{\delta [{u_i}]}^{(n)} = {{2\pi \Delta {t^{(n)}}} \over {\sqrt 3 {T_c}}}{u_c}.$$

    Here it is assumed that the typical correction of the mean flow, δ[u i ](n) = [u i ](n) − [u i ](n−1), is of the order 2πΔt(n)uc/(√3Tc), where Δt(n) = t n t n −1

  2. 2.

    The Kalman gain is then given by

    $$\alpha _i^{(n + 1)} = K_i^{(n + 1)} = {{P_i^{(n + 1)\ast}} \over {P_i^{(n + 1)\ast} + \sigma _{\delta {u_i}}^{2\;(n)}}},$$
    (123)

    where

    $$\sigma _{\delta {u_i}}^{2\;(n)} = \max \left({\left\vert {\delta u_i^{(n)}} \right\vert,\;0.1{u_c}} \right){u_c}$$

    is the contribution of the fluctuating component \(\delta u_i^{(n)} \equiv u_i^{^{\prime}(n)} = u_i^{(n)} - {[{u_i}]^{(n)}}\) to the error variance. The lower bound on \(\sigma _{\delta {u_i}}^{2\;(n)}\) is necessary to obtain non-vanishing fluctuations from an initially smooth flow with [u i ] = u i .

  3. 3.

    The corrected error variance for the next step is given by

    $$P_i^{(n + 1)} = \left({1 - K_i^{(n + 1)}} \right)P_i^{(n + 1)\ast}.$$
    (124)

    In a statistically stationary state, the velocity fluctuations should be of the order \(\sigma _{\delta {u_i}}^{(n)} \simeq {u_{\rm{c}}}\). In this case the Kalman filter corresponds to simple exponential smoothing:

    $$\alpha _i^{(n + 1)} \simeq {{\sigma _{\delta [{u_i}]}^{(n)}} \over {\sigma _{\delta {u_i}}^{(n)}}} \simeq {{2\pi \Delta {t^{(n)}}} \over {\sqrt 3 {T_c}}} \ll 1.$$

    The two filter parameters, Tc and uc, have to be chosen such that uc is roughly the integral velocity of turbulence if the flow enters a steady state and Tc is the characteristic time scale over which the flow evolves. In Cahuzac et al. (2011), LES of turbulence produced by the flow past a cylinder were shown to agree well with experimental data if Kalman filtering is applied with Tc and uc being set to the inverse of the expected vortex-shedding frequency and upstream velocity, respectively.

Figure 8:
figure 8

Slices of the mass density (left), squared vorticity (middle), and specific SGS turbulence energy (right) or AMR simulations of a gravitationally bound gas cloud in a wind after 3 Gyr of evolution. The box size is 4 Mpc. Image reproduced with permission from Schmidt et al. (2014), copyright by the authors.

A similar computational problem in astrophysics is an isothermal spherical cloud, which is bound by a static gravitational potential, in a homogeneous wind. The initial density profile of the cloud is determined by hydrostatic equilibrium. As mass is stripped from the cloud by the wind, a turbulent wake forms in the downstream direction. This problem was originally investigated as a simple model for the infall of small subcluster into the ICM of a much more massive cluster, computed in the frame of reference attached to the center of mass of the subcluster (Iapichino et al., 2008). LES with the shear-improved SGS model are presented in Schmidt et al. (2014), where the Kalman filter parameters are given by the velocity of the wind and the turn-over time of the largest eddies. Figure 8 compares the gas density and flow structure for two runs, one with the standard SGS model and the other one with the shear-improved model. In both cases, a 1283 root grid and three levels of refinement are used. The AMR control variables are the squared vorticity and the compression rate (Iapichino et al., 2008; Schmidt et al., 2009). The latter is defined by the substantial time derivative of the velocity divergence and tracks down the bow shock in front of the cloud. The top panels in Figure 8 show that the strain associated with the bow shock produces substantial SGS turbulence energy if the standard SGS model is applied. One can also discern the down-scaling of the SGS turbulence energy that is transported with the wind, as it enters the refined regions around the cloud. This is based on the algorithm explained in the previous section. Kelvin-Helmholtz instabilities between the low-density wind and the high-density gas in the cloud cause vortex shedding, which in turn produces turbulence. This becomes manifest in large values of ω and K in the turbulent wake that extends from the cloud towards the right. In contrast to the turbulent wake, however, the shear experienced by the gas when it is passing the bow shock is not associated with turbulence. As a result, the steep increase of K is largely a spurious effect induced by the eddy-viscosity closure (49). Indeed, the SGS turbulence energy in the shocked wind is substantially reduced if the shear-improved closure (118) is applied (bottom right panel in Figure 8). In comparison to the standard SGS model, production is also suppressed at the front side of the cloud and, hence, vortex shedding is the dominant source of turbulence production. This has a clearly visible impact on the structure of the turbulent wake.

In general, an appropriate prior choice of the filter parameters might not be obvious. Nevertheless, they can be calibrated a posteriori by performing low-resolution test runs. This is shown for cosmological simulations in Schmidt et al. (2014). In this case, the Kalman filter also has the merit of separating the turbulent flow in clusters from the gravity-driven bulk flows (see Section 6.3.3).

6 Astrophysical Applications

6.1 Thermonuclear combustion in white dwarfs

Among the various possible scenarios for supernovae of type Ia are thermonuclear explosions of gas accreting white dwarfs in close binary systems (Hillebrandt and Niemeyer, 2000; Röpke et al., 2011; Hillebrandt et al., 2013). If the mass of a white dwarf approaches the Chandrasekhar limit, explosive carbon and oxygen burning is ignited (Nonaka et al., 2012). Owing to the degeneracy of white dwarf matter, the thermonuclear reaction zones propagate as thin flame fronts, whose thickness δ f and propagation speed sf are determined by the very high thermal conductivity of the fuel and the nuclear reaction rates. This mode of burning is called deflagration. Since the burned material has lower density than the fuel, it rises because of its buoyancy. Consequently, the energy released by thermonuclear deflagration drives convection. Since eddies exert strong shear on rising bubbles of burning material, they are deformed into mushroom-like shapes and Kelvin-Helmholtz instabilities at the surface are rapidly producing turbulence (Malone et al., 2014). Eventually, this results in a very complex flame front with a fractal structure that cannot be resolved in numerical simulations over the full range of length scales.

This problem was addressed by performing LES with an effective propagation speed of the turbulent flame front (Reinecke et al., 2002; Schmidt and Niemeyer, 2006; Röpke et al., 2007, to mention just a few examples). In these simulations, the flame front is tracked by means of the level set method (Osher and Sethian, 1988; Reinecke et al., 1999), which is able to follow complex topological changes by determining the interface between fuel and burned material as the spatial surface for which a signed distance function is zero. The evolution of the distance function from given initial conditions is given by the advection velocity relative to the fluid plus the flame propagation speed. In a fully resolved simulation, the flame propagation speed would be given by the microscopic flame speed sf (also called laminar flame speed). If the flame front is underresolved, however, the wrinkling and folding of the front by turbulent eddies below the grid scale leads to an enhanced rate of energy release. In this case, the relevant time scale is not the microscopic diffusion time scale but the turn-over time of numerically unresolved eddies. Simple dimensional reasoning implies that sf has to be replaced by the turbulent flame speed \({s_{\rm{t}}} = \sqrt {2K} \) (Niemeyer and Hillebrandt, 1995; Niemeyer and Kerstein, 1997; Peters, 1999), where K is the SGS turbulence energy given by Eq. (64). In Schmidt et al. (2006c), the formula of Pocheau (1994),

$${s_{\rm{t}}} = {s_{\rm{f}}}\sqrt {1 + {C_t}{{2K} \over {s_{\rm{f}}^2}}},$$
(125)

is adopted for a smooth transition between laminar and turbulent flame propagation. A further complication comes from the pronounced anisotropy of turbulence at the flame front (Ciaraldi-Schoolmann et al., 2009; Malone et al., 2014). While the burned material inside the flame is highly turbulent, there is little or no turbulence in the fuel just outside the flame. In a way, this is similar to walls in terrestrial flows. For this reason, it is important to apply the dynamic procedure explained in Section 4.2 to locally calculate the eddy-viscosity coefficient, which determines the rate of production of K.

The minimal length scale for which the flame front is affected by turbulence is given by the Gibson scale G. If υ′ () is the mean turbulent velocity fluctuation on the length scale , then G is implicitly given by the condition

$$v^{\prime}({\ell _{\rm{G}}}) = {s_{\rm{f}}},$$
(126)

i.e., the turbulent velocity fluctuation on the Gibson scale equals the microscopic flame speed. For G, the turn-over time associated with υ′(G) is much longer than the crossing time of the flame front through an eddy of size . On these scales, the flame front is virtually unaffected by turbulence. If G, on the other hand, eddies will significantly distort the flame front on length scales between G and . In LES of a thermonuclear supernova, the Gibson scale G is much smaller than computationally feasible grid resolutions Δ during most of the deflagration phase and, consequently, a turbulent flame speed model has to be applied.Footnote 20

As an illustration of the level set method with the turbulent flame speed (125), Figure 9 shows snapshots of the flame front for a thermonuclear supernova simulation from Schmidt and Niemeyer (2006). The turbulent velocity fluctuations on subgrid scales, \(\sqrt {2K} \), are shown as color shades at the flame surface. The asymptotic value of turbulent flame speed is \({s_{\rm{t}}} \simeq {C_{\rm{t}}}\sqrt {2K} \) if stsf. In this simulation, a Poisson process is used to randomly place small ignition spots in the high-density core of the white dwarf. The statistics of this process is based on a simple model for temperature fluctuations produced by convection in pre-ignition phase. Although the number of ignition resulting from this model appears to be by far too large in the light of recent numerical studies of the ignition process (Nonaka et al., 2012), the simulation nevertheless demonstrates in an exemplary manner how the turbulent deflagration progresses. At early time (a), one can see a large number of small bubbles generated by the stochastic ignition process at distances of the order 100 km from the center of the white dwarf. As the burning bubbles are rising from the centre, they begin to form the typical Rayleigh-Taylor mushroom shapes (b). At this point, the effective flame propagation speed is already dominated by turbulence. After about half a second (c), the turbulent flames mostly have merged into a single structure of about 1000 km diameter. Then the burning front rapidly expands to much larger radii and causes the white dwarf to explode after roughly one second (d).

Since the the asymptotic rising velocity of a perturbation of size due to its buoyancy is given by the Sharp-Wheeler relation (Sharp, 1984)

$${v_{{\rm{RT}}}}(\ell) = 0.5\sqrt {\ell {g_{e{\rm{ff}}}}},$$
(127)

where \({g_{{\rm{eff}}}}\) is the effective gravity associated with the density contrast at the interface between burned and unburned material, \({s_{\rm{t}}} \simeq {v_{{\rm{RT}}}}(\Delta)\) was used as a simple turbulent flame speed model (e.g., Gamezo et al., 2003). In Schmidt et al. (2006c), buoyancy is included as an ad-hoc source term in the SGS turbulence energy equation (27):

$$\Gamma \propto \rho {g_{{\rm{eff}}}}\sqrt {2K}.$$
(128)

Here, Γ is defined such that it vanishes everywhere except for the close vicinity of the flame front. Moreover, Γ is assumed to be non-zero only if the so-called fire polishing length \({\lambda _{{\rm{fp}}}} = 4\pi s_{\rm{f}}^2/{g_{{\rm{eff}}}}\) is smaller than the grid resolution Δ. The fire polishing length is the smallest length scale on which perturbations in the flame front are Rayleigh-Taylor unstable. If the buoyancy term dominates over the turbulent cascade, \(\sqrt {2K} \sim {v_{{\rm{RT}}}}(\Delta)\) is obtained as asymptotic solution. However, different numerical studies indicate Kolmogorov scaling on small scales (Zingale et al., 2005; Röpke et al., 2007; Ciaraldi-Schoolmann et al., 2009), corresponding to Δ ≪ Σ in Eq. (27). Although the impact of the turbulent flame speed model on the energy release was demonstrated in various simulations (Niemeyer and Hillebrandt, 1995; Schmidt et al., 2006c), it tends to become less significant if the burning is resolved down to very small scales. This has become possible with the application of AMR and the power of contemporary computing facilities (Ma et al., 2013).

However, since the predictions from pure deflagration models are not consistent with observational properties of the majority of type Ia supernovae, a transition from the deflagration phase to a supersonic detonation is now considered as the most likely explosion mechanism of Chandrasekharmass white dwarfs (Röpke et al., 2011; Hillebrandt et al., 2013). A theoretical explanation of such a transition is a long standing problem. One possibility is that very strong local velocity fluctuations at the onset of distributed burning might trigger the transition to a detonation (Niemeyer and Kerstein, 1997; Khokhlov et al., 1997; Lisewski et al., 2000; Woosley, 2007). The statistical distribution of turbulent velocity fluctuations following from the SGS turbulence energy K in a high-resolution simulation (see Figure 10) shows that such events might indeed occur (Röpke, 2007). The very wide tail indicates the strong intermittency of turbulence at the flame front. On the basis of this result, a criterion was recently proposed to set off detonations in LES of thermonuclear supernovae (Ciaraldi-Schoolmann et al., 2013).

Figure 9:
figure 9

Evolution of the flame front in a thermonuclear supernova simulation. The colour scale indicates the logarithm of the magnitude of turbulent velocity fluctuation predicted by the SGS model in units of cm/s (i.e., 6 and 8 corresponds to 10 and 1000 km/s, respectively). A similar simulation is presented in Schmidt and Niemeyer (2006). The spatial scale is adjusted to the bulk expansion of the white dwarf by means of a co-moving grid technique (Röpke and Hillebrandt, 2005b).

Figure 10:
figure 10

Flame front in a high-resolution simulation of a thermonuclear supernova. High turbulent velocity fluctuations occur in the reddish regions. Extremely strong fluctuations that could trigger a transition to a detonation are indicated by the green arrows. Image reproduced with permission from Röpke (2007), copyright by AAS.

6.2 Galaxy simulations

Numerical simulations of galaxies, particularly from cosmological initial conditions, cannot fully resolve processes such as star formation and feedback from supernova (SN) explosions. As suggested in Joung et al. (2009), Scannapieco and Brüggen (2010), and Braun et al. (2014), a production term Σ* due to feedback can be incorporated as an additional source term in the SGS turbulence energy equation (27). The simplest model for this term is

$${\Sigma _ \star } = {C_ \star }{{\rho {e_{{\rm{SN}}}}} \over {{\tau _{{\rm{ff}}}}}},$$
(129)

where eSN is the typical energy released by a core-collapse supernova and τffρ−1/2 is the free-fall time scale. The parameter C * controls the effective time scale of the feedback. The above expression for Σ* follows from the simple Kennicutt-Schmidt relation \({\dot \rho _\star} \propto {\rho ^{ - 3/2}}\) for the star formation rate, where the factor C*ρ/τff is some fraction of \({\dot \rho _ \star }\). Apart from stellar feedback, turbulence is produced by instabilities in the ISM, such as gravitational, thermal, and hydrodynamical instabilities. Typically, energy is injected on numerically resolved length scales by these instabilities and transferred to smaller scales by the turbulent cascade. If we simply assume that a turbulent velocity of magnitude V is produced on the length scale L, the energy flux through the turbulent cascade is of the order Σ ∼ ρV3/L (Braun and Schmidt, 2012). In LES, Σ is defined in terms of the grid scale and the shear of the numerically resolved flow via the closure (55). As a very crude model, let us consider the local equilibrium between production and dissipation, i.e., Σ + Σ*ρε. By substituting Eqs. (55), (129), and ερK3/2/Δ, the equilibrium condition reads

$${C_1}\Delta {(2K)^{1/2}}\vert S\ast{\vert ^2} - 4{C_2}K{{{u_{i,k}}{u_{j,k}}S_{ij}^\ast} \over {\vert \nabla \otimes u{\vert ^2}}} - {2 \over 3}Kd + {C_ \star }{{{e_{{\rm{SN}}}}} \over {{\tau _{e{\rm{ff}}}}}}\sim{{{K^{3/2}}} \over \Delta }.$$

. This relation implies the turbulent pressure

$${P_{{\rm{eq}},{\rm{tot}}}}\sim\rho {\Delta ^{2/3}}{\left({{K \over \tau } + {C_ \star }{{{e_{{\rm{SN}}}}} \over {{\tau _{{\rm{ff}}}}}}} \right)^{2/3}},$$
(130)

where

$$\tau = {{\rho K} \over \Sigma },$$
(131)

is the dynamical time scale associated with energy transfer through the cascade. Depending on the energy ratio K/eSN and the ratio of the dynamical and feedback time scales, τ/τff, the production of SGS turbulence energy is dominated by the turbulent cascade or by supernovae.

The star formation rate, which in turn determines the feedback rate, can be parameterized by the gas density, temperature, and turbulence intensity (Krumholz and McKee, 2005; Padoan and Nordlund, 2011; Hennebelle and Chabrier, 2011; Federrath and Klessen, 2012; Padoan et al., 2012). The two main parameters appearing in these parameterizations are the turbulent Mach number \({\mathcal{M}_ \star }\) and the virial parameter α* (Bertoldi and McKee, 1992), which correspond to the ratios of the turbulence energy to the internal and gravitational energies:

$${\mathcal{M}_ \star } = \sqrt 3 {{\sigma ({\ell _ \star })} \over {{c_{\rm{s}}}}}\quad {\rm{and}}\quad {\alpha _ \star } = {{15{\sigma ^2}({\ell _ \star })} \over {\pi G\rho \ell _ \star ^2}}.$$
(132)

Here, σ(ℓ) is the turbulent velocity dispersion on the length scale , cs the speed of sound, and G Newton’s constant. A major problem in galaxy simulations is the choice of a characteristic scale ℓ* of star formation. In Braun and Schmidt (2012) it is argued that * is given by the Jeans length for gravitational instability, i.e.,

$${\ell _ \star } = {\lambda _{\rm{J}}}: = {c_{\rm{s}}}{\left({{\pi \over {\gamma G\rho }}} \right)^{1/2}},$$
(133)

where γ is the adiabatic exponent of the gas. Since the SGS model predicts the specific turbulence energy K = 3σ2(Δ)/2 on the grid scale Δ, it is possible to estimate the turbulent velocity dispersion in gravitationally unstable, star-forming clouds as

$$3\sigma _ \star ^2 = 2K{\left({{{{l_ \star }} \over \Delta }} \right)^{2\eta }},$$
(134)

where the scaling exponent is in the range between 1/3 and 1/2, depending on the compressibility of the gas and the intermittency of turbulence (Kritsuk et al., 2007; Schmidt et al., 2008, 2009; Hennebelle and Falgarone, 2012). The star formation rate is then given by

$${\dot \rho _ \star } = {\rm{SFR}}({\alpha _ \star },\;{\mathcal{M}_ \star }){\rho \over {{\tau _{{\rm{ff}}}}}}.$$
(135)

The mass fraction SFR(α*, \({\mathcal{M}_ \star }\)) that is converted into stellar mass per free-fall time is model-dependent. Basically, this factor accounts for the gravo-turbulent fragmentation of star-forming clouds.Footnote 21 For example, a power-law function based on data from numerical simulations is proposed in Krumholz and McKee (2005). In Padoan and Nordlund (2011) it is assumed that SFR \(({\alpha _\star},\;{\mathcal{M}_\star})\) is determined by the log-normal probability density function, whose width depends on the turbulent Mach number \({\mathcal{M}_\star}\) (Padoan et al., 2007; Federrath et al., 2010b), and a critical density that is controlled by the virial parameter α * . The different parameterizations are compared and calibrated by numerical data in Federrath and Klessen (2012). A simple star formation law, for which SFR (α*) is solely determined by the virial parameter, is proposed in Padoan et al. (2012).

A star formation rate proportional to the density of molecular hydrogen, \({f_{{{\rm{H}}_2}}}\rho \) instead of the total gas density ρ, was suggested on grounds of the observed tight correlation between the star formation rate and the column density of molecular hydrogen (Krumholz et al., 2009; Gnedin et al., 2009). Since molecular hydrogen forms only in the cold phase of the interstellar medium, the star formation law (135) is further modified by replacing ρ and cs with the mean mass density and speed of sound, respectively, in the cold phase (Braun and Schmidt, 2012). Since the separation into cold and warm phases cannot be fully resolved in galaxy simulations, a model such as (Springel and Hernquist, 2003) has to be employed to estimate the fractional densities of the two phases. This entails additional turbulence energy production by cooling instabilities on length scales below the grid resolution (Braun and Schmidt, 2012; Iwasaki and Inutsuka, 2014). If the resulting internal driving due to feedback and cooling instabilities, Σint = Σ* + ΣTI, along with the closure (55) for the compressible turbulence energy flux Σ are incorporated into Eq. (64), a full model for the numerically unresolved turbulence energy budget in galaxy simulations is obtained (Braun et al., 2014).

Figure 11:
figure 11

The SGS turbulence energy ρK has several functions in the multiphase model for the turbulent ISM in disk galaxies (Braun et al., 2014). It determines the shape and width of the mass density PDF of cold clumps and the effective pressure equilibrium between the cold and warm phases. Via the density PDF, it influences the star formation rate, which acts back on ρK through feedback. Other sources of ρK are the turbulent cascade and cooling instabilities.

In combination with the two-phase model for the gas and a particular flavor of the star formation and feedback models outlined above, this SGS turbulence energy model has recently been applied to adaptively refined LES of isolated disk galaxies. As schematically shown in Figure 11, the computation of the SGS turbulence energy ρK plays a central role. In Figure 12 (left), an equatorial slice of the gas density illustrates the fragmentation of the inner gaseous disk into dense star-forming clumps (bluish regions) in one of these simulations. The dilute gas with densities below 1 cm−3 (reddish regions) is partially heated by supernovae. The fragmentation of the disk produces turbulence, as indicated by the slice of ρK in Figure 12 (right). Particularly large turbulent energies are found in regions with strong feedback, which suggest that the non-thermal feedback from supernovae plays a significant role in the production of SGS turbulence energy. This is indeed confirmed by the statistics of the production terms plotted in Figure 13 (left). This plot shows profiles of

$$\frac{{\Sigma _{\operatorname{int} } }} {{\Sigma _{tot} }} = \frac{{\Sigma _* + \Sigma _{TI} }} {{\Sigma + \Sigma _* + \Sigma _{TI} }} and \frac{\Sigma } {{\Sigma _{tot} }} = \frac{\Sigma } {{\Sigma + \Sigma _* + \Sigma _{TI} }},$$

which are calculated by averaging Σ, Σint, and Σtot over bins of SGS turbulence energy per unit mass. The maximum of Σ/Σtot around K ∼ 100 (km/s)2 implies that the turbulent cascade maintains a ground level corresponding to a turbulent velocity dispersion around 10 km/s, while the internal driving caused by supernovae excites much stronger turbulent velocity fluctuations. Remarkably, the average turbulence energy is quite close to the equilibrium value implied by the balance between production and dissipation. This can be seen in Figure 13 (right), where the ratio of the equilibrium pressure to the dynamical turbulent pressure,

$${{{P_{{\rm{eq}},{\rm{tot}}}}} \over {{P_{{\rm{sgs}}}}}}: = {1 \over K}{\left({{{\Delta {\Sigma _{{\rm{tot}}}}} \over \rho }} \right)^{2/3}},$$
(136)

is plotted against Psgs = 2/3ρK. Also shown is the asymptotic equilibrium pressure

$$\frac{{P_{eq,\operatorname{int} } }} {{P_{sgs} }}: = \frac{1} {K}\left( {\frac{{\Delta \Sigma _{\operatorname{int} } }} {\rho }} \right)^{2/3} $$
(137)

if Σtot ≃ Σint ≫ Σ, i.e., internal driving dominates. This is the case for high turbulence intensity. Since Peq,int/Psgs is only slightly above unity, turbulence is close to equilibrium in this regime. For lower SGS turbulence energies, Eq. (130) tends to overestimate the turbulent pressure, which indicates that turbulence production exceeds dissipation. In this case, the turbulent cascade contributes significantly to the production. For comparison, also profiles of the thermal pressures of the cold and warm phases are plotted. On the one hand, Peq,tot is small compared to the warm-phase pressure, except for the highly turbulent regions with strong feedback, where PwPeq,totPsgs. On the other hand, Peq,totPc. The pressure equilibrium between the cold and warm phases, which is one of the basic assumptions of the model, is therefore dominated by the turbulent pressure in the cold clumps and the thermal pressure in the warm medium.

Figure 12:
figure 12

Slices of the number density (left) and SGS turbulence energy (right) on the equatorial plane in an isolated disk galaxy simulation (Braun et al., 2014).

6.3 Cosmological simulations

Since cosmological scale structure formation produces a strongly clumped medium through gravitational contraction or collapse, the numerical simulation of turbulence on cosmological scales is particularly challenging. Potential production mechanisms of turbulent flows in the baryonic gas are mergers between dark-matter halos and the accretion of gas into halos, but also feedback from active galactic nuclei and winds produced by strongly star-forming galaxies (Dolag et al., 2008; Kravtsov and Borgani, 2012). A largely open question concerns whether the gaseous component of halos is in a state of developed turbulence. While there is no doubt about turbulence in the interior of galaxies, there is no direct observational evidence yet for turbulence on larger scales, particularly in the ICM. Theoretical and numerical studies suggest that magnetic fields play a key role in the physical dissipation mechanism and, possibly, the onset of instabilities, but the associated length scales are highly uncertain (Ferrari et al., 2008; Brüggen, 2013). Notwithstanding these uncertainties, turbulent flows resulting from cosmological structure formation are numerically investigated by computing the evolution of an ideal fluid subject to the gravitational potential of dark matter, which is modeled as a collisionless N-body system (Borgani and Kravtsov, 2011). To follow gravitational collapse, AMR is essential if Eulerian grid codes are used. As pointed out in Section 5, turbulence is generally underresolved in a clumpy medium, even at very high refinement levels.

Figure 13:
figure 13

Left: profiles of the rate of production through the turbulent cascade, Σ, and internal driving due to numerically unresolved feedback from supernovae and thermal instabilities, Σint, relative to Σtot = Σ + Σint for the disk galaxy simulation shown in Figure 12 (Braun et al., 2014). Right: profiles of the pressure ratios PX/Psgs, where Psgs = 2/3ρK is the numerically unresolved turbulent pressure, for the thermal pressures Pc and Pw of the cold and warm phases and the equilibrium pressures Peq,tot and Peq,int defined by Eqs. (136) and (137), respectively. By courtesy of Harald Braun.

This problem was addressed for the first time by combining LES and AMR in Maier et al. (2009). The main idea is to solve the filtered Euler equations (4042) and the SGS turbulence energy equation (44) in co-moving coordinates. For grid refinement and de-refinement, Maier et al. (2009) apply the power-law relation (108) between the SGS turbulence energies at different refinement levels. In contrast to the algorithms outlined in Section 5.1, however, their implementation is not conservative by construction. The SGS turbulence energy decrement \( - \Delta (\overline {\rho K})\) in the case of refinement from a coarser level, for example, is simply compensated by adjusting the interpolated momenta \((\rho U)_n^{\ast}\) such that

$${1 \over N}\sum\limits_n {{1 \over 2}{{(\rho U)_n^2} \over {{\rho _n}}}} = {1 \over N}\sum\limits_n {{1 \over 2}{{{{[(\rho U)_n^\ast]}^2}} \over {{\rho _n}}}} + \Delta (\overline {\rho K}),$$
(138)

which follows by setting

$${(\rho K)_n} = {r^{ - 2\eta }}(\rho K)_n^\ast\quad \;\;{\rm{and}}\quad \;\;{(\rho U)_n} = (\rho U)_n^\ast\sqrt {1 + {{2{\rho _n}(\rho K)_n^\ast} \over {{{[(\rho U)_n^\ast]}^2}}}(1 - {r^{ - 2\eta }})}.$$

. While the first relation is identical to Eq. (116), the second relation, when substituted into the total energy balance, implies

$$\Delta (\overline {\rho e}) = {1 \over 2}{{{{[(\rho U)_n^\ast]}^2}} \over {{\rho _n}}} - {1 \over 2}{{(\rho U)_{{\rm{crs}}}^2} \over {{\rho _{{\rm{crs}}}}}}.$$
(139)

This is just the standard method of compensating the difference of the kinetic energies between refinement levels entirely by internal energy, which is applied in addition to the kinetic energy transfer (138). Analogous relations are applied to average the data from a finer to a coarser level. The rationale of Eq. (138) is that the resolved turbulent velocity fluctuations should increase from a coarser to a finer level. However, as argued in Section 5.1 and in Schmidt et al. (2014), conservative interpolation of the momenta to a finer level entails already an increase of the resolved kinetic energy that encompasses and, in some cases, even overestimates the scale-dependence of the numerically resolved turbulence energy. Both energy and momentum conservation is guaranteed by compensating this increase with \(\Delta (\overline {\rho K})\) and \(\Delta (\overline {\rho e})\) defined by Eqs. (114) and (115), respectively.

Despite its shortcomings, the method based on Maier et al. (2009) has its merits as a first approximation. For example, Figure 14 shows the density of the baryonic gas and the small-scale turbulence in an adiabatic cosmological simulation of a cluster from Maier et al. (2009). As indicator of turbulence, the magnitude of the numerically unresolved turbulent velocity fluctuation \(\sqrt K \) is scaled down from refinement level l to the minimal cell size at the maximal level lmax via the power law factor \({r^{(l - {l_{\max }})\eta }}\). It is particularly interesting that the infall of a subhalo, which is marked by the small box in slice (a), produces a pronounced turbulent wake, quite similar to the idealized scenario discussed in Section 5.2. However, a shear-improved model was not used in this simulation. Even at z = 0, the turbulent velocity slice (d) in Figure 14 shows a clear trace of the minor merger, which is not discernible in the density slice (c).

6.3.1 Turbulence production and support against gravity

In adiabatic simulations of cosmological structure formation with Enzo (Bryan et al., 2014), statistics of the SGS turbulence energy indicate different production mechanisms of turbulence in the ICM and the warm-hot intergalactic medium (WHIM) (Iapichino et al., 2011). While the former is found at high densities and temperatures and is dominated by weakly compressible, subsonic turbulence, the WHIM is constituted by gas of lower density, which is undergoing compression by shocks. The plot in Figure 15 shows different histories of the mean thermal and SGS turbulence energies in the ICM and WHIM. While the SGS turbulence energy in the ICM peaks around redshift z = 1, which can be interpreted as a consequence of mergers between galaxy clusters, the energy in the WHIM gradually rises as a result of the continuous turbulence production by accretion shocks. These are caused by the infall of low-density gas into the potential wells of clusters and filaments, which accelerates the gas to supersonic speed.

Figure 14:
figure 14

Slices of baryonic overdensity (left) and the turbulent velocity at the scale of the highest-resolution level in km/s (right) in a cluster simulation at redshifts z = 0.05 (top) and z = 0.0 (bottom). The comoving size of the shown region is 6.4 Mpc h−1. Image reproduced with permission from Maier et al. (2009), copyright by AAS.

In addition to other indicators, Iapichino et al. (2011) analyses the relative importance of the turbulent and thermal pressures for the support of the gas against gravity. As shown in Schmidt et al. (2013), the local support of self-gravitating gas due to the thermal pressure P is given by

$${\Lambda _{{\rm{therm}}}} = - {1 \over \rho }{\Delta ^2}P + {1 \over {{\rho ^2}}}\nabla \rho \cdot \nabla P.$$
(140)

The above expression is derived by taking the divergence of the equation for the gas velocity, which results in an equation for the rate of compression of a fluid parcel, −Dd/Dt, where D/Dt signifies the substantial time derivative and d the divergence of the flow. In the weakly compressible limit, the thermal support reduces to Λtherm ≃ −(∇2P)/ρ. By considering Eq. (23), it immediately follows that the influence of turbulence below the grid scale can be expressed as

$${\Lambda _{{\rm{sgs}}}} = - {1 \over \rho }{\nabla ^2}{P_{{\rm{sgs}}}} + {1 \over {{\rho ^2}}}\nabla \rho \cdot \nabla {P_{{\rm{sgs}}}},$$
(141)

where Psys = 2/3ρK is the turbulent pressure on the grid scale. However, Λsgs does not account for the effect of the non-diagonal stresses \(\tau _{ij}^{\ast}\). The contribution of numerically resolved turbulence, on the other hand, is given by

$${\Lambda _{{\rm{turb}}}} = {1 \over 2}({\omega ^2} - \vert S{\vert ^2}),$$
(142)

where ω = Λ × u is the vorticity and ∣S∣ the rate of strain (see equation 15). The first term on the right-hand side is associated with the pressure-like support caused by turbulent eddies, the second term with the compression of the gas by shocks. Figure 16 shows a mass-weighted histogram of the ratio

$${r_{{\rm{tp}}}} = {{\rho ({\omega ^2} - \vert S{\vert ^2}) - 2{\nabla ^2}{P_{{\rm{sgs}}}}} \over {2{\nabla ^2}P}}$$
(143)

for different overdensities of the baryonic gas. The ratio rtp is an approximation to the ratio of the turbulent and thermal support functions, (Λturb + Λsgs)/Λtherm, which was used for the sake of comparability with Zhu et al. (2010). The distribution of rtp plotted in Figure 16 shows two peaks, one at intermediate densities and the other at higher densities, which can be associated with the WHIM and the ICM. There is a trend of decreasing rtp toward high densities. This reflects the smaller Mach numbers of turbulence in the ICM, but a caveat is that only positive contributions to the support are taken into account. Since it is demonstrated in Schmidt et al. (2013) and Latif et al. (2013b) that Λturb can be predominantly negative in the presence of shocks, the question of the turbulent support of the gas in clusters needs to be revisited.

Figure 15:
figure 15

Evolution of the mean internal (solid) and SGS turbulence (dashed) energies in the ICM and the warm-hot intergalactic medium (WHIM) of a galaxy cluster. Image reproduced with permission from Iapichino et al. (2011), copyright by the authors.

6.3.2 Gravitational collapse of gas in primordial halos

The Enzo implementation of the SGS model for cosmological fluid dynamics was also applied to the gravitational collapse of gas clouds in primordial atomic cooling halos. This scenario is potentially relevant for the formation of intermediate-mass black holes through direct collapse and subsequent accretion (Latif et al., 2013a, c). The resulting black holes are possibly the seeds for supermassive black holes at the centers of galaxies. Direct collapse into a supermassive star, which subsequently collapses into an intermediate-mass black hole, can occur if the molecular hydrogen formation is suppressed through photodissociation by a Lyman-Werner radiation background. To follow the collapse of the primordial gas down to AU scales, deep-zoom in simulations with many level of refinements have to be performed. In Latif et al. (2013c), both LES and ILES of several different halos of masses of the order 107 M are presented. One of the most remarkable results of this study is that the additional turbulent viscosity produced by the SGS model favors disk-like structures around collapsed objects in LES. This can be seen from the comparison of density slices in the innermost regions around the density peaks. Figure 17 shows that more or less compact collapsed structures are produced in ILES. In contrast, more extended and disk-like structures are found in the corresponding LES, as shown in Figure 18.

Figure 16:
figure 16

Mass-weighted histogram of the turbulent-to-thermal support ratio vs. baryon overdensity (the solid line indicates the mean value for given overdensity). Image reproduced with permission from Iapichino et al. (2011), copyright by the authors.

Apart form morphological differences, the SGS model also has a significant influence on the accretion of mass by the protostars that are formed through the collapse of gas clouds in the halos. As mentioned in Section 3.4, it is possible to follow the evolution of gravitationally bound dense objects, such as protostars, by inserting sink particles (Federrath et al., 2010a; Wang et al., 2010). This method is applied in Latif et al. (2013e) to follow the accretion history of protostars in atomic cooling halos. The resulting time evolution of the accretion rate is plotted in Figure 19 for simulations of three different halos, using both ILES and LES. As one can see, the accretion rates reach peak values around 10 M yr−1 roughly within 104 yr. This value agrees with the theoretical expectation for Bondi-Hoyle accretion. A comparison of ILES and LES suggests a systematically higher accretion rate for LES. This trend is confirmed by calculating the cumulative masses of the sink particles (see right plot in Figure 19), which reach masses above 105 M. As a result, the SGS model favors the formation of higher black hole masses.

6.3.3 Turbulent velocity dispersion

A key question that arises in connection with analyzing turbulence production by cosmological structure formation concerns the turbulent velocity dispersion σturb. In Iapichino et al. (2013), σturb is defined as \({\sigma _{{\rm{turb}}}} = \sqrt {2K} \), which implies that σturb is identified with the magnitude of the turbulent velocity fluctuations on the grid scale of the simulation. Although this variable obviously depends on numerical resolution, it is nevertheless useful to infer the dependence on various factors that influence the production of turbulence. For example, Figure 20 (left plot) shows the ratio of the turbulent and thermal pressures, where Pt = Psgs = 2/3ρK, and the Doppler broadening parameter \({b_{\rm{t}}}\) (Evoli and Ferrara, 2011) for a cosmological simulation with radiative background and cooling in a box of 10 Mpc h−1 comoving size. The analysis is carried out for data cubes at redshift z = 2.0. For the intergalactic medium (IGM), which is usually defined by moderate baryonic overdensities ρ/ρ0 in the range from 1 to about 50, bt is found to increase with density from roughly 1 to 10 km/s. The WHIM, on the other hand, has a flat turbulent velocity dispersion, corresponding to Pt/P ∼ 0.1. The phase plot in Figure 20 shows that the ratio bt/b varies over several orders of magnitude for different densities and temperatures and reaches peak values ∼ 1. The WHIM is associated with gas that is heated by accretion shocks to temperatures between 105 and 107 K. For this reason, the WHIM tends to be more turbulent then the diffuse gas in the IGM. The phase diagram also shows that the distinction between gas in the IGM (1 ≤ ρ/ρ0 ≤ 50) and the WHIM (105 K ≤ T ≤ 107 K) is not mutually exclusive and somewhat arbitrary.

Figure 17:
figure 17

Gas density in the central 300 AU of collapsing atomic cooling halos computed with ILES. Image reproduced with permission from Latif et al. (2013c), copyright by the authors.

Figure 18:
figure 18

The same halos as in Figure 17 computed with an explicit SGS model. Image reproduced with permission from Latif et al. (2013c), copyright by the authors.

Figure 19:
figure 19

Left plot: accretion rates of the most massive sink particles in simulations of three different halos. Right plot: Comparison of the mass distributions of the sink particles in LES and ILES. Image reproduced with permission from Latif et al. (2013e), copyright by the authors.

A resolution-dependent turbulent velocity dispersion is avoided by the — hybrid model for Rayleigh-Taylor-driven turbulence (see Section 3.4), however, at the cost of smearing out turbulent fluctuations on numerically resolved length scales. This model was used to simulate the production of turbulence by AGN feedback in galaxy clusters (Scannapieco and Brüggen, 2008; Brüggen and Scannapieco, 2009). Another caveat of both the K-L model and the SGS model based on Maier et al. (2009) are the constant closure coefficients. Strictly speaking, constant-coefficient models are applicable only to statistically homogeneous and stationary turbulence. Turbulence produced by cosmological structure formation and AGNs, however, is highly inhomogeneous and non-stationary.

To ameliorate this problem as well as the dependence of σturb on the grid scale, the shear-improved SGS model outlined in Section 5.2 has recently been applied to cosmological structure formation (Schmidt et al., 2014). By running simulations of the Santa Barbara cluster (Heitmann et al., 2005) with the Nyx code (Almgren et al., 2013), it is demonstrated that the application of the standard SGS model for homogeneous turbulence can produce biases, particularly in regions of active turbulence production, such as the WHIM. In Figure 21, the growth of the cluster is illustrated by slices at different redshifts from a simulation performed with the shear-improved SGS model. There clearly is a correlation between the high vorticity and SGS turbulence energy in the cluster and a sharp drop at the outer accretions shocks. Although a meaningful comparison to the results in Iapichino et al. (2013) cannot be made at this point because the Santa Barbara cluster is based on a matter-dominated universe without cosmological constant and the gas dynamics is adiabatic, some general conclusions regarding the nature of turbulence in clusters can be drawn. In Schmidt et al. (2014), the turbulent velocity dispersion is defined by

$${\sigma _{{\rm{turb}}}} = \sqrt {{U^{\prime 2}} + 2K},$$
(144)

where U is the fluctuating component of the velocity computed with the Kalman filter (see Section 5.2). The average values of σturb for logarithmic bins of the baryonic overdensity as well as the corresponding radial profiles are plotted for different numerical resolutions in Figure 22. Except for the lowest-resolution case, the radial profiles of σturb show little sensitivity to numerical resolution. This is an important property of the turbulent velocity dispersion defined by Eq. (144). While the profiles are nearly flat for the ICM, there is a sharp drop around 10 Mpc, which is about the radius of the outer accretion shocks in Figure 21. Since the gravitational pull of the cluster causes a large non-turbulent bulk flow toward the center, the total flow velocity beyond the accretion shocks is much higher than σturb. As a function of the overdensity ρ/ρ0, σturb gradually increases towards the center of the cluster. Although the simulations in Iapichino et al. (2013) differ in important aspects, a roughly similar trend can be seen for the Doppler broadening parameter bt in Figure 20.

Figure 20:
figure 20

Left plot: Average ratios of turbulent and thermal pressures vs. baryonic overdensity in a 10 Mpc h−1 box with heating and cooling at z = 2.0. Right plot: Ratio of turbulent to thermal Doppler broadening for bins of temperature and baryonic overdensity. Image reproduced with permission from Iapichino et al. (2013), copyright by the authors.

The turbulent kinetic energy \({1 \over 2}\rho \sigma _{{\rm{turb}}}^2\) and the energy flux S through the turbulent cascade (see Section 2.2) define the dynamical time scale

$$\tau = {{\rho \sigma _{{\rm{turb}}}^2} \over {2\Sigma }}.$$
(145)

The radial profiles of 1/τ plotted in Figure 23 show that the dynamical time scale in the ICM is several Gyr. The pronounced peak at 10 Mpc radius is a further indication of turbulence production by the accretions shocks. Analogous to Eq. (145), the dissipation time scale can be defined as

$${\tau _\epsilon } = {{\sigma _{{\rm{turb}}}^2} \over {2\epsilon }}.$$
(146)

For statistically stationary turbulence, the balance between turbulence production and dissipation, Σ ∼ ρε, implies τ ∼ τ ε . This can indeed be seen in Figure 23 for the central region of the cluster. In the vicinity of the accretion shocks and outside the cluster, however, ετ ε . In this case, the flow is far from equilibrium. These results demonstrate how (nearly) scale-invariant quantities computed with the shear-improved SGS model can be utilized to investigate statistical properties of turbulence.

Figure 21:
figure 21

Slices of the baryonic mass density (left), vorticity modulus (middle), and specific SGS turbulence energy (right) at different redshifts for a simulation of the Santa Barbara cluster (Schmidt et al., 2014). The box size is 64 Mpc h−1. Image reproduced with permission from Schmidt et al. (2014), copyright by the authors.

Figure 22:
figure 22

Mean turbulent velocity dispersion defined by Eq. (144) vs. baryonic overdensity (left) and radius from the center (right) for simulations of the Santa Barbara cluster with different numerical resolutions. Starting from a uniform root-grid with N0 grid cells, refinement by overdensity and vorticity modulus up to lmax levels is applied. Image reproduced with permission from Schmidt et al. (2014), copyright by the authors.

Figure 23:
figure 23

Radial profiles of the inverse dynamical (left) and dissipation (right) time scales for the same simulations as in Figure 22. Image reproduced with permission from Schmidt et al. (2014), copyright by the authors.

6.4 Concluding remarks

It is a particular difficulty of validating LES in astrophysics that neither DNS nor sufficiently accurate experimental data exist. For virtually all instances of astrophysical turbulence, DNS are infeasible because of the very high Reynolds numbers (maybe, with the exception of turbulence in the ICM) and observations are usually not sensitive enough to clearly discriminate between different numerical models. As a consequence, the standards for model validation in astrophysics are very different compared to engineering or atmospheric sciences, where direct comparisons with measurements or DNS data are commonly made.

Statistical properties of homogeneous turbulence such as two-point statistics, turbulence energy spectra, or the distribution of mass density fluctuations are not very sensitive on the application of an explicit SGS model, provided that the numerical resolution is high enough. This is necessarily so, because the scale locality of hydrodynamical turbulence implies that properties related to inertialrange dynamics on some scale must be largely independent of effects on much smaller scales. However, the determination of second-order scaling exponents of compressible turbulence might benefit from the application of an explicit SGS model by reducing the bottleneck effect in the energy spectra (Woodward et al., 2006). Moreover, the SGS turbulence energy has its merits as a buffer variable, which helps to reduce artificial manipulations of the internal energy in AMR simulations (Schmidt et al., 2014). Much more difficult, however, is the assessment of systematic differences in the flow structure compared to ILES. An example are the preferentially disk-like structures in LES of collapsing primordial halos (Latif et al., 2013c). Although the correct solution is unknown, we know for sure that there are couplings between numerically resolved scales and subgrid scales in simulations, which are caused by the turbulent stresses on the grid scale. The calculation of these stresses on the basis of scale separation (Section 2) and testable closures (Sections 3 and 4) is at least physically better motivated and provides a more accurate approximation than purely numerical truncation errors.

Most likely, however, the main utility of SGS models in astrophysics is going to be the treatment of complex sub-resolution physics. The calculation of the turbulent burning speed in supernova simulations, which was pioneered by Niemeyer and Hillebrandt (1995), is just one example. Recently, Braun et al. (2014) demonstrated that feedback from supernovae in galaxy simulations can be implemented as a source term in the SGS turbulence energy equation in addition to heating. This offers an alternative to the often employed kinetic feedback, which produces resolved gas motions. Particularly in cosmological simulations, however, the production of turbulence by supernova blast waves is mostly an unresolved process. The application of LES in cosmological simulations therefore not only improves hydrodynamics and the consistency of AMR, but it allows for a treatment of feedback and mixing processes below the local grid scale.

A potential application of LES, which has not been exploited yet in astrophysics, are three-dimensional simulations of convection in stars, including the convection zone of our Sun. Arnett et al. (2014) investigate the performance of RANS in comparison to ILES for the cores of massive stars. Their model, which encompasses nucleosynthesis, reproduces statistics from three-dimensional simulations with high resolution. LES would be an intermediate approach that captures part of the turbulent convective flow, while explicitly modeling mixing processes on smaller scales. Many problems in astrophysics also involve magnetic fields. In this case, there is an analogy to feedback, because small scales can significantly contribute to the amplification of magnetic fields in turbulent plasmas. As demonstrated, for example, by Sur et al. (2010) and Latif et al. (2013d), this causes a pronounced dependence on resolution in ILES of collapsing gas clouds. Currently, it is completely open whether SGS models for MHD turbulence could amend this problem, but the potential of treating dynamo effects by an SGS model is an exciting prospect. Owing to the anisotropy of MHD turbulence and the potential influence of kinetic effects, however, this is indeed a very challenging problem.