A&A 450, 265-281 (2006)
DOI: 10.1051/0004-6361:20053617

A localised subgrid scale model for fluid dynamical simulations in astrophysics

I. Theory and numerical tests

W. Schmidt^1,2 - J. C. Niemeyer¹ - W. Hillebrandt²

1 - Lehrstuhl für Astronomie, Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Am Hubland, 97074 Würzburg, Germany
2 - Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany

Received 11 June 2005 / Accepted 19 December 2005

Abstract
We present a one-equation subgrid scale model that evolves the turbulence energy corresponding to unresolved velocity fluctuations in large eddy simulations. The model is derived in the context of the Germano consistent decomposition of the hydrodynamical equations. The eddy-viscosity closure for the rate of energy transfer from resolved toward subgrid scales is localised by means of a dynamical procedure for the computation of the closure parameter. Therefore, the subgrid scale model applies to arbitrary flow geometry and evolution. For the treatment of microscopic viscous dissipation a semi-statistical approach is used, and the gradient-diffusion hypothesis is adopted for turbulent transport. A priori tests of the localised eddy-viscosity closure and the gradient-diffusion closure are made by analysing data from direct numerical simulations. As an a posteriori testing case, the large eddy simulation of thermonuclear combustion in forced isotropic turbulence is discussed. We intend the formulation of the subgrid scale model in this paper as a basis for more advanced applications in numerical simulations of complex astrophysical phenomena involving turbulence.

Key words: hydrodynamics - turbulence - methods: numerical

1 Introduction

In the last decade, the significance of turbulence in various astrophysical phenomena from stellar to cosmological scales has been recognised. In retrospect, this is hardly surprising, since virtually all the matter in the Universe is fluid, whereas the solid state is encountered as a rare exception. Moreover, both the integral length L and the characteristic velocity V in astrophysical systems are quite large compared to terrestial standards, while the viscosity $\nu$ is comparable to what is found for liquids or gas on the Earth. For this reason, the dimensionless Reynolds number

$$% {\rm Re}=LV/\nu$$ (1)

becomes very large. A typical figure is ${\rm Re}\sim 10^{14}$ for turbulent stellar convection.

From the computational point of view, the number of degrees of freedom in a fluid dynamical system is given by the relation (see Landau & Lifshitz 1987)

$$% N\sim\left(L/\eta_{{\rm K}}\right)^{3}\sim{\rm Re}^{9/4}.$$ (2)

The length scale $\eta_{{\rm K}}$ is called the Kolmogorov scale and specifies the smallest dynamically relevant length scale. Due to the restriction of the CFL time step for compressible flow, the total number of operations required to compute the evolution over one sound crossing time is of the order

$$% N_{{\rm cr}}\sim N^{4/3}\sim{\rm Re}^{3}.$$ (3)

Hence, $N_{{\rm cr}}\sim 10^{42}$ operations would be required to solve the problem of stellar convection in a direct numerical simulation.

The relevance of the Landau criterion 2 has been questioned recently on the grounds of the intermittency of turbulence (Kritsuk et al. 2006). In fact, the number of degrees of freedom N refers to the ensemble of turbulent flow realisations. At any particular instant of time, however, turbulent dynamics is concentrated in regions of fractal dimension D<3. Topologically, these regions can be either vortex filaments in subsonic flow or shocklets in the case of supersonic turbulence. The fractal dimension of vortex filaments is D=1, whereas D=2 for shocklets. According to the $\beta$ model (see Frisch 1995, Sect. 8.5), the effective number of degrees of freedom is then given by

$$% N_{{\rm eff}}\sim {\rm Re}^{3D/(D+1)}.$$ (4)

If an adaptive numerical scheme with maximal efficiency in tracking the intermittent turbulent regions were applied, the total computational cost could be lowered substantially in comparison to a direct numerical simulation on a static grid. For example, it would appear feasible to treat subsonic turbulence up to a Reynolds number of 10¹⁰ with an anelastic adaptive code on high-end platforms in the near future. On account of the limited efficiency of adaptive schemes, however, the actual constraint might be lower. Apart from that, gravity, magnetohydrodynamic effects and, possibly, reaction networks increase the work load dramatically in astrophysical applications. Moreover, we have $N_{{\rm cr}}\sim N_{{\rm eff}}N^{1/3}\sim{\rm Re}^{11/4}$ for the particularly interesting case of supersonic turbulence. As a consequence, even with sophisticated adaptive schemes, it remains intractable to resolve completely the turbulent fluid dynamics encountered in astrophysics.

In large eddy simulations (LES), on the other hand, only a limited number of degrees of freedom, which correspond to the largest scales of the system, is treated explicitly. For the turbulent dynamics on smaller scales, a so-called subgrid scale model is utilised. Among astrophysicists, the most often used subgrid scale (SGS) model is numerical dissipation. This means that all fluctuations on length scales smaller than the resolution $\Delta$ of the numerical grid are smoothed out and it is assumed that the dynamics on length scales larger than $\Delta$ are more or less independent of the smaller scales. This point of view is motivated by the second similarity hypothesis of Kolmogorov (1941), which holds that the actual mechanism of dissipation is insignificant, provided that there is sufficient scale separation. In other words, on length scales $l\ga\Delta$, it is unimportant whether energy dissipation is caused by the microscopic viscosity $\nu$ at the length scale  $\eta_{{\rm K}}$ or by numerical effects at the cutoff length $\Delta$. The notion of numerical dissipation has been exhaustively investigated for the piece-wise parabolic method (PPM) proposed by Colella & Woodward (1984), which is one of the most popular finite-volume schemes applied in astrophysics. At least for statistically stationary isotropic turbulence, the numerical dissipation of the PPM appears to work well as an implicit SGS model (Schmidt et al. 2006; Sytine et al. 2000).

For the treatment of transient or inhomogeneous turbulent flow, however, an explicit SGS model becomes mandatory. One of the most prominent examples in contemporary theoretical astrophysics are numerical simulations of turbulent combustion in type Ia supernovae (Hillebrandt & Niemeyer 2000). In this case, the velocity scale associated with SGS turbulence determines the effective propagation speed of flame fronts (Niemeyer & Hillebrandt 1995). Part II of this paper will be dedicated to the problem of type Ia supernova simulations. The exchange of energy between resolved and subgrid scales is expected to become dynamically significant in the case of highly intermittent turbulence, for instance, in collapsing turbulent gas clouds in the interstellar medium (Mac Low & Klessen 2004; Larson 2003).

In this paper (Paper I), we present a general framework for the formulation of SGS models based upon the filtering approach of Germano (1992). The mathematical operation of filtering smoothes the flow on length scales smaller than the prescribed numerical resolution $\Delta$. Consequently, a scale separation is introduced, where the smoothed density, velocity, etc. are identified with the resolved quantities computed in LES. We have extended this formalism to compressible flows, using the Reynolds stress model proposed by Canuto (1997) as a guideline. Next we discuss the one-equation SGS turbulence energy model (Schumann 1975; Sagaut 2001). For the energy transfer across the numerical cutoff, we introduce a localised eddy-viscosity closure which makes use of the dynamical procedures introduced by Germano et al. (1991). Hence, the SGS model becomes independent of a priori structural assumptions, in particular, whether the resolved flow is homogeneous or stationary. However, it remains a necessary condition that turbulent regions become locally isotropic on length scales comparable to the numerical cutoff length $\Delta$, because the linear eddy-viscosity closure presumes alignment between the turbulence stress and the rate-of-strain tensors. In the future, multi-parameter closures for the turbulent energy transfer which are not subject to this restriction might be adapted. For the still more complicated non-local transport of kinetic energy on subgrid scales, we use a simple gradient-diffusion closure. In contrast to what has been suggested in the literature (Pope 2000, Sect. 10.3), we find that the optimal diffusivity parameter is larger by about one order of magnitude than the SGS viscosity parameter, i.e. the turbulent kinetic Prandtl number is large compared to unity. Both the localised eddy-viscosity and the statistical gradient-diffusion closure, respectively, were tested by means of data from simulations of forced compressible turbulence. As a case study, we present the LES of turbulent combustion in a periodic box. Although gravitational effects on subgrid scales, in principle, can be incorporated into the model as well, in Paper I we restrict the detailed formulation and application to the case of negligible gravity. However, a simple closure which accounts for unresolved buoyancy effects in simulations of thermonuclear supernovae will be discussed in Paper II.

2 Decomposition of the hydrodynamical equations

Large eddy simulations pose the problem of scale separation. The numerically computed flow can conceptually be defined by a set of smoothed fields which correspond to the low-pass filtered physical flow realisation, where $k_{{\rm c}}=\pi/\Delta$ is the cutoff wavenumber for a numerical grid of resolution $\Delta$. A low-pass filter is a convolution operator which is defined by

$$% q(\vec{x},t) = \left\langle\ensuremath{\stackrel{\infty}{q}... ...right)\ensuremath{\stackrel{\infty}{q}}\left(\vec{x}',t\right)$$ (5)

for a particular kernel $G(\vec{x}-\vec{x}',t)$. The Fourier transform, the so-called transfer function  $\hat{G}(k,t)$, drops to zero for wavenumbers  $k\ga k_{{\rm c}}$. Consequently, only modes of wavenumbers less than  $k_{{\rm c}}$ contribute significantly to the filtered field  $q(\vec{x},t)$. The exact, unfiltered variable  $\ensuremath{\stackrel{\infty}{q}} (\vec{x},t)$corresponds to the limit $k_{{\rm c}}\rightarrow\infty$. The filter operation smoothes out the fluctuations  $q'=\ensuremath{\stackrel{\infty}{q}} -q$ on length scales l smaller than $\Delta$. Albeit being mathematically determined by some dynamical equation, $\ensuremath{\stackrel{\infty}{q}} (\vec{x},t)$ is generally not computable and therefore will be referred to as an ideal quantity.

The dynamical equation for the ideal velocity field is the generalisation of the Navier-Stokes equation for compressible fluids (see Landau & Lifshitz 1987):

$$% \frac{\ensuremath{\stackrel{\infty}{\ensuremath{{\rm D}}}} ... ...\infty}{\vec{v}}} = \ensuremath{\stackrel{\infty}{\vec{f}}} .$$ (6)

The differential operator on the left hand side is the Lagrangian time derivate

$$% \frac{\ensuremath{\stackrel{\infty}{\ensuremath{{\rm D}}}} ... ...t} + \ensuremath{\stackrel{\infty}{\vec{v}}}\cdot\vec{\nabla},$$ (7)

and the effective force density acting upon the fluid is defined by

$$% \ensuremath{\stackrel{\infty}{\vec{F}}} ~ =~ \ensuremath{\s... ... + \ensuremath{\stackrel{\infty}{\vec{F}}} {^{~({\rm ext})}}.$$ (8)

The first term on the right hand side is the pressure gradient, the second term accounts for viscous dissipation and the third term is the total external force per unit volume which encompasses gravitational and, possibly, stirring forces:

$$% \ensuremath{\stackrel{\infty}{\vec{F}}} {^{~({\rm ext})}} =... ...}}} + \ensuremath{\stackrel{\infty}{\vec{F}}} {^{~({\rm s})}}.$$ (9)

The viscous dissipation tensor is defined by

$$% \ensuremath{\stackrel{\infty}{\sigma}} _{ij}~ = 2\ensuremat... ...}{3}\ensuremath{\stackrel{\infty}{{\rm d}}}\delta_{ij}\right),$$ (10)

where $\nu$ is the microscopic viscosity of the fluid,

$$% \ensuremath{\stackrel{\infty}{S}} _{ij}~ = \frac{1}{2}\le... ...ensuremath{\stackrel{\infty}{v}} _{j}}{\partial x_{i}}\right),$$ (11)

are the components of the rate-of-strain tensor and $\ensuremath{\stackrel{\infty}{{\rm d}}} =\partial_{i}\ensuremath{\stackrel{\infty}{v}} _{i}$ is the divergence of the ideal flow.

Equation (6) can be written in a conservative form as

$$% \frac{\partial}{\partial t}\ensuremath{\stackrel{\infty}{\r... ...fty}{\vec{v}}} ~ =~ \ensuremath{\stackrel{\infty}{\vec{F}}} .$$ (12)

The equality of the left hand side in both equations follows from the continuity equation which expresses the conservation of mass

$$% \frac{\partial}{\partial t}\ensuremath{\stackrel{\infty}{\r... ...5mm}{$\rho$ }}}\ensuremath{\stackrel{\infty}{\vec{v}}} ~ = 0.$$ (13)

The goal of the filtering approach is the formulation of dynamical equations for smoothed quantities which are amenable to the numerical computation. We define the Favre or mass-weighted filtered velocity field

$$% \vec{v}(\vec{x},t) \frac{\left\langle\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}... ...gle\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rho$ }}}\right\rangle_{G}}$$ =

$$\frac{1}{\rho(\vec{x},t)} \int\ensuremath{{\rm d}} ^{3}x'~G\left(... ...ec{x}',t\right) \ensuremath{\stackrel{\infty}{\vec{v}}}\left(\vec{x}',t\right),$$ = (14)

where $\rho = \langle\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rho$ }}}\rangle_{G}$ is the filtered mass density. For a homogeneous and time-independent kernel, the filter operation commutes with the Lagrangian derivative. Favre filtering the dynamical Eq. (12), one obtains

$$% \frac{\partial}{\partial t}\rho\vec{v} + \vec{\nabla}\cdo... ...emath{\stackrel{\infty}{\vec{v}}}\right\rangle_{G} = \vec{F},$$ (15)

where $\vec{F}=\langle\ensuremath{\stackrel{\infty}{\vec{F}}}\rangle_{G}$. The filtered equation can be cast into a form analogous to Eq. (6),

$$% \rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}\ve... ...ty}{\vec{v}}} ,\ensuremath{\stackrel{\infty}{\vec{v}}}\right),$$ (16)

by virtue of the generalised turbulence stress tensor

$$% \tau\left(\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$... ...infty}{\vec{v}}}\right\rangle_{G} + \rho\vec{v}\otimes\vec{v}$$ (17)

which Germano (1992) introduced without mass-weighing for incompressible turbulence. We will use $\tau_{ij}=\tau(\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rho$ }}}\ensuremath{\stackrel{\infty}{v}} _{i},\ensuremath{\stackrel{\infty}{v}} _{j})$ as a shorthand notation for the components of the turbulence stress tensor. The term  $\partial_{j}\tau_{ij}$ in the momentum equation stems from the non-linear advection term in the Lagrangian time derivate and can be interpreted as the stress exerted by turbulent velocity fluctuations smoothed out by the filter. In the following, we will encounter a variety of $\tau$-terms. For this reason, we define $\tau(\cdot,\cdot)$ as a generic bilinear form which maps any pair of ideal fields to a mass-weighted smoothed field which is called the generalised second moment. The resulting field can be scalar, vectorial or tensorial. Of course, the notion of a generalised moment applies to moments of higher order as well.

A dynamical equation for the specific kinetic energy, $\ensuremath{\stackrel{\infty}{k}}$ = $\frac{1}{2}\vert\ensuremath{\stackrel{\infty}{\vec{v}}} \vert^{2}$, is readily obtained from the contraction of Eq. (6) with the ideal velocity field  $\ensuremath{\stackrel{\infty}{\vec{v}}}$:

$$% \ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rho$ }}}\... ...infty}{\vec{v}}}\cdot\ensuremath{\stackrel{\infty}{\vec{F}}} .$$ (18)

The mass-weighted filtered kinetic energy  $k(\vec{x},t)$ is defined by

$$% k(\vec{x},t) = \frac{\left\langle\ensuremath{\stackrel{\in... ...el{\infty}{\raisebox{0.05mm}{$\rho$ }}}\right\rangle_{G}}\cdot$$ (19)

Filtering Eq. (18) results in the following equation for  $k(\vec{x},t)$:

$$% \rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}k -... ...c{v}}} ,\ensuremath{\stackrel{\infty}{\vec{v}}}\right)\right].$$ (20)

In addition to the turbulence stress term on the right hand side, there is a non-local transport term which is given by the divergence of the turbulent kinetic flux  $\vec{\mathcal{F}}^{({\rm kin})}$. The flux is defined by the contraction of the completely symmetric generalised third-order moment  $\tau(\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rho$ }}}\ensuremath{\sta... ...suremath{\stackrel{\infty}{\vec{v}}} ,\ensuremath{\stackrel{\infty}{\vec{v}}} )$. In component notation, we have

 

$$% 2\mathcal{F}_{i}^{({\rm kin})} \tau_{ijj} \equiv \tau\left(\ensuremath{\stackrel{\infty}{\raiseb... ...uremath{\stackrel{\infty}{v_{j}}} ,\ensuremath{\stackrel{\infty}{v_{j}}}\right)$$ =

$$-\left\langle\ensuremath{\stackrel{\infty}{\raisebox{0.05mm}{$\rh... ...}\ensuremath{\stackrel{\infty}{v_{j}}}\right\rangle_{G}v_{j} - 2\tau_{ij}v_{j}.$$ = (21)

Since the filtered kinetic energy k is a second-order moment of the ideal velocity field, contributions from velocity fluctuations on all scales are included. For this reason, k differs from the specific kinetic energy of the smoothed flow, $\frac{1}{2}\vert\vec{v}\vert^{2}$, and

$$% k_{{\rm turb}} = k - \frac{1}{2}\vert\vec{v}\vert^{2} = -\frac{1}{2}\tau_{ii}$$ (22)

can be identified with the generalised turbulence energy associated with scales smaller than the characteristic length of the filter  $\langle\ \rangle_{G}$. The dynamical equation for  $k_{{\rm turb}}$ is obtained by subtracting

$$% \rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}\le... ...ec{v}}} ,\ensuremath{\stackrel{\infty}{\vec{v}}}\right)\right]$$ (23)

from Eq. (20). The result, in component notation, reads

$$% \rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}k_{... ...infty}{v}} _{i},\ensuremath{\stackrel{\infty}{F}} _{i}\right),$$ (24)

where $\tau(\ensuremath{\stackrel{\infty}{v}} _{i},\ensuremath{\stackrel{\infty}{F}} _{i})$ is the contraction of the tensor

$$% \tau\left(\ensuremath{\stackrel{\infty}{\vec{v}}} ,\ensurem... ...l{\infty}{\vec{F}}}\right\rangle_{G} + \vec{v}\otimes\vec{F}.$$ (25)

In order to put Eq. (24) into a form which is more adequate for modelling purposes, flux terms are split off on the right hand side. Let us first substitute the definition of the effective force (8):

$$% \tau\left(\ensuremath{\stackrel{\infty}{v}} _{i},\ensuremat... ...suremath{\stackrel{\infty}{\vec{F}}} {^{~({\rm ext})}}\right).$$ (26)

The three resulting generalised moments respectively correspond to pressure, viscous and external forces. Because stirring forces usually supply energy on the integral length L of the flow only, it follows that  $\tau(v_{i},F_{i}^{~({\rm s})})$ is negligibly small for  $\Delta\ll L$. The first and the second term on the right hand side of Eq. (26) can be split as follows:

  

$$-\tau\left(\ensuremath{\stackrel{\infty}{v}} _{i},\partial_{i}\en... ...ft(\ensuremath{\stackrel{\infty}{d}} ,\ensuremath{\stackrel{\infty}{P}}\right),$$ (27)

$$\tau\left(\ensuremath{\stackrel{\infty}{v}} _{i},\partial_{j}\ens... ...stackrel{\infty}{S}} _{ij},\ensuremath{\stackrel{\infty}{\sigma}} _{ij}\right).$$ (28)

The new $\tau$-terms on the right hand side are, respectively, the pressure-dilatation and the rate of viscous dissipation. Substituting the expression for the viscous dissipation tensor (10), it follows that

$$% \tau\left(\ensuremath{\stackrel{\infty}{S}} _{ij},\ensurema... ...h{\stackrel{\infty}{S}} {_{ij}^{\ast}}\right\rangle_{G}S_{ij}.$$ (29)

The norm $\vert\ensuremath{\stackrel{\infty}{S}} {^{\ast}}\vert$ is defined by the total contraction of the trace-free part of the rate of strain tensor:

$$% \left\vert\ensuremath{\stackrel{\infty}{S}} {^{\ast}}\right... ...\vert^{2}-\frac{2}{3}\ensuremath{\stackrel{\infty}{d}} {^{2}}.$$ (30)

The norm of the rate-of-strain tensor is a probe of the velocity fluctuations on the smallest dynamical length scales which are of the order of the Kolmogorov scale  $\eta_{{\rm K}}$. We shall assume that the characteristic length of the filter is much greater than the Kolmogorov scale. In this case, the rate of strain of the filtered velocity field is much less than the rate of strain of the ideal velocity field, i.e. $\vert\ensuremath{\stackrel{\infty}{S}} {^{\ast}}\vert^{2}\la\vert S^{\ast}\vert^{2}$. Consequently, the first term on the right hand side of Eq. (29) dominates the second term, and we can set

$$% \tau\left(\ensuremath{\stackrel{\infty}{S}} _{ij},\ensurema... ...tackrel{\infty}{S}} {^{\ast}}\right\vert^{2}\right\rangle_{G}.$$ (31)

Summarising, Eq. (24) can be written in the form

$$% \rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}k_{... ...}} - \mathfrak{D} = \Gamma + \Sigma - \rho(\lambda+\epsilon),$$ (32)

where the source contributions on the right hand side are

  

$$\Gamma = \left\langle\ensuremath{\stackrel{\infty}{\raisebox{0.05... ... _{i}\ensuremath{\stackrel{\infty}{g}} _{i}\right\rangle_{G} - \rho v_{i}g_{i},$$ (33)

$$\Sigma = \tau_{ij}S_{ij},$$ (34)

$$\rho\lambda = -\left\langle\ensuremath{\stackrel{\infty}{d}}\ensuremath{\stackrel{\infty}{P}}\right\rangle_{G} + d P,$$ (35)

$$\rho\epsilon = \left\langle\ensuremath{\stackrel{\infty}{\nu}}\en... ...ert\ensuremath{\stackrel{\infty}{S}} {^{\ast}}\right\vert^{2}\right\rangle_{G},$$ (36)

and the transport term $\mathfrak{D}$ is given by

$$% \mathfrak{D} = \frac{\partial}{\partial x_{i}} \left[\fra... ...ensuremath{\stackrel{\infty}{v}} _{j}\right\rangle_{G}\right].$$ (37)

The generalised moment

$$% \vec{\mu} = -\left\langle\ensuremath{\stackrel{\infty}{\vec{v}}}\ensuremath{\stackrel{\infty}{P}}\right\rangle_{G}+\vec{v}P.$$ (38)

accounts for the transport of turbulence energy due to pressure fluctuations.

For the internal energy, the filtered dynamical equation is

 

$$\rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}e_{{\rm i... ...l{\infty}{P}}\ensuremath{\stackrel{\infty}{d}}\right\rangle_{G} - \rho\epsilon.$$ (39)

The source term Q accounts for heat generation by chemical or nuclear reactions. The transport of heat due to turbulent temperature fluctuations gives rise to the generalised conductive flux,

$$% \vec{\mathcal{F}}^{({\rm cond})} = -\tau\left(\ensuremath{\... ...tackrel{\infty}{T}}\right\rangle_{G} - \rho\chi\vec{\nabla}T,$$ (40)

for fluid of thermal conductivity $\ensuremath{\stackrel{\infty}{\chi}}$. The additional transport term on the left hand of Eq. (39) side arises from the transport of heat by turbulent advection. In the case of buoyancy-driven turbulence, this transport mechanism is known as convection. The generalised convective flux is defined by

$$% \vec{\mathcal{F}}^{({\rm conv})} = \tau\left(\ensuremath{\... ...uremath{\stackrel{\infty}{e}} _{{\rm int}}\right) + \vec{\mu},$$ (41)

where $h=e_{{\rm int}}+P/\rho$ is the specific enthalpy.

Adding the budgets of the specific kinetic energy $\frac{1}{2}\vert\vec{v}\vert^{2}$ and internal energy  $e_{{\rm int}}$, we obtain the total energy per unit mass on length scales $l\ga\Delta$, i.e. $e_{{\rm tot}}=e_{{\rm int}}+\frac{1}{2}\vert\vec{v}\vert^{2}$. The dynamical equation governing the evolution of  $e_{{\rm tot}}$ is

 

$$\rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}e_{{\rm t... ...ckrel{\infty}{\vec{v}}} ,\ensuremath{\stackrel{\infty}{\vec{v}}}\right)\right].$$ (42)

This conservation law in combination with Eq. (32) extends the Germano consistent decomposition to compressible turbulence. The sum of internal energy and kinetic energy on all scales is $e_{{\rm tot}} + k_{{\rm turb}} = e_{{\rm int}} + k$. Comparing Eqs. (32) and (42), one can see that  $\rho(\lambda + \epsilon)$ accounts for the dissipation of turbulence energy into internal energy by compression effects and viscous dissipation, respectively. The turbulence production term $\Sigma$ is related to the energy transfer through the turbulence cascade across the characteristic length of the filter. The injection of energy due to buoyancy and the action of stirring forces on length scales larger than $\Delta$ is given by $\vec{v}\cdot[\vec{g}+\vec{f}^{({\rm s})}]$, whereas the interaction of gravitational potential energy fluctuations and turbulence on length scales smaller than $\Delta$ is accounted for by the term $\Gamma$.

From the discussion in this section, it should become clear that the presumed scale separation in LES is essentially based upon the disentanglement of a variety of dynamical effects. This task is considerably complicated by the non-linear interactions across the cutoff scale, which become manifest in the various generalised moments occurring in Eqs. (20), (32), (39) and (42). Hence, one faces the problem of finding closures for the generalised moments in the decomposed dynamical equations.

In the simplest of all cases, a closure is a sufficiently convincing argument for neglecting a certain term. This kind of closure is applied in many cases. In the proper sense, a closure is a more or less tentative approximation which is made on grounds of heuristic physical arguments. Two major categories of closures can be distinguished: an algebraic closure is some function of filtered quantities only. Usually, algebraic closures contain at least one free parameter. Depending on whether this parameter is a constant or varies in space and time, the closure is either statistical or localised. On the other hand, dynamical closures determine generalised moments from additional dynamical equations. However, these equations, in turn, entail closures for higher-order generalised moments. This is the problem of the infinite hierarchy of equations for filtered quantities. Inevitably, the hierarchy must be truncated at some point with the help of algebraic closures.

3 The subgrid scale turbulence energy model

In LES, numerical solutions for the filtered quantities $\rho$, P, T, $\vec{v}$ and  $e_{{\rm tot}}$ have to be computed from the continuity equation for $\rho$, the momentum Eq. (16), the energy conservation law (42) and the equation of state. The filter operation naturally introduces a cutoff which is related to the numerical resolution. Since the filtering in LES is not necessarily explicit but sometimes inherent to the numerical scheme, we will subsequently use the generic notation $\langle\ \rangle_{{\rm eff}}$. For example,

$$% \vec{v}(\vec{x},t) = \frac{\left\langle\ensuremath{\stackr... ...infty}{\raisebox{0.05mm}{$\rho$ }}}\right\rangle_{{\rm eff}}},$$ (43)

is the mass-weighted velocity field which is to be computed numerically. For finite precision, the numerical solution actually corresponds to a whole ensemble of exact flow realisations. In this regard, one can think of a reduction of the number of degrees of freedom due to filtering.

The length scales smaller than the characteristic length  $\Delta_{{\rm eff}}$ of the effective filter are the subgrid scales (SGS). The scales $l\ga\Delta_{{\rm eff}}$, on the other hand, are numerically resolved. A complete set of closures for the generalised moments in the dynamical equations constitutes the subgrid scale model. A general SGS model which includes dynamical equations for the moments of second and third order was formulated by Canuto (1994). Unfortunately, the computational cost of solving the whole set of dynamical equations for the generalised moments is considerable. Moreover, the problem of stability appears to be non-trivial.

The SGS model which will be discussed in the following involves the solution of the dynamical equation for the subgrid scale turbulence energy only. The definition of the density of SGS turbulence energy is as follows:

$$% K_{{\rm sgs}} \equiv \frac{1}{2}\rho q_{{\rm sgs}}^{2} = -... ...\right\rangle_{{\rm eff}} - \rho\vert\vec{v}\vert^{2}\right].$$ (44)

The magnitude of SGS velocity fluctuations is given by  $q_{{\rm sgs}}$. The equation governing the evolution of the specific turbulence energy $k_{{\rm sgs}}=K_{{\rm sgs}}/\rho$ is just Eq. (32) with the filter $\langle\ \rangle_{{\rm eff}}$. For the various second-order moments in this equation we will invoke standard algebraic closures. Hence, the implementation of the SGS model requires the solution of only one additional dynamical equation. The inherent limitations of the simple closures are in part compensated by localisation. Thereby, the SGS model becomes basically independent of a priori model parameters which presume certain flow properties such as homogeneity. In essence, the SGS model which will be formulated is robust and particularly well suited for complex flow geometries and transients, although requiring relatively little computational resources.

There are two different sources of SGS turbulence production. The first one is the SGS gravity term $\Gamma_{{\rm sgs}}$ (33) which accounts for the conversion of potential into kinetic energy and vice versa due to correlations between SGS fluctuations of the velocity and gravity. A putative closure for SGS buoyancy in reactive flows will be presented in Paper II. The other production term is the rate of energy transfer $\Sigma_{{\rm sgs}}=\tau_{ij}S_{ij}$across the length scale  $\Delta_{{\rm eff}}$ due to non-linear turbulent interactions. In general, energy transfer is the primary source of SGS turbulence. A common closure is based on the eddy-viscosity hypothesis for the the trace-free part of the SGS turbulence stress tensor (cf. Pope 2000, Sect. 10.1.):

$$% \tau_{ij}^{\ast} \circeq 2\rho\nu_{{\rm sgs}}S_{ij}^{\ast}... ...o\nu_{{\rm sgs}}\left(S_{ij} - \frac{1}{3}d\delta_{ij}\right),$$ (45)

where

$$% \tau_{ij}^{\ast} = \tau_{ij} - \frac{1}{3}\tau_{ii}\delta_{ij} = \tau_{ij} + \frac{2}{3}K_{{\rm sgs}}\delta_{ij}.$$ (46)

This closure is formulated analogously to the viscous stress tensor in a Newtonian fluid. The eddy viscosity  $\nu_{{\rm sgs}}$ is assumed to be proportional to the product of the characteristic length  $\Delta_{{\rm eff}}$ of the numerical scheme and the characteristic velocity of SGS turbulence (cf. Sagaut 2001, Sect. 4.3; and Pope 2000, Sect. 13.6.3), i.e.,

$$% \nu_{{\rm sgs}} \circeq C_{\nu}\Delta_{{\rm eff}}k_{{\rm sgs}}^{1/2} = \ell_{\nu}q_{{\rm sgs}}.$$ (47)

The length scale $\ell_{\nu}=C_{\nu}\Delta_{{\rm eff}}/\!\sqrt{2}$is thus associated with SGS turbulence production.

Among the dissipation terms, the rate of viscous dissipation  $\epsilon_{{\rm sgs}}$ defined in Eq. (36) dominates in subsonic turbulent flows. Assuming a Kolmogorov spectrum, the mean SGS turbulence energy corresponding to a sharp spectral cut-off can be related to the mean rate of dissipation:

$$% \left\langle k_{{\rm sgs}}\right\rangle = \int_{\pi/\Delta... ...ght\rangle^{2/3} \left(\frac{\pi}{\Delta}\right)^{-2/3}\cdot$$ (48)

Hence, assuming that $C\approx 1.65$ (Yeung & Zhou 1997),

$$% \left\langle\epsilon_{{\rm sgs}}\right\rangle = \Sigma\lef... ...ac{\left\langle k_{{\rm sgs}}\right\rangle^{3/2}}{\Delta}\cdot$$ (49)

Conjecturing that the above relation also holds locally (cf. Pope 2000, Sect. 13.6.3), we have

$$% \epsilon_{{\rm sgs}} \circeq C_{\epsilon} \frac{k_{{\rm s... ...lta_{{\rm eff}}} = \frac{q_{{\rm sgs}}^{3}}{\ell_{\epsilon}},$$ (50)

where $\ell_{\epsilon}=2\sqrt{2}\Delta_{{\rm eff}}/C_{\epsilon}$ and  $C_{\epsilon}\sim 1$. Basically, Eq. (50) implies that SGS eddies of kinetic energy $\sim$ $q_{{\rm sgs}}^{2}$ are dissipated on a time scale $\sim$ $\ell_{\epsilon}/q_{{\rm sgs}}$.

Pressure dilatation poses severe difficulties because one needs to find the correlations between pressure fluctuations and compression or rarefaction of the fluid. The first-order hypothesis is that kinetic energy is dissipated if the fluid is contracting (d < 0). In the opposite case (d > 0), internal energy is converted into mechanical energy which produces turbulence. This line of reasoning leads to the closure proposed by Deardorff (1973):

$$% \lambda_{{\rm sgs}} \circeq C_{\lambda}k_{{\rm sgs}}d.$$ (51)

Unfortunately, numerical tests reveal that this closure is extremely crude. Since compressibility effects tend to diminish toward smaller length scales, the above closure will do for the LES of subsonic turbulence. In the case of supersonic turbulence, however, $\lambda_{{\rm sgs}}$ becomes more significant. Alternative closures for pressure-dilatation are described in Canuto (1997).

A customary algebraic closure for the transport term in Eq. (32) is the gradient-diffusion hypothesis (cf. Sagaut 2001, Sect. 4.3)

$$% \mathfrak{D}_{{\rm sgs}} \circeq \frac{\partial}{\partial ... ... sgs}}^{2} \frac{\partial q_{{\rm sgs}}}{\partial x_{k}}\cdot$$ (52)

The characteristic length scale of diffusion is defined by $\ell_{\kappa}=C_{\kappa}\Delta_{{\rm eff}}/\sqrt{2}$ and the SGS diffusivity is given by $\kappa_{{\rm sgs}} = \ell_{\kappa}q_{{\rm sgs}}$. The notion of a turbulent diffusivity of kinetic energy stems from the analogy to the thermal diffusion of internal energy on microscopic scales. This analogy also suggests the definition of a  kinetic Prandtl number,

$$% {\rm Pr}_{~{\rm kin}} = \frac{\nu_{{\rm sgs}}}{\kappa_{{\rm sgs}}} = \frac{C_{\nu}}{C_{\kappa}}\cdot$$ (53)

Summarising, if gravitational effects on subgrid scales are negligible, we obtain the following dynamical equation for the SGS turbulence energy:

 

$$\frac{\ensuremath{{\rm D}} }{\ensuremath{{\rm D}} t}k_{{\rm sgs}}... ..._{{\rm sgs}}d - C_{\epsilon}\frac{k_{{\rm sgs}}^{3/2}}{\Delta_{{\rm eff}}}\cdot$$ (54)

The remaining problem is the determination of the closure parameters  $C_{\kappa}$, $C_{\nu}$, $C_{\lambda}$ and  $C_{\epsilon}$ which are dimensionless similarity parameters, i.e. the values become asymptotically scale invariant in statistically stationary isotropic turbulence.

4 Closure parameters

In this section, methods for the calculation of closure parameters will be discussed. In particular, we will present a so-called dynamical procedure for the computation of the eddy-viscosity parameter $C_{\nu}$. Originally introduced by engineers in order to improve the performance of simple SGS models such as the Smagorinsky model for turbulent flows near walls, the application of dynamical procedures for the localised computation of closure parameters has turned out to be a powerful tool for the treatment of inhomogeneous and non-stationary turbulence. For this reason, we adapted a procedure proposed by Kim et al. (1999) for the LES of turbulent combustor flows. Using data from simulations of forced isotropic turbulence, we found that this procedure yields a significantly better match with the rate of production than the statistical closure with a constant parameter. For the parameter of dissipation, $C_{\epsilon}$, we propose a semi-statistical solution: A time-dependent value is determined from the energy budget of the resolved flow in extended spatial regions. Regarding the non-local transport, the gradient-diffusion closure produces satisfactory results if the parameter  $C_{\kappa}$ is determined appropriately.

   
4.1 Production

The rate of production $\Sigma_{{\rm sgs}}$ corresponds to dissipation of kinetic energy on resolved scales due to the effect of subgrid scale turbulence. Pictorially, unresolved eddies drain energy from larger eddies at the rate  $\Sigma_{{\rm sgs}}$. This idea motivated the eddy-viscosity closure (45) for  $\Sigma_{{\rm sgs}}$. Extending further the analogy between viscous and turbulent dissipation, an experimental assertion known as the law of finite dissipation could be carried over to the production of SGS turbulence: if, in a large eddy simulation of turbulent flow, all the control parameters are kept the same except for  $\Delta_{{\rm eff}}$, which is lowered as much as possible, the energy dissipation per unit mass, $\Sigma_{{\rm sgs}}$, behaves in a way consistent with a finite positive limit. This suggests that the parameter $C_{\nu}$ in the definition of the turbulent viscosity (47) becomes asymptotically scale-invariant in the limit  $\Delta_{{\rm eff}}/L\rightarrow 0$ and assumes a universal value in the stationary limit.

We verified this hypothesis by analysing data from numerical simulations of forced isotropic turbulence. The driving force which supplies energy on the characteristic length scale L is modelled by a stochastic process in Fourier space (Eswaran & Pope 1988; Schmidt 2004). Under the action of this force, the flow evolves on the characteristic time T which is called the large-eddy time scale. In the statistically stationary limit, the flow velocity is of the order V=L/T. In addition, the weight of solenoidal (divergence-free) relative to dilatational (rotation-free) components of the force field can be varied by setting the control parameter $\zeta$ in the range between 1 and 0. Choosing different characteristic Mach numbers V/c₀, where c₀ is the initial sound speed, and values of $\zeta$, we performed several simulations using the piece-wise parabolic method (PPM) with N=432³ grid cells (Colella & Woodward 1984). A realistic equation of state for electron-degenerate matter was used in these simulations (see Reinecke 2001) and the numerical dissipation of PPM provided an implicit subgrid scale model.

It is possible to evaluate generalised moments from the simulation data on a length scale which is large compared to the cutoff scale $\Delta$. To that end, let us introduce new smoothed fields $\rho^{<}$ and  $\vec{v}^{<}$ which are associated with a basis filter  $\langle\ \rangle_{<}$ of characteristic length  $\Delta _{<}$:

$$% \rho^{<} = \left\langle\ensuremath{\stackrel{\infty}{\raise... ...$ }}}\ensuremath{\stackrel{\infty}{\vec{v}}}\right\rangle_{<}.$$ (55)

If $\Delta _{<}$ is sufficiently large compared to  $\Delta_{{\rm eff}}$, then $q\equiv \langle\ensuremath{\stackrel{\infty}{q}}\rangle_{<} \simeq \langle q\rangle_{<}$ for any quantity q. This property of low-pass filters becomes immediately apparent in spectral space in which the filter operation is just a multiplication of Fourier modes with the transfer function. Consequently, the turbulence stress tensor at the level of the basis filter is approximately given by

$$% \tau_{<}\left(\ensuremath{\stackrel{\infty}{\raisebox{0.05m... ...es\vec{v}\rangle_{<} + \rho^{<}\vec{v}^{<}\otimes\vec{v}^{<}.$$ (56)

Evaluating $\tau_{<}(\rho\vec{v},\vec{v})$, it is possible to determine the eddy-viscosity parameter from the rate of energy transfer across the length scale  $\Delta _{<}$:

$$% C_{\nu}^{<} = \frac{\Sigma_{<}} {\rho^{<}\Delta_{<}\sqrt{k_{{\rm turb}}^{<}}~\left\vert S^{<~\ast}\right\vert},$$ (57)

where $\Sigma_{<} = \tau_{<}^{\ast}(\rho v_{i},v_{j})S_{ij}^{<}$ and $\rho^{<}k_{{\rm turb}}^{<}\simeq -\frac{1}{2}\tau_{<}(\rho v_{i},v_{i})$. We computed $C_{\nu}^{<}$from several flow realisations. A sample of average values for the whole domain is listed in Table 1. We found $\langle C_{\nu}^{<}\rangle\approx 0.05$ for developed turbulence, in agreement with the literature (Pope 2000). Only in the case of transonic flow with a mostly compressive driving force the eddy-viscosity parameters appears to be systematically lower. Figure 1a shows a visualisation of the turbulence energy isosurfaces given by $k_{{\rm turb}}^{<}=0.25V^{2}$ with the contour sections of the dimensionless rate of energy transfer $\tilde{\Sigma}_{<}=(T/\rho_{0}V^{2})\Sigma_{<}$ for V/c₀=0.66 at time t=4T. The corresponding contours obtained with the eddy-viscosity closure and $C_{\nu}^{<}=0.0476$ are plotted in panel (b). Clearly, the rate of energy transfer is not well reproduced. Although there is a significant correlation of about 0.8, the magnitude of spatial variations is greatly reduced.

Table 1: Closure parameters for selected flow realisations from three different simulations of forced isotropic turbulence.

	Figure 1: Isosurfaces of the turbulence energy $k_{{\rm turb}}^{<}=0.25V^{2}$ with contour sections of the dimensionless rate of energy transfer. The flow realisation is taken from a simulation with characteristic Mach number  V/c0=0.66 at time t=4T.
Open with DEXTER

In fact, $C_{\nu}^{<}$ exhibits spatiotemporal variations comparable to the mean value. In consequence, the assumption of a constant eddy-viscosity closure parameter is not valid. However, the information about the variation of $C_{\nu}$ is not available in a LES. A solution to this problem can be found by means of a similarity hypothesis which relates the energy transfer across different scales. Let us consider a length scale  $\Delta_{{\rm T}}$ which is somewhat larger than  $\Delta _{<}$. Introducing a suitably defined filter operation  $\langle\ \rangle_{{\rm T}}$ of characteristic length  $\Delta_{{\rm T}}$, the turbulence stress  $\tau_{{\rm T}}(\rho\vec{v},\vec{v})$ is given by an expression analogous to the right hand side of Eq. (56):

$$% \tau_{{\rm T}}\left(\rho^{<}\vec{v}^{<},\vec{v}^{<}\right) ... ...\rho^{({\rm T})}\vec{v}^{({\rm T})}\otimes\vec{v}^{({\rm T})}.$$ (58)

Here $\rho^{({\rm T})}=\langle\rho^{<}\rangle_{{\rm T}}$ and $\vec{v}^{({\rm T})}=\langle\rho^{<}\vec{v}^{<}\rangle_{{\rm T}}/ \langle\rho^{<}\rangle_{{\rm T}}$. The stress tensors associated with the length scales  $\Delta_{{\rm T}}$ and  $\Delta _{<}$, respectively, are related by an identity which Germano (Germano 1992) originally formulated for incompressible turbulence:

$$% \tau_{{\rm T}}(\rho\vec{v},\vec{v})= \langle\tau_{<}(\rho\... ...+ \tau_{{\rm T}}\left(\rho^{<}\vec{v}^{<},\vec{v}^{<}\right).$$ (59)

The first term on the right-hand side is the filtered turbulence stress tensor associated with the length scale  $\Delta _{<}$, whereas the second term accounts for the turbulence stress on intermediate length scales in between  $\Delta _{<}$ and  $\Delta_{{\rm T}}$. For small scaling ratios $\gamma_{{\rm T}}=\Delta_{{\rm T}}/\Delta_{<}$, there is significant correlation not only between  $\tau_{{\rm T}}(\rho\vec{v},\vec{v})$ and $\langle\tau_{<}(\rho\vec{v},\vec{v})\rangle_{{\rm T}}$, but also between $\tau_{{\rm T}}(\rho\vec{v},\vec{v})$ and $\tau_{{\rm T}}(\rho^{<}\vec{v}^{<},\vec{v}^{<})$. In particular, this was demonstrated by Liu et al. (1994) from velocity measurements in round jets.

Based upon these experimental findings, Kim et al. (1999) proposed a similarity hypothesis for the eddy-viscosity parameter:

$$% C_{\nu}^{<} = C_{\nu}^{({\rm T})} = \frac{\Sigma^{({\rm T})... ...T}}k_{{\rm T}}^{1/2}\left\vert S^{\ast~({\rm T})}\right\vert},$$ (60)

where $\vert S^{\ast~({\rm T})}\vert$ is the norm of

$$% S_{ij}^{({\rm T})} = \frac{1}{2}\left[\partial_{j}v_{i}^{({\rm T})}+\partial_{i}v_{j}^{({\rm T})}\right],$$ (61)

and $\Sigma^{({\rm T})} = \tau_{{\rm T}}^{\ast}(\rho^{<}v_{i}^{<},v_{j}^{<})S_{ij}^{({\rm T})}$. The specific turbulence energy  $k_{{\rm T}}$ corresponding to intermediate velocity fluctuations in between the basis and the test filter level is defined by

$$% \rho^{({\rm T})}k_{{\rm T}} = -\frac{1}{2}\tau_{{\rm T}}\le... ...\left\langle\rho^{<}k_{{\rm turb}}^{<}\right\rangle_{{\rm T}}.$$ (62)

The second expression for  $k_{{\rm T}}$ follows from the contraction of the Germano identity (59). Thus, the parameter  $C_{\nu}^{<}$ for the eddy-viscosity closure at the level of the basis filter is determined by probing the flow at the length scale  $\Delta_{{\rm T}}>\Delta_{<}$. This is why  $\langle\ \rangle_{{\rm T}}$ is called a test filter.

	Figure 2: Probability density functions for the rate of energy transfer  $\tilde{\Sigma}_{<}$ across the length scale  $\Delta _{<}$ in dimensionless scaling at two different instants of time.
Open with DEXTER

Using data from the simulations of forced isotropic turbulence, we tested the proposition made above by computing explicitly the rate of energy transfer across a certain length scale  $\Delta _{<}$ and comparing it to the eddy-viscosity closure with the closure parameter calculated at test filter levels for different scaling ratios  $\gamma_{{\rm T}}$. In order to apply approximation (56), we had to choose a basis filter length  $\Delta _{<}$ which was at least an order of magnitude larger than the resolution $\Delta$ in the simulations. On the other hand, a sufficient range of inertial length scales greater than  $\Delta _{<}$is required for the test filter operation. These requirements substantially constrained the choice of  $\Delta _{<}$. Further complications come from the so-called bottleneck effect which causes a distortion of the energy spectrum function for wave numbers close to the cutoff at $k_{{\rm c}}=\pi/\Delta$ (Dobler et al. 2003; Haugen & Brandenburg 2004). A detailed discussion of the kinetic energy spectrum functions and, particularly, the bottleneck effect in turbulence simulations with PPM is given in Schmidt et al. (2006). As one can see in Fig. 2, the match between the probability density functions of the dimensionless rate of energy transfer and the corresponding localised closure is substantially better than for the closure with constant eddy-viscosity parameter. This is highlighted by the statistical moments listed in Table 2. In particular, the variance of the energy transfer is largely underestimated by the statistical closure. This also becomes apparent from the three-dimensional visualisations in Figs. 1b and c, respectively, which suggest that variations of the energy transfer are flattened by a wide margin in the case of a constant eddy-viscosity parameter, while the localised closure reproduces local extrema quite well. On the other hand, it appears that the characteristic length  $\Delta_{{\rm T}}$ of the test filter should not be chosen too large in relation to  $\Delta^{<}$. Otherwise the mean of the energy transfer will be systematically underestimated (Fig. 2 and Table 2).

Table 2: Statistical moments of the dimensionless rate of energy transfer for the probability density functions plotted in Fig. 2b.

The variability of $C_{\nu}^{({\rm T})}$ is illustrated by the probability density functions plotted in Fig. 3. The similarity of the functions suggest a fairly robust behaviour of  $C_{\nu}^{({\rm T})}$ for driven isotropic turbulence. In a fraction of roughly 15 to $20\%$ of the domain, negative values of the closure parameter are found which are commonly interpreted as inverse energy transfer from length scales smaller than  $\Delta_{{\rm T}}$ toward larger scales. This phenomenon, which is also know as backscattering, is predicted by turbulence theory. However, as we shall argue in Sect. 5, backscattering introduces numerical difficulties in combination with PPM. But panel (d) in Fig. 1 demonstrates that the localised closure is superior even when negative values of the eddy-viscosity are suppressed.

	Figure 3: Probability density functions for the localised eddy viscosity parameter calculated from different flow realisations.
Open with DEXTER

For the application in LES, the basis filter corresponds to the effective filter introduced in the previous section, and the test filter is applied to the computed fields  $\rho(\vec{x},t)$ and  $\vec{v}(\vec{x},t)$. Then we have

$$% C_{\nu} = \frac{\tau_{{\rm T}}^{\ast}\left(\rho v_{i},v_{j}... ..._{{\rm T}}^{1/2}\left\vert S^{\ast~({\rm T})}\right\vert}\cdot$$ (63)

The characteristic length scale of SGS turbulence production, $\ell_{\nu}=$, depends on the scaling ratio  $\gamma_{{\rm T}}=\Delta_{{\rm T}}/\Delta_{{\rm eff}}$ and, consequently, $\ell_{\nu}C_{\nu}\Delta_{{\rm eff}}/\!\sqrt{2}$ is proportional to the parameter $\beta=\Delta_{{\rm eff}}/\Delta$. Schmidt et al. (2006) determined $\beta\approx 1.6$ for the statistically stationary turbulent regime in simulations with PPM.

4.2 Dissipation

The localised closure for the rate of production works because the energy transfer across a certain cutoff wavenumber is mostly determined by interactions between Fourier modes within a narrow band around the cutoff. Concerning the rate of dissipation  $\epsilon_{{\rm sgs}}$, we encounter an entirely different problem. In fact, viscous dissipation takes place on length scales which are of the order of the Kolmogorov scale $\eta_{{\rm K}}\la\Delta_{{\rm eff}}$. There is no obvious similarity between the dissipation on resolved scales (due to SGS turbulence) and the dissipation on subgrid scales (due to microscopic viscosity). The simplest of all SGS models, which is known as the Smagorinsky model, assumes a local equilibrium between the dissipation on resolved and subgrid scales, respectively. However, it is the very point of the SGS turbulence energy model that such a balance does not hold locally. Nevertheless, the mean rate of energy transfer can be related to the rate of viscous dissipation in the case of homogeneous turbulence. If the flow is inhomogeneous, equilibrium might be assumed at least for some nearly homogeneous regions. Thus, we attempt to determine the closure parameter  $C_{\epsilon}$ from the averaged energy budget on the test filter level for a suitably chosen flow region.

The method is loosely based on the variational approach of Ghosal et al. (1995). They subtracted the test-filtered SGS turbulence energy Eq. (54) from the corresponding equation for the turbulence energy at the level of the test filter in order to determine  $C_{\epsilon}$ as a function of both space and time. Our approach is an intermediate one, where spatially averages energy equations are considered. For the mean SGS turbulence energy, averaging Eq. (54) yields

$$% \left\langle\rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\... ...(\lambda_{{\rm sgs}}+\epsilon_{{\rm sgs}}\right)\right\rangle.$$ (64)

Here it is assumed that the surface contributions from the transport term cancel out or at least can be neglected. Since

$$% \left\langle\rho\frac{\ensuremath{{\rm D}} }{\ensuremath{{\... ...ensuremath{{\rm d}} t}\left\langle K_{{\rm sgs}}\right\rangle,$$ (65)

we also neglect the effect of advection by the resolved flow. The turbulence energy density associated with the characteristic scale of the test filter is defined by the trace of the Germano identity (59)

$$% -\frac{1}{2}\tau_{{\rm T}}\left(\ensuremath{\stackrel{\inft... ...eft\langle K_{{\rm sgs}}\right\rangle_{{\rm T}} + K_{{\rm T}},$$ (66)

where $K_{{\rm T}}=\rho^{({\rm T})}k_{{\rm T}}$ (62). The spatial average of the turbulence energy (66) is given by the following dynamical equation:

$$% \frac{\partial}{\partial t}\left\langle K_{{\rm sgs}}+K_{{\rm T... ...rho^{({\rm T})}\left(\lambda_{{\rm T}}+\epsilon_{~{\rm T}}\right)\right\rangle.$$ (67)

Equations (64) and (66) in combination with the Germano identity (59) imply

 

$$% \frac{\ensuremath{{\rm d}} }{\ensuremath{{\rm d}} t}\left\langl... ...rho^{({\rm T})}\left(\lambda_{{\rm T}}+\epsilon_{~{\rm T}}\right)\right\rangle.$$ (68)

Substituting the eddy-viscosity closures for the various production terms on the right-hand side, the above equation becomes

$$% \frac{\ensuremath{{\rm d}} }{\ensuremath{{\rm d}} t}\langle K_{{\rm T}}\rangle \simeq \underbrace{\left\langle\rho^{({\rm T})} C_{\nu}\Delta_{{\rm T}}k... ...e - \frac{2}{3}\left\langle K_{{\rm T}} d^{({\rm T})}\right\rangle}_{{\rm (I)}}$$

$$- \left\langle\rho^{({\rm T})}\lambda_{{\rm T}}\right\rangle + \left\langle\rho^{({\rm T})}\epsilon_{~{\rm T}}\right\rangle$$

$$+ \underbrace{\left\langle\left\langle\rho\nu_{{\rm sgs}}S_{ik}^{... ...rho\nu_{{\rm sgs}}\left\vert S^{\ast}\right\vert^{2}\right\rangle}_{{\rm (II)}}$$

$$- \frac{2}{3}\underbrace{\left\langle\left\langle K_{{\rm sgs}}\right\rangle_{{\rm T}}d^{({\rm T})}- K_{{\rm sgs}}d\right\rangle}_{{\rm (III)}}.$$ (69)

Due to the large number of filtered quantities, the complete numerical computation of the source terms in the above equation would be rather demanding. For this reason, we drop the contributions (II) and (III) while only retaining (I), which is presumably the most significant contribution to the rate of energy transfer across  $\Delta_{{\rm T}}$. Then  $\langle K_{{\rm T}}\rangle$ is approximately given by

 

$$% \frac{\ensuremath{{\rm d}} }{\ensuremath{{\rm d}} t}\left\langle K_{{\rm T}}\right\rangle \left\langle\rho^{({\rm T})} C_{\nu}\Delta_{{\rm T}}k_{{\rm T}}^{1/2} \left\vert S^{\ast~({\rm T})}\right\vert^{2}\right\rangle -\frac{2}{3}\left\langle\rho^{({\rm T})}\left(k_{{\rm T}}d^{({\rm... ...t)\right\rangle - \left\langle\rho^{({\rm T})}\epsilon_{~{\rm T}}\right\rangle.$$ (70)

Invoking the closure dimensional closure (50) both for the SGS rate of dissipation  $\epsilon_{{\rm sgs}}$ and the rate of dissipation at the test filter level, we obtain the following expression for the mean rate of dissipation on length scales in between  $\Delta_{{\rm eff}}$ and  $\Delta_{{\rm T}}$:

$$% \left\langle\rho^{({\rm T})}\epsilon_{~{\rm T}}\right\rangl... ...3/2} - \gamma_{{\rm T}}\rho k_{{\rm sgs}}^{3/2}\right\rangle.$$ (71)

Furthermore, setting

$$% \lambda_{{\rm T}}\circeq C_{{\rm\lambda}}k_{{\rm T}}d^{({\rm T})},$$ (72)

the closure parameter $C_{\epsilon}$ is determined by

 

$$% C_{\epsilon} = -\left[\frac{\ensuremath{{\rm d}} }{\ensuremath{... ... T}}\right)^{3/2} - \gamma_{{\rm T}}\rho k_{{\rm sgs}}^{3/2}\right\rangle^{-1}.$$ (73)

Contrary to the eddy-viscosity parameter $C_{\nu}$ which varies both in space and time  $C_{\epsilon}$ is a time-dependent constant for a suitably chosen spatial region. For homogeneous turbulence, there is only one region encompassing the whole domain of the flow. In a stratified medium, it is appropriate to average horizontally. Then  $C_{\epsilon}$ varies with depth. For turbulent combustion problems, such as type Ia supernova explosions, one can distinguish fuel, the burning zone and the burned material within. For each of these three regions a value of the dissipation parameter is calculated as a function of time. For the pressure-dilatation parameter  $C_{\lambda}$, on the other hand, we preliminarily assume the constant, time-independent value $C_{\lambda}=-\frac{1}{5}$ for subsonic turbulence (see Fureby et al. 1997).

	Figure 4: Turbulence energy isosurfaces as in Fig. 1 with contour sections of the dimensionless flux magnitude of turbulent transport.
Open with DEXTER

4.3 Diffusion

As in the case of the energy transfer, we shall first consider the problem of non-local transport at the level of a basis filter of characteristic length  $\Delta _{<}$ which is large compared to the numerical cutoff length. The generalised kinetic flux (21) is given by

$$% \mathcal{F}_{i}^{(\rm kin)~<} = -\frac{1}{2}\left\langle\rho v_... ...\rho v_{i}v_{j}\right\rangle_{<}v_{j}^{<} - \rho^{<}v_{i}^{<}v_{j}^{<}v_{j}^{<}$$ (74)

and the pressure diffusion flux 38 reads

$$% \vec{\mu} = -\langle P\vec{v}\rangle_{<}+P^{<}\vec{v}^{<}.$$ (75)

Assuming that the total flux vector $\vec{\mathcal{F}}^{(\rm kin)~<} + \vec{\mu}^{<}$is aligned with the turbulence energy gradient  $\vec{\nabla}k_{{\rm turb}}^{<}$, the gradient-diffusion closure can be written as follows:

$$% \vec{\mathcal{F}}^{(\rm kin)~<} + \vec{\mu}^{<} \circeq C_... ...a_{<}\sqrt{k_{{\rm turb}}^{<}}~\vec{\nabla}k_{{\rm turb}}^{<}.$$ (76)

Contracting the above relation with $\vec{\nabla}k_{{\rm turb}}^{<}$ and averaging over the domain of the flow, one obtains

$$% C_{\kappa}^{({\rm vec})~<} = \frac{\left\langle\vec{\mathc... ...ec{\nabla}k_{{\rm turb}}^{<}\right\vert^{2}\right\rangle}\cdot$$ (77)

A sample of values for $C_{\kappa}^{({\rm vec})<}$ is listed in Table 1. In agreement with a turbulent kinetic Prandtl number of the order unity, $C_{\kappa}^{({\rm vec})<}$ is of the same order of magnitude as the closure parameter  $C_{\nu}^{<}$(see Pope 2000, Sect. 10.3). Contour sections of the flux magnitude $\vert\vec{\mathcal{F}}^{(\rm kin)~<}+\vec{\mu}^{<}\vert$ and the corresponding closure at the $k_{{\rm turb}}^{<}=0.25$ isosurfaces for V/c₀=0.66 at time t=4T are shown in Fig. 4. However, as one can see from a comparison of the panels (a) and (b), the closure underestimates the diffusive flux by about an order of magnitude. Even more clearly, this is demonstrated by the probability distribution functions plotted in Fig. 5. We also investigated the hypothesis of setting the turbulent diffusivity parameter equal to the localised eddy-viscosity parameter (see Kim et al. 1999; Sagaut 2001, Sect. 4.3). Since negative diffusivity would induce numerical instability, we truncated the diffusivity parameter at zero, i.e. $C_{\kappa}^{<}=C_{\nu}^{({\rm T})~+}$. The resulting visualisation in panel (c) of Fig. 4 and the corresponding graph in Fig. 5, however, show very little if any improvement compared to the statistical closure.

The reason for the discrepancies is the flawed assumption of alignment between the turbulent flux vector and the energy gradient. Setting

$$% C_{\kappa}^{({\rm scl})<} = \frac{\left\langle\left\vert\v... ...urb}}^{<}\right\vert \sqrt{k_{{\rm turb}}^{<}}\right\rangle},$$ (78)

where an equality of the flux magnitude but not the direction is presumed, results in significantly larger turbulent diffusivity (see Table 1). In particular, panel (d) in Fig. 4 and the probability distribution functions shown in Fig. 5 reveal a very close match between the explicitly evaluated turbulent flux and the gradient-diffusion closure with the parameter $C_{\kappa}^{({\rm scl})<}=0.422$. Remarkably, the implied turbulent kinetic Prandtl number is of the order of ten rather than unity.

	Figure 5: Probability distribution functions for $\vert\vec{\mathcal{F}}^{(\rm kin)~<}+\vec{\mu}^{<}\vert$ and the corresponding gradient-diffusion closures with different turbulent diffusivity parameters.
Open with DEXTER

It appears that the gradient-diffusion closure provides a diffusive mechanism which accounts for the intensity of turbulent transport but fails to reproduce anisotropic properties of third-order generalised moment. This is why advanced statistical theories of turbulence abandon the gradient-diffusion closure and introduce dynamical equations for the third-order moments or make use of other, more sophisticated closures (Canuto & Dubovikov 1998; Canuto 1997). Such equations have been suggested for the application in SGS models as well (Canuto 1994). On account of the difficulties solving these equations, however, we prefer the simple algebraic closure (52) with a constant diffusivity parameter

$$% C_{\kappa} \approx 0.4$$ (79)

corresponding to the turbulent diffusivity

$$% \kappa_{{\rm sgs}} = 0.4\rho\Delta_{{\rm eff}}k_{{\rm sgs}}^{1/2}.$$ (80)

According to our numerical investigation, $C_{\kappa}\approx 0.4$ is representative for stationary isotropic turbulence of Mach number $\la$1. In the case of developing turbulence, the effects of turbulent transport are rather marginal, and the deviations introduced by the statistical diffusivity (80) are not overly important for the subgrid scale dynamics. For higher Mach numbers, however, there appears to be a trend towards systematically larger diffusivity.

In a similar fashion as the gradient-diffusion hypothesis, a turbulent conductivity  $\chi_{{\rm sgs}}$ for the generalised conductive flux in fluid of heat capacity  $c_{{\rm P}}$ and thermal conductivity $\chi$ can be introduced:

$$% \vec{\mathcal{F}}^{({\rm cond})}\circeq \rho c_{{\rm P}}\left(\chi+\chi_{{\rm sgs}}\right)\vec{\nabla}T.$$ (81)

For the generalised convective flux  $\vec{\mathcal{F}}^{({\rm cond})}$, a closure might be based upon the super-adiabatic gradient (Canuto 1994). Moreover, in some combustion problems or in simulations of multi-phase media the turbulent mixing of particle species is yet another challenge. These problems are left for future work.

   
5 Turbulent burning in a box

As a simple testing scenario, we performed LES of turbulent thermonuclear deflagration in degenerate carbon and oxygen. In these simulations, we utilised a greatly simplified reaction scheme, where the products of thermonuclear fusion are nickel and alpha particles. The thermonuclear burning zones propagate in a fashion similar to premixed chemical flames. For the chosen mass density, $\rho_{0}\approx 2.9$ $\times$ $10^{8}~{\rm g~cm^{-3}}$, the width of the flames is $\delta_{{\rm F}}\approx 0.006~{\rm cm}$(cf. Timmes & Woosley 1992). Hence, the flame fronts are appropriately represented by discontinuities for the numerical resolution $\Delta=2$ $\times$ $10^{3}~{\rm cm}$ in the simulations we run. The front propagation is numerically implemented by means of the level set method (Reinecke et al. 1999; Osher & Sethian 1988). The domain of the flow is cubic with periodic boundary conditions (BCs). In this scenario, the burning process consumes all nuclear fuel within finite time. We set $X=216\Delta=4.32~{\rm km}$ for the size of the domain, which is comparable to the resolution of the large scale supernova simulations to be discussed in Paper II. Since self-gravity is insignificant on length scales of the order of a few kilometres, we apply an external solenoidal force field in order to produce turbulent flow. Each Fourier mode of the force field is evolved as a distinct stochastic process of the Ornstein-Uhlenbeck type. The characteristic wavelength L of the forcing modes is half the size of the domain. L can be interpreted as integral length scale of the flow. An detailed description of the methodology and a discussion of numerous simulations is given in Schmidt et al. (2005).

Figure 6: Evolution of statistical quantities in a LES of thermonuclear deflagration in a cubic domain subject to periodic boundary conditions with N=2163 numerical cells. In panel  a) the spatially averaged rate of nuclear energy generation in combination with the mass fractions of unprocessed material (carbon and oxygen), alpha particles and nickel are plotted. The mean as well as the standard deviation of the SGS turbulence velocity  $q_{{\rm sgs}}$ are shown in panel  b).

Open with DEXTER

	Figure 7: LES of thermonuclear deflagration in a cubic domain with periodic BCs. Shown are snapshots of the flame fronts with contour sections of the SGS turbulence velocity in logarithmic scaling.
Open with DEXTER

The LES of turbulent combustion is a particularly appropriate case study for the performance of subgrid scale models because the evolution of the system is strongly coupled to the SGS turbulence energy via the turbulent flame speed relation. For the notion of a turbulent flame speed see Pocheau (1994), Niemeyer & Hillebrandt (1995) and Peters (1999). In the framework of the filtering formalism, the underlying hypothesis is the following: if the flow is smoothed on a certain length scale $\Delta$, then the effective propagation speed  $s_{{\rm turb}}(\Delta)$ of a burning front is of the order of the turbulent velocity fluctuations $v^{~\prime}\sim k_{{\rm turb}}^{1/2}$, provided that $\Delta\gg l_{{\rm G}}$. The length scale  $l_{{\rm G}}$ is called the Gibson scale. It specifies the minimal size of turbulent eddies affecting the flame front propagation. In the context of a LES, we have $s_{{\rm turb}}\sim q_{{\rm sgs}}$ for the turbulent flame speed. Consequently, the SGS model determines the propagation speed of turbulent flames. If $l_{{\rm G}}\la \Delta$, on the other hand, the front propagation is determined by the microscopic conductivity of the fuel. The corresponding propagation speed is called the laminar flame speed and is denoted by  $s_{{\rm lam}}$. Since  $s_{{\rm lam}}$ is determined by the balance between thermal conduction and thermonuclear heat generation, conduction effects are implicitly treated by the level set method. For this reason, we do not include the conduction terms in Eq. (42) for the total energy.

Both limiting cases of turbulent and laminar burning, respectively, are accommodated in the flame speed relation proposed by Pocheau (1994):

$$% \frac{s_{{\rm turb}}}{s_{{\rm lam}}} = \left[1 + C_{{\rm ... ...c{q_{{\rm sgs}}}{s_{{\rm lam}}}\right)^{2} \right]^{1/2}\cdot$$ (82)

The coefficient $C_{{\rm t}}$ is of the order unity and determines the ratio of  $s_{{\rm turb}}$ and  $q_{{\rm sgs}}$ in the turbulent burning regime. Peters (1999) proposes $C_{{\rm t}}=4/3$, while $C_{{\rm t}}\approx20/3$ is suggested by Kim et al. (1999). Here we set $C_{{\rm t}}=1$ corresponding to the asymptotic relation $s_{{\rm turb}}\simeq q_{{\rm sgs}}$assumed by Niemeyer & Hillebrandt (1995). For a study of the influence of  $C_{{\rm t}}$ see Schmidt et al. (2005). The laminar flame speed for an initial mass density of 2.90 $\times$ $10^{8}~{\rm g~cm^{-3}}$ is $s_{{\rm lam}}\approx 9.78$ $\times$ $10^{5}~{\rm cm~s^{-1}}$. Choosing a characteristic velocity $V=L/T=100s_{{\rm lam}}$, where T is the autocorrelation time of the stochastic force driving the flow, the Gibson scale becomes $l_{{\rm G}}\sim 10^{-6}L\sim 0.1~{\rm cm}$. Note that the Gibson length is still large compared to the flame thickness. Therefore, the internal structure of the burning zones is not disturbed by turbulent velocity fluctuations, i.e. the flamelet regime of turbulent combustion applies (Peters 1999).

	Figure 8: Figure 7 continued.
Open with DEXTER

	Figure 9: Evolution of the mean dimensionless rate of nuclear energy generation  a) and the ratio of the mass-weighted mean SGS turbulence velocity to laminar burning speed  b) in a sequence of LES with varying resolution.
Open with DEXTER

Running a LES with the parameters outlined above and setting eight small ignition spots on a numerical grid of N=216³ cells, the expectation was that the burning process would initially proceed slowly, but as turbulence was developing due to the action of the driving force, $q_{{\rm sgs}}$ would eventually exceed the laminar flame speed and substantially accelerate the flame propagation. Indeed, this is what can be seen in Fig. 6 which shows plots of statistical quantities as functions of time. The corresponding flame evolution is illustrated in the sequence of three-dimensional visualisations in Figs. 7 and 8, where the colour shading indicates the contour sections of  $q_{{\rm sgs}}$ in logarithmic scaling. Initially, the spherical blobs of burning material are expanding slowly and become gradually elongated and folded by the onsetting flow which is produced by the driving force. As the SGS turbulence velocity  $q_{{\rm sgs}}$ exceeds the laminar burning speed  $s_{{\rm lam}}$ in an increasing volume of space, the spatially averaged rate of nuclear energy release, $\langle P_{{\rm nuc}}\rangle$, is increasing rapidly (Fig. 6). Eventually, $\langle P_{{\rm nuc}}\rangle$ assumes a peak value at dimensionless time $\tilde{t}=t/T\approx 1.8$ which coincides with the maximum of turbulence energy. Subsequently, the flow approaches statistical equilibrium between mechanical production and dissipation of kinetic energy. Thus, the greater part of the fuel is burned within one large-eddy turn-over time of the turbulent flow. This observation in combination with the tight correlation between the growth of the mean rate of nuclear energy release and the SGS turbulence velocity verifies that the burning process is dominated by turbulence.

As a further indicator for the reliability of the SGS model, we varied the resolution in a sequence of LES, while maintaining the physical parameters unaltered. The resulting global statistics is shown in Fig. 9. In particular, the time evolution of  $\langle P_{{\rm nuc}}\rangle$ appears to be quite robust with respect to the numerical resolution. The deviations which can be discerned in the height, width and location of the peak are mostly a consequence of the different flow realisations due to the random nature of the driving force. Actually, even if we had used identical sequences of random numbers to compute the stochastic force field in each simulation, the dependence of the time steps on the numerical resolution nevertheless would have produced different discrete realisations. Thus, we initialised the random number sequences differently and restricted the resolution study to statistical comparisons. The evolution of the mass-weighted SGS turbulence velocity which is plotted in panel (b) of Fig. 9 reveals that turbulence is developing slightly faster in the case N=192³. This can be attributed to a somewhat larger root mean square force field during the first large-eddy turn-over in this simulation. Consequently, the burning process proceeds systematically faster. Note, however, that the level of SGS turbulence becomes monotonically lower with increasing resolution for the almost stationary flow at time $\tilde{t}=3.0$. The deviations for the LES with the lowest resolution (N=120³), on the other hand, are likely to be spurious. For this reason, it would appear that the minimal resolution for sufficient convergence has to be set in between N=120³ and N=160³.

This conclusion is also supported by the turbulence energy spectra plotted in Fig. 10. We computed the normalised energy spectrum functions for the transversal modes of the velocity fields after two integral time scales have elapsed. Details of the computation of discrete spectrum functions are discussed in Schmidt et al. (2006). One can clearly discern maxima in the vicinity of the normalised characteristic wave number $\tilde{k}=Lk/2\pi=1.0$of the driving force. For the LES with N greater than 120³, an inertial subrange emerges in the interval $2\la \tilde{k}\la 6$. The dimensionless cutoff wave number in the case N=216³ is $\tilde{k}=54$. As demonstrated in Schmidt et al. (2006), the numerical dissipation of PPM, which was used to solve the hydrodynamical equations, noticeably smoothes the flow for wavenumbers $\tilde{k} \ga 54/9=6$. This is exactly what is observed in Fig. 10. For N=120³, on the other hand, virtually all wavenumbers not directly affected by stochastic forcing are subject to numerical dissipation, i.e. there is no inertial subrange at all. Considering the more common power-of-two numbers of cells, a grid of N=128³ cells will provide only marginally sufficient resolution, whereas one will be on the safe side with N=256³. In Paper II, however, it is shown that still higher resolutions might be required for LES of non-stationary inhomogeneous turbulence such as in the case of thermonuclear supernova simulations.

	Figure 10: Transversal kinetic energy spectrum functions at time t=2T for the same sequence of LES as in Fig. 9.
Open with DEXTER

It is also argued in Schmidt et al. (2006) that the intrinsic mean rate of dissipation produced by PPM closely agrees with the prediction of the Smagorinsky model for stationary isotropic turbulence. This suggests that the numerical dissipation can be utilised as an implicit SGS model with regardto the velocity field. In fact, the LES presented in Figs. 6-8 was computed without including the SGS stress term in the dynamical Eq. (16), while the total energy  $e_{{\rm tot}}$, which is conserved by PPM, was coupled to the SGS turbulence energy  $k_{{\rm sgs}}$. One can think of  $k_{{\rm sgs}}$ as a buffer between the resolved kinetic energy  $\frac{1}{2}\vert\vec{v}\vert^{2}$ and the internal energy  $e_{{\rm int}}$. Apart from the energy budget, the SGS model influences the resolved dynamics via the turbulent flame speed. For the LES with varying resolution (Figs. 9 and 10), on the other hand, we applied complete coupling of the SGS model, i.e. the turbulent stress term in the momentum equation was included as well. Comparing Figs. 6a and 9a for N=216³, it appears that the burning process is slightly delayed in the latter case. As is discussed at length in Schmidt et al. (2005), the discrepancy can be attributed to a difficulty related to inverse energy transfer. Since backscattering injects energy on the smallest resolved scales, which are sizeably affected by numerical dissipation, the kinetic energy added to the flow is more or less instantaneously converted into internal energy. Thus, the backscattering of energy from subgrid scales to the resolved flow results in an artificially enhanced dissipation which depletes turbulence energy. Using partial coupling, this unwanted effect is simply ignored. For consistency, one must then introduce a cutoff for the eddy-viscosity parameter $C_{\nu}$ in order to dispose of negative viscosities. Mending the shortcoming of the treatment of inverse energy transfer is the subject of ongoing research. For the time being, the partial coupling of the SGS model with backscattering suppressed serves as a pragmatic solution in hydrodynamical simulations with PPM.

6 Conclusion

The localised SGS turbulence energy model offers robustness and flexibility at relatively low computational cost. For this reason, it is particularly suitable for the application in LES of astrophysical fluid dynamics. The energy transfer from resolved toward subgrid scales is modelled with the standard eddy-viscosity closure, where the closure parameter is computed from local properties of the flow. Hence, there are no a priori assumption about the resolved flow incorporated in the model. Non-local transport is treated with the down-gradient closure, using a constant statistical parameter. With a turbulent kinetic Prandtl number significantly larger than unity, it is possible to reproduce the magnitude of diffusive flux quite well. The rate of viscous dissipation appears to be particularly challenging. We found that a semi-statistical approach yields satisfactory results.

The SGS model was implemented in a code for the LES of turbulent thermonuclear combustion in a periodic box using the piece-wise parabolic method (PPM) for the resolved hydrodynamics and the level set method for the flame front propagation. Since PPM produces significant numerical dissipation, we found it favourable to decouple the SGS model form the momentum equation and suppressing inverse energy transfer from unresolved toward resolved scales. In this kind of application, the SGS turbulence energy serves as a buffer between the resolved kinetic energy and the internal energy and supplies a velocity scale for calculation of the turbulent burning speed.

Furthermore, gravitational and thermal effects can be included in the SGS model, although closures specific to a certain physical system have to be formulated. An example is presented in Paper II, where the application of the SGS model to Rayleigh-Taylor-driven thermonuclear combustion in type Ia supernova is discussed. Adapting the model to other applications, possibly with different numerical techniques, is the goal of on-going research.

Acknowledgements

The simulations of forced isotropic turbulence were run on the Hitachi SR-8000 of the Leibniz Computing Centre. For the LES of turbulent combustion we used the IBM p690 of the Computing Centre of the Max-Planck-Society in Garching, Germany. The research of W. Schmidt and J. C. Niemeyer was supported by the Alfried Krupp Prize for Young University Teachers of the Alfried Krupp von Bohlen und Halbach Foundation.

References

• Canuto, V. M. 1994, ApJ, 428, 729 [NASA ADS] [CrossRef] (In the text)
• Canuto, V. M. 1997, ApJ, 482, 827 [NASA ADS] [CrossRef] (In the text)
• Canuto, V. M., & Dubovikov, M. 1998, ApJ, 493, 834 [NASA ADS] [CrossRef]
• Colella, P., & Woodward, P. R. 1984, J. Comp. Phys., 54, 174 [NASA ADS] (In the text)
• Deardorff, J. W. 1973, ASME J. Fluids Engng., 429 (In the text)
• Dobler, W., Haugen, N. E., Yousef, T. A., & Brandenburg, A. 2003, Phys. Rev. E, 68, 026304 [NASA ADS] [CrossRef]
• Eswaran, V., & Pope, S. B. 1988, J. Comp. Physics, 16, 257
• Frisch, U. 1995, Turbulence (Cambridge Universtiy Press) (In the text)
• Fureby, C., Tabor, G., Weller, H. G., & Gosman, A. D. 1997, Phys. Fluids, 9, 3578 [NASA ADS] [CrossRef] (In the text)
• Germano, M. 1992, J. Fluid Mech., 238, 325 [NASA ADS] (In the text)
• Germano, M., Piomelli, U., Moin, P., & Cabot, W. H. 1991, Phys. Fluids, 3, 1760 [NASA ADS] [CrossRef] (In the text)
• Ghosal, S., Lund, T. S., Moin, P., & Akselvoll, K. 1995, J. Fluid Mech., 286, 229 [NASA ADS] (In the text)
• Haugen, N. E., & Brandenburg, A. 2004, Phys. Rev. E, 70, 026405 [NASA ADS] [CrossRef]
• Hillebrandt, W., & Niemeyer, J. C. 2000, ARA&A, 38, 191 [NASA ADS] (In the text)
• Kim, W., Menon, S., & Mongia, H. C. 1999, Combust. Sci. and Tech., 143, 25 (In the text)
• Kolmogorov, A. N. 1941, C. R. Acad. Sci. URSS (In the text)
• Kritsuk, A. G., Norman, M. L., & Padoan, P. 2006, ApJ, 638, L25 [NASA ADS] [CrossRef] (In the text)
• Landau, L. D., & Lifshitz, E. M. 1987, Course of Theoretical Physics, Vol. 6, Fluid Mechanics, 2nd Ed. (Pergamon Press) (In the text)
• Larson, R. B. 2003, Rept. Prog. Phys., 66, 1651 [NASA ADS] [CrossRef]
• Liu, S., Meneveau, C., & Katz, J. 1994, J. Fluid Mech., 275, 83 [NASA ADS] (In the text)
• Mac Low, M., & Klessen, R. S. 2004, Rev. Mod Phys., 76, 125 [NASA ADS] [CrossRef]
• Niemeyer, J. C., & Hillebrandt, W. 1995, ApJ, 452, 769 [NASA ADS] [CrossRef] (In the text)
• Osher, S., & Sethian, J. A. 1988, J. Comp. Phys., 79, 12
• Peters, N. 1999, Journal of Fluid Mechanics, 384, 107 [NASA ADS] [CrossRef] (In the text)
• Pocheau, A. 1994, Phys. Rev. E, 49, 1109 [NASA ADS] [CrossRef] (In the text)
• Pope, S. B. 2000, Turbulent Flows (Cambridge University Press) (In the text)
• Reinecke, M. A. 2001, Ph.D. Thesis, Department of Physics, Technical University of Munich, online available from http://tumb1.biblio.tu-muenchen.de/publ/diss/ph/2001/reinecke.html (In the text)
• Reinecke, M., Hillebrandt, W., Niemeyer, J. C., Klein, R., & Gröbl, A. 1999, A&A, 347, 724 [NASA ADS]
• Sagaut, P. 2001, Large Eddy Simulation for Incompressible Flows (Springer)
• Schmidt, W. 2004, Ph.D. Thesis, Department of Physics, Technical University of Munich, online available from http://tumb1.biblio.tu-muenchen.de/publ/diss/ph/2004/schmidt.html
• Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2005, Combust. Theory Modelling, 9, 693 [CrossRef] (In the text)
• Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comp. Fluids., 35, 353 [CrossRef]
• Schumann, U. 1975, J. Comp. Phys., 18, 376 [NASA ADS]
• Sytine, I. V., Porter, D. H., Woodward, P. R., Hodson, S. W., & Winkler, K. 2000, J. Comp. Phys., 158, 225 [NASA ADS]
• Timmes, F. X, & Woosley, S. E. 1992, ApJ, 396, 649 [NASA ADS] [CrossRef] (In the text)
• Yeung, P. K., & Zhou, Y. 1997, Phys. Rev. E, 56, 1746 [NASA ADS] [CrossRef] (In the text)