Low-Dissipation Simulation Methods and Models for Turbulent Subsonic Flow

The simulation of turbulent flows by means of computational fluid dynamics is highly challenging. The costs of an accurate direct numerical simulation (DNS) are usually too high, and engineers typically resort to cheaper coarse-grained models of the flow, such as large-eddy simulation (LES). To be suitable for the computation of turbulence, methods should not numerically dissipate the turbulent flow structures. Therefore, energy-conserving discretizations are investigated, which do not dissipate energy and are inherently stable because the discrete convective terms cannot spuriously generate kinetic energy. They have been known for incompressible flow, but the development of such methods for compressible flow is more recent. This paper will focus on the latter: LES and DNS for turbulent subsonic flow. A new theoretical framework for the analysis of energy conservation in compressible flow is proposed, in a mathematical notation of square-root variables, inner products, and differential operator symmetries. As a result, the discrete equations exactly conserve not only the primary variables (mass, momentum and energy), but also the convective terms preserve (secondary) discrete kinetic and internal energy. Numerical experiments confirm that simulations are stable without the addition of artificial dissipation. Next, minimum-dissipation eddyviscosity models are reviewed, which try to minimize the dissipation needed for preventing sub-grid scales from polluting the numerical solution. A new model suitable for anisotropic grids is proposed: the anisotropic minimum-dissipation model. This model appropriately switches off for laminar and transitional flow, and is consistent with the exact sub-filter tensor on anisotropic grids. The methods and models are first assessed on several academic test cases: channel flow, homogeneous decaying turbulence and the temporal mixing layer. As a practical application, accurate simulations of the transitional flow over a delta wing have been performed. Mathematics Subject Classification 65M08 · 65M12 · 76F65 · 76F06 · 76G25


Introduction
Reduction of the aerodynamic drag of aircraft is a formidable task, because viscous friction forces are subject to the chaotic process of turbulence, which engineers would like to better understand. Although the Navier-Stokes equations for turbulent flow have been known since 1822, and can be written in a few lines of mathematics, it seems that the origin and evolution of turbulence can only be understood either through detailed physical experiments or through computer simulation of the Navier-Stokes equations. This paper focuses on the latter. 1 3

Convection Versus Diffusion
Flows of air, whether they are laminar or turbulent, are governed by the compressible Navier-Stokes equations. The momentum equation for compressible flow is [6] where is the mass density, ∈ ℝ 3 is the velocity field, p is the pressure, and ∈ ℝ 3×3 is a tensor which models the viscous friction in a fluid. Note that throughout the paper, the symbol ⋅ will denote an inner product in ℝ 3 with corresponding norm | ⋅ | 3 .
The common physical explanation of turbulence is that it is a cascade of progressively smaller and more complex flow structures. The driving force of the turbulent cascade is the non-linear convective term ∇ ⋅ ( ⊗ ) . This term models the transfer of momentum to smaller scales (energy cascade) and conserves both momentum and kinetic energy, i.e. it only redistributes kinetic energy over the scales of motion.
The diffusion term of the Navier-Stokes equations, ∇ ⋅ , models the viscous friction in a fluid. This term conserves momentum, but dissipates kinetic energy. The diffusion of the viscous terms is stronger for smaller scales of motion. For flows at turbulent Reynolds numbers, the viscous diffusion hardly affects the energy of the larger scales. The cascade of kinetic energy to smaller flow structures continues until they become so small that the viscous diffusion becomes noticeable: the so-called Kolmogorov scales.

Turbulence Modeling
The challenging part of flow simulation is the modeling of turbulence. Depending on the amount of flow details that are being resolved, several modeling levels can be distinguished, ranging from RaNS to DNS.

Direct Numerical Simulation
A direct numerical simulation (DNS) of a turbulent flow intends to capture the cascade of kinetic energy from the largest to the smallest Kolmogorov scales. To give a rough indication, the computational complexity of a DNS of homogeneous isotropic turbulence scales with the Reynolds number as Re 11∕4 [11,16,65]. Thus, increasing the Reynolds number by a factor 10 increases the computational complexity of a DNS by a factor of approximately 1000.
In spite of its high computational costs, DNS is already feasible for flows at moderate Reynolds numbers. As an example, below an accurate simulation of the transitional flow over a delta wing at a chord Reynolds number of (1) t ( ) + ∇ ⋅ ( ⊗ ) + ∇p = ∇ ⋅ , Re c = 150,000 will be described. For practical engineering purposes at higher Reynolds numbers, however, a DNS is currently often too expensive. Therefore, engineers usually resort to cheaper coarse-grained models of turbulent flow: for example the Reynolds-averaged Navier-Stokes model, the large-eddy simulation model, or hybrid models.

Reynolds-Averaged Navier-Stokes Models
A cheap way of modeling turbulence is by means of Reynolds averaging the Navier-Stokes equations (RaNS). In a RaNS all turbulent flow scales are being modeled. Successful models for aerodynamic applications are Menter's SST k-model [48], and the Spalart-Allmaras model [82]. In spite of their simplicity, such models can give satisfactory results for attached flow in aircraft cruise conditions. However, RaNS models have insufficient accuracy in simulations in off-design conditions. Especially flows with large separation regions are challenging to capture [83]. Also, RaNS models are not appropriate if time-dependent quantities, such as vibrations or acoustic waves, are of primary interest [37].

Large-Eddy Simulation
A more sophisticated coarse-grained flow model is used in large-eddy simulation (LES). In a LES, the larger scales in a flow are computed explicitly, but the effect of the smaller unresolved scales is modeled [47]. The basic ingredient of LES modeling is a model for the sub-grid dynamics in terms of the resolved velocity field : a socalled eddy-viscosity model. An important example of an eddy-viscosity model for LES is the classical Smagorinsky model [80]. The Smagorinsky model gives good results for homogeneous isotropic turbulence, but removes too much kinetic energy from laminar and transitional flows. This can delay the transition to turbulence of a shear layer. LES models give more accurate results than RaNS models for flows with massive separation. However, this increased accuracy comes at the cost of a higher computational complexity. Therefore, hybrids between LES and RaNS have been developed, of which the most popular one is detached-eddy simulation (DES) in various variants [15,56,81,83]. Attempts to improve the Smagorinsky model have led to the development of a class of low-dissipation models, such as the dynamic Smagorinsky model [22], the wall-adapting local eddy-viscosity WALE model [58], the Vreman model [106], several regularization models [24,29,89,95] and the QR model [96]. The theoretical background and performance of such models is one of the topics of this paper.

Energy-Conserving Discretization
In CFD the coarse-grained flow models and the full Navier-Stokes equations are solved by a computational method. It is difficult to identify a best all-round CFD method, as it may differ significantly from one application to the other. A criterion specific to LES is that energy errors of the numerical method should not overwhelm the energy dissipation of the LES model [7,39,51]. Too much artificial dissipation can delay transition to turbulence of shear layers, and can also cause bad predictions of the near-wall shear stress. In computational aeroacoustic simulations, too much dissipation can inadequately damp the radiated acoustic waves, which causes predictions of lower sound levels [13,37].
These numerical errors cannot be addressed by only increasing the order of accuracy of a method: even in fifthorder accurate upwind methods the numerical dissipation can still overwhelm the dissipation of the LES model [51]. The need for good accuracy and low energy errors resulted in a quest for higher-order accurate central discretization methods without numerical energy dissipation.
The favorable influence of energy-conserving discretizations was already recognized half-way the 20th century. In 1956, the meteorologist Phillips [63] observed that simulations of the vorticity equation became unstable due to spurious generation of energy and enstrophy. He concluded that the numerical instability is related to the spatial discretization of the non-linear convective term, and called it a nonlinear numerical instability. To stabilize the simulations, Phillips proposed to damp the instability by application of a numerical smoothing process [63]. This initiated the study of the influence of discrete convection on stability.
An early example of an energy-conserving, and hence stable, finite-volume discretization of the two-dimensional incompressible Navier-Stokes equations is the staggered discretization for uniform rectangular grids by Harlow and Welch, developed in the 1960s [28]; see also [64]. Closely related is Arakawa's method from 1966 [4], in which he adapted Phillips' spatial discretization of the vorticity equation such that boundedness of the discrete solution could be proved, without the need for numerical smoothing. He also emphasized the conservation of enstrophy next to energy [5]. Generalization of these ideas to non-uniform and curvilinear grids took several decades.
As a starting point for this paper, we mention the higherorder accurate energy-conserving discretization methods developed in the 1990s. Notable examples of fourth-order accurate methods with small energy errors for three-dimensional incompressible flow are the fourth-order accurate method of Morinishi et al. [55], which conserves momentum and kinetic energy to fourth-order accuracy on staggered uniform rectangular grids, and the fourth-order symmetry-preserving method of Verstappen and Veldman [98][99][100] which exactly conserves momentum and kinetic energy on non-uniform staggered rectangular grids. The discrete kinetic energy can only decrease in the latter method, as the numerical solution satisfies a discrete energy bound t || || 2 h ≤ 0 (where || ⋅ || h denotes the discrete L 2 -norm over the flow domain). I.e. the method is stable in a discrete energy norm, and it preserves the mathematical skew-symmetry of the convective terms at the discrete level, as we will see below.

Outline of the Paper
After formulating the equations for compressible flow in Sect. 2, the symmetry principles behind energy-conserving discretization are explained in Sect. 3. Therafter, the energy-conserving discretization methods for incompressible flow are generalized to compressible flow in Sect. 4. A new framework is developed for the analysis of energyconserving methods for compressible flow. This will lead to discretization methods that conserve the primary variables mass, momentum and total energy; but also they convectively conserve the secondary variables kinetic and internal energy. Also, the time integration method can be designed to conserve total energy. As these methods do not numerically dissipate kinetic energy by convection, energy dissipation is exclusively given by molecular diffusion or the sub-grid dissipation of an LES model.
As argued above, also the dissipation of LES models should be confined. In Sect. 6.2 a number of low-dissipation LES models is being reviewed. Also, a new anisotropic minimum-dissipation (AMD) model is proposed. This model generalizes the earlier proposed minimum dissipation QRmodel (Sect. 7.2) to anisotropic grids which are often used in engineering applications.
The proposed methods and models are extensively validated in simulations of academic test cases: DNS results are presented in Sects. 5.1 and 5.2, whereas LES results are presented in Sects. 6.3 and 7.3. To also obtain practical experience with the developed low-dissipation methods, as 'proof-of-the-pudding', DNS simulations of the transitional flow over a triangular delta wing at Re ≈ 150,000 are presented in Sect. 5.3.

Primary Conservation
Our objective is the development of simulation methods for aerodynamic flow, governed by the Navier-Stokes equations for compressible flow [21]: Here is the mass density, the flow velocity, E = 1 2 ⋅ + e the total energy density, e the internal energy density, p the pressure, the tensor with viscous stresses, and the heat diffusion flux. The above equations are called the continuity equation, momentum equation, and total energy equation, respectively.
The second term at the left-hand sides is the convective term, which models the transport of a quantity with the local flow velocity. The third term at the left-hand sides describes the effects of pressure differences. The first term at the righthand sides models how viscous friction resists local flow velocity differences. Finally the last term in the total energy equation models how heat is diffused.
The Navier-Stokes equations are closed by the standard thermodynamical relations for a calorically perfect gas: the pressure is related to the temperature T by the equation of state p = RT , where R is the gas constant; the internal energy is given by e = c v T , where c v is the specific heat at constant volume; the latter is related to the specific heat at constant pressure c p by the specific heat ratio ≡ c p ∕c v ≈ 1.4 for air at room temperature.
The viscous stress is closed by assuming that air is a Newtonian fluid where (T) is the dynamic viscosity which depends on the temperature (according to Sutherland's law), S is the symmetric rate-of-strain tensor, and I is the identity tensor. The heat diffusion is modeled using Fourier's law q = −k∇T , where k is the thermal conductivity. The thermal conductivity is usually set by imposing the Prandtl number Pr ≡ c p ∕k ≈ 0.72.
In Eq. (2), the Navier-Stokes equations for compressible flow have been expressed in conservative form. This form directly expresses conservation of mass, momentum, and total energy in a flow, because all the terms in the equations of motion are either in divergence or gradient form. Therefore they are called primary conservation laws.

Secondary Conservation
Some terms of the compressible Navier-Stokes equations also conserve quantities which are not given by primary conservation laws. Such conservation laws are called secondary conservation properties. In this paper, the secondary conservation properties of the convective terms are of special interest. The convective terms do not just conserve mass, momentum, and total energy, but also kinetic energy and internal energy separately. This is not a property of just one of the equations, but of a combination of both the continuity and the momentum equation. For example, ignoring the viscous terms, the evolution equation for the kinetic energy density follows as where the product rule has been applied multiple times to obtain the third equality (*). The first term on the righthand side of this equation is due to the convective terms, and is a divergence form. Thus, kinetic energy is conserved by convective transport. Because the total energy is also conserved by convective transport, this implies that internal energy is also conserved by convective transport. The second and third terms in the above equations model the effects of pressure differences. The pressure terms conserve kinetic energy in the incompressible limit ∇ ⋅ = 0 , but not in general. Through compression, they are responsible for exchange between kinetic and internal energy.
Conservation of kinetic and internal energy are important secondary conservation properties of convective transport, and they will be preserved by the discretization without artificial dissipation that is developed in the next section.

Energy Stability and Mathematical Symmetries
Studies of the energy errors and numerical stability of simulation methods for compressible flow, including attempts to control the energy error have been published, often inspired by an early reformulation of the flow equations by Feiereisen et al. [19]. Some authors have proposed to discretize the entropy or internal energy equation instead of the total energy equation [30]. However, this does not yet allow for simultaneous conservation of other physical quantities. From 2008 on, methods have been proposed that preserve both primary as well as secondary conservation properties of the convective term [31,34,37,54,67,86]. Below, such methods are summarized.

Incompressible Flow and Symmetries
We will start the quest for energy-conserving discretizations for compressible flow, by first describing similar methods for incompressible flow. A difference is that the latter are usually defined on staggered rectangular computational grids, needed to avoid numerical odd-evendecoupling of the Poisson equation. A disadvantage of rectangular grids is that they do not allow boundary-fitted simulations of flow around practical complex geometries. However, recently also practical energy-conserving methods for collocated curvilinear computational grids have been proposed [37,67], so that this aspect does not pose fundamental problems (Sect. 4.1). Also, energy-conserving variants for unstructured grids, staggered as well as collocated, have been developed [32,91]. The symmetry-preserving finite-volume method for incompressible flow by Verstappen and Veldman [98][99][100], developed in the 1990s, provides energy stability on non-uniform grids. A symmetry-preserving spatial discretization of the incompressible Navier-Stokes equations can be expressed in matrix-vector notation as where is a diagonal matrix with the grid cell volumes, is a grid function of the velocity field, ( ) is the discretization matrix of the non-linear convective term,  of the divergence, − T of the gradient, and  of the viscous terms [100].
If a discretization preserves the skew-symmetry of the convective operator in Eq. (4) at the discrete level, then T ( ) = − T ( ) , so that A symmetry-preserving discretization matrix  of the viscous friction is positive-definite, so that the discrete norm of the numerical solution decreases. Thus a symmetry-preserving discretization is stable in the discrete energy norm. The property of skew-symmetry is closely related to the summation-by-parts (SBP) property introduced by Strand [84], which includes the influence of boundaries. Some (mixed and spectral) finite-element methods also take discrete kinetic energy conservation into account; see for example [61].

Generalization to Compressible Flow
As mentioned above, several attempts to achieve primary and secondary conservation started from a reformulation of the momentum equation. Either from a skew-symmetric form in primitive variables (e.g. [18,54,93]), or by transformation into other variables (such as the entropy variables in [30]). In all cases, the conservation of mass plays an essential role in showing the analytical equivalence of these formulations. Consistency between discrete conservation of mass and discrete conservation of momentum is then a key ingredient for success. The present approach follows a different route, in which a skew-symmetric basic formulation leads to all desired (primary and secondary) conservation properties. As we will see below, this generates not only the desired properties of the spatial discretization, but it will also allow an energy-conserving time integration (Sect. 4.2) as well as an energy-neutral regularization turbulence model (Sect. 6.2).

Square-Root Variables
Similar to the above reasoning for incompressible flow, the energy conservation properties of compressible flow can be encoded in a skew-symmetry of the convective terms if the solution is expressed in square-root variables instead of conservative variables. Thus, instead of using the state vector with the conserved quantities , , and e , the state vector of a compressible flow is a real-valued vector of the square-root variables where it is assumed that all quantities have been non-dimensionalized. The L 2 (V)-norm of the square-root state vector h is equal to the sum of the mass and total energy in a periodic domain V: Because mass and total energy are conserved in a compressible flow, the norm of the vector h is bounded from above on a periodic domain V: Thus, from physical arguments, square-root variables of the compressible Navier-Stokes equations are bounded in an energy norm, just as a solution of the incompressible

3
Navier-Stokes equations. The aim is to preserve this energy stability after discretization.

Skew Symmetry and Conservation of Products
Indeed, the conservation identified above is reflected in skew-symmetry of the convective term in the Navier-Stokes equations when written in square-root variables. Let denote one of the square-root variables, then some elementary analysis shows that the convective part of its evolution equation reads It is easily seen that the adjoint of the convective operator ( ) is equal to −( ) . In other words, the convective operator ( ) in (9) is skew-symmetric. This mathematical property is very important, because it expresses many conservation properties in a single mathematical statement. To see this, consider the evolution of products of two square-root variables and . If the convective transport of and is given by Eq. (9) and non-convective terms are ignored, then the evolution of the product is given by Because the right-hand side in (10) is a divergence and we consider a periodic domain, it immediately follows that the integrated product vanishes: All quantities of interest (mass, momentum, ...) can be written as a product of two of the square-root variables This demonstrates that all interesting quantities are conserved by convective transport because mathematically ( ) is skew-symmetric. Now that all interesting conservation properties have been expressed as a mathematical property, it is possible to develop discretization methods with similar properties. It will be no surprise that the methods to be developed correspond with finite-volume methods with more conservation properties than usual (i.e. supra-conservative). We will describe this relation in Sect. 4.1.

Energy-Conserving Discretization
In this section, energy-conserving spatial discretizations of the compressible Navier-Stokes equations are derived. These energy-conserving discretizations preserve the conservation properties of 1. mass, 2. (linear) momentum, 3. kinetic energy, 4. internal energy, and (therefore) 5. total energy of (the terms of) the Navier-Stokes equations at the discrete level. The choice of preserving these conservation properties is motivated by practical experience with simulations of incompressible flow, and is not guided by deep mathematical ideas. E.g. conservation of vorticity, helicity or enstrophy [60] can also be interesting, but is out of the scope of the present research. Further research will have to reveal which quantity is more essential in this respect. The proposed analysis of energy-conserving methods is theoretical, but one of the objectives is to implement the discretizations in the finite-volume method Enflow of the Netherlands Aerospace Centre NLR [8,37]. Firstly, the discretization of the convective terms is discussed. The use of the square-root variables is analyzed, and new second-order accurate discretizations on collocated grids are proposed. Also, the road from second-order accurate to higher-order accurate energy-conserving discretizations is discussed.

Skew-Symmetric Discretization on a Collocated Curvilinear Grid
In this section a symmetry-preserving discretization of the convective term on collocated grids is proposed. First, a discrete inner product on the space of grid functions is defined as where k is the discrete volume of grid cell k , k and k denote the numerical solution in the grid cell k , is a matrix with grid cell volumes on the diagonal, and and denote grid functions of the numerical solution. This inner product defines a natural discretization of global physical quantities. For example, is a natural discretization of the amount of mass inside the flow domain V.
The discrete counterpart of a continuous skew-symmetric differential operator  is a discrete skew-symmetric matrix: C = −C T . Just as at the continuous level, a skew-symmetric discrete operator can be related to conservation of inner products of square-root variables. Suppose that the convective part of the equations in square-root variables in Eq. (9) is discretized as [compare (4)] where C is the discretization matrix for the convective operator ( ) in Eq. (9), and non-convective terms have been ignored. Compare this notation of the convective operator for compressible flow with the convective operator for incompressible flow in (4). Then the evolution of the discrete integral of the product of and is Thus, if the matrix C is skew-symmetric, then products of square-root variables do not change by convection.

Discretization of the Convective Terms
A spatial discretization on a curvilinear grid can be derived most easily using a transformation from physical space coordinates = (x 1 , x 2 , x 3 ) to computational space coordinates = ( 1 , 2 , 3 ) . The curvilinear grid in physical space is considered to be a continuously differentiable and invertible image ( ) of a uniform grid in computational space , and the order of accuracy is expressed in terms of the mesh spacing in computational space. Conservation laws are revealed upon spatial integration of the conserved quantities in physical space. Integration of physical quantities in computational space requires multiplication of (9) with the determinant  of the Jacobian matrix J ≡ [101,102]. By setting j ≡  j x i , which can be viewed as the area vector in the j-direction in computational space, and using j j = (see "Appendix 1"), with Einstein convention the convective part of the flow equations can be written as This form expresses the analytic skew-symmetry in computational space, cf. (9), and we want this property to hold in a discrete setting too. In first instance, its discrete form is only skew-symmetric if the area vectors j are discretized at cell centers. Because the discretization should be implemented in a finite-volume method, the area vectors should be located at cell faces of a grid cell. Thus some consistent form of interpolation between cell centers and cell faces is required.
After multiplication by 3 , the second term in Eq. (13) can be discretized at the cell center k by a central second-order finite-difference/volume discretization in computational space. In terms of face variables: where F k is the set of faces of the grid cell with index k , whereas nb(f ) is the neighbor of cell k which shares the face f. Further, f is some second-order accurate interpolation of the velocity vector to the faces (see below) and f is the face area vector.
The third term of the skew-symmetric form in Eq. (13) can be discretized by second-order discretization at the cell face f of the cell k normal to the direction j in computational space after which interpolation of these cell-face discretizations to the cell center gives Summation of the two proposed spatial discretizations (14) and (15) gives Observe how the terms with the central k cancel, and how the coefficient of neighbor nb(f ) only depends on face values. Hence, starting from a skew-symmetric analytic operator, we created a skew-symmetric discrete operator. As shown in Sect. 3.2.2, it conserves discrete products of the basic square-root variables (6). This then, ultimately, leads to the desired discrete primary and secondary conservation properties (Sect. 4.1.2).
The proposed discretization is skew-symmetric for any interpolation of the velocity vector f at face f, and any discretization of the grid cell volumes k and area vectors f . It is second-order accurate if the interpolation of the velocity vector f to the cell faces, the discretization of the cell volumes k at cell centers, and the area vector f at cell faces, are secondorder accurate in computational space. Discretization of the cell volumes k and area vectors f is addressed in "Appendix 1". Some possible second-order accurate interpolations of the velocity vector to the cell faces are (14) 1 2

3
and In Sect. 4.1.2 it is shown that the first two interpolation laws correspond to existing finite-volume methods. The third interpolation law is used in simulations with the symmetrypreserving regularization models [72][73][74]76].
For further technical details we refer to Rozema's PhD thesis [71]. Also, details about the implementation of the pressure terms, the viscous terms and the heat diffusion can be found there.

Relation to Existing Finite-Volume Discretizations
By its construction, a skew-symmetric discretization of the convective terms in square-root variables conserves mass, momentum, kinetic energy, internal energy, and total energy. A finite-volume discretization of the convective terms in standard variables conserves only mass, momentum, and total energy. Thus, it will be no surprise that a skew-symmetric discretization of the convective terms leads to a finitevolume discretization of the convective terms. This can be shown by constructing the discrete counterpart of Eq. (10). If the evolution of discrete square-root variables is governed by the proposed spatial discretization in Eq. (16) and exact time-integration is assumed, then products of discrete square-root variables satisfy where the dots denote non-convective terms. Because of its symmetric nature, k nb(f ) + nb(f ) k ∕2 is a local flux function. Therefore, this is a finite-volume discretization of the convective terms for all the possible products of squareroot variables. The corresponding finite-volume discretizations of the mass, momentum, and total energy equations can be expressed as , where point-wise equalities are assumed in cell centers, so that k = √ k √ k . The interpolation of the velocity vector to the cell faces f is a freedom. Indeed, for the interpolation in (18), the finite-volume discretization proposed by Kok [37] is obtained. The interpolation in (17) gives the finitevolume discretization of the mass and momentum equations proposed by Jameson [31], Subbareddy and Candler [86], and Morinishi [54].

Higher-Order Accurate Discretizations
The energy-conserving discretizations proposed in the previous sections are second-order accurate. As for many computational problems, higher-order accurate methods can be more efficient [37,87,103], below we will create fourth-order accurate discretizations of the inviscid terms of the Navier-Stokes equations.
Richardson extrapolation can be used to create such higherorder accurate discretizations, while at the same time it preserves both the primary and secondary conservation properties of the original discretization. It is typically applied in computational space. To streamline the presentation, the proposed second-order accurate discretizations are expressed as where the denotes that the proposed discretizations are related to a single grid cell.
To achieve higher-order accuracy, Richardson extrapolation combines discretizations with different stencils: where R 2 k and R 3 k are obtained by applying the discretization R k on stencils two and three times as big as the original stencil (see Fig. 1).
By nulling the leading terms in the Taylor expansion of (22), a one-parameter family of discretizations is obtained [37]: Here is a free parameter, which can be used to optimize accuracy in wave number space. For the discretization of a first derivative, the error in wave number space is small if the free parameter is set to = − 0.6668 [37]. For this value, in computational space the derivative is equivalent to the dispersion-relation-preserving discretization proposed by Tam and Webb [88].

Conservative Time-Integration Methods
A symmetry-preserving spatial discretization eliminates conservation errors from the discretization of spatial derivatives, assuming exact time integration. However, in practice a discrete time-integration method is used, thereby possibly introducing discrete conservation errors. We will show that the square-root variables allow time-integration methods that discretely preserve conservation laws upon time stepping. Time-integration methods that preserve energy are called symplectic methods. Recently, Runge-Kutta variants have been applied as conservative time-integration methods for the incompressible Navier-Stokes equations by Sanderse [77]. An important example is the second-order accurate mid-point rule. For variable-density incompressible flow, an early observation of conservation of discrete energy using mid-point integration combined with the square root √ was made by Guermond and Quartapelle [26]. They, however, sacrificed momentum conservation, while the present approach conserves both primary and secondary quantities.
If symplectic mid-point integration is applied in squareroot variables, then the continuity equation in the grid cell k is discretized as where the bar n+ 1 2 = ( n+1 + n )∕2 denotes averaging with equal weights to the time t n+ 1 2 . This can be expressed in conservative variables through multiplication by 2 √ n+ 1 2 k : where the skew-symmetric spatial discretization in Eq. (16) was used. This is a mid-point time integration of a general finite-volume discretization of the continuity equation if the velocity interpolation [compare (18)] where the tilde denotes square-root-weighted temporal averaging of the velocity field A similar averaging was used in the conservative secondorder time-integration methods proposed by Morinishi [54] and Subbareddy and Candler [86]. Thus, it seems that conservation of mass, momentum and kinetic energy upon temporal time integration requires computations with squareroot variables. The above mid-point integration generates a secondorder accurate conservative time-integration method. Also higher-order accurate symplectic time-integration methods exist [77], and so the proposed square-root variables allow for straightforward derivation of time-integration methods of arbitrary order. This has also been observed, independently, by Reiss et al. [9,68].
It must be stressed, however, that symplectic Runge-Kutta methods are implicit, and therefore conservative time integration can be considerably more expensive than using an explicit 1 3 time-integration method; their efficiency depends on the application. In this paper, implicit time integration is not used. Instead, a small time step is used to keep the conservation errors due to explicit time integration sufficiently small [38].

DNS Results
The symmetry-preserving discretization for incompressible flow is known to be a stable and accurate method for simulations of channel flow [100]. To test if these properties generalize to compressible flow, simulations of (subsonic) compressible channel flow are performed. Also, simulations of decaying grid turbulence at a high Reynolds number are performed. These test cases also assess the sensitivity of the results to the grid resolution, in particular with respect to under-resolved grids. Here, we present an overview of some results; a more detailed description can be found in [73].

Turbulent Channel Flow at Re ≈ 180
To examine the applicability of the symmetry-preserving discretization to wall-bounded flow, simulations of a turbulent channel flows at relatively low Reynolds numbers have been performed. The channel flow is nearly incompressible with a bulk Mach number of M b = 0.2 and a friction Reynolds number Re ≈ 180 . As in the DNS by Moser et al. [57], the bulk Reynolds number is fixed at Re b = 2800 . For completeness, the bulk Reynolds number and the bulk Mach number are defined by where the angled brackets denote spatial averaging over the channel.

Set Up of Simulations
The channel is rectangular with a half-height H. The dimensions of the channel are denoted L x = 4 , The channel is periodic in the stream-wise and span-wise directions, and the boundaries in the wall-normal directions are isothermal walls at temperature T w . The initial condition is a Poiseuille flow with uniform density 0 and temperature T w : where w is the fixed molecular viscosity at the wall, and Re b is the prescribed bulk Reynolds number.
The channel flow is driven by a uniform body force in the stream-wise direction with a magnitude that fixes the bulk Reynolds number at the prescribed value. After transition to turbulence, the channel flow is expected to become statistically stationary. Important outputs of a channel flow simulation are the friction velocity and the friction Reynolds number in the turbulent regime, defined by where now the angled brackets denote spatio-temporal averaging over the walls of the channel. These quantities reflect the friction at the walls of the channel as a result of imposing a bulk mean flow.

Numerical Parameters
The simulations in this section have been performed with the fourth-order accurate dispersion-relation-preserving discretization proposed in Sect. 4.1.3 [36,37]. Time integration is performed with a low-storage Runge-Kutta method, and the time step size is set so that the Courant number is below unity. All the simulations have been performed without artificial dissipation.
The grid spacing in the stream-wise and span-wise directions is uniform. The locations of the vertices in the lower half of the channel are given by a hyperbolic sine distribution and reflected with respect to the center plane [100]. Simulations have been performed on two fine grids A and B, a medium grid C, and a coarse grid D (Table 1).
At this low bulk Mach number, the channel flow is practically incompressible, and therefore the incompressible DNS by Moser et al. [57] can be used for validation. The flow has transitioned to turbulence within ⟨u⟩ xyz,0 t∕H = 800 dimensionless time units, and flow statistics are recorded from 800 to 1600 dimensionless time units.

Simulation Results
The simulations are stable without artificial dissipation on all the grids. The obtained friction Reynolds numbers are listed in Table 2. The predicted friction Reynolds numbers become more accurate as the grid is refined, and practically agree with the DNS on the fine grids A and B. where primes denote residuals with respect to the temporal averaging. The = ⟨ ⟩ w ∕u is the viscous length scale, where the angled brackets denote averaging over the walls of the channel. The results obtained on the fine grids A and B accurately agree with the DNS.
The coarser grids C and D do not completely resolve the turbulent energy cascade in this channel flow. In general such simulations can give unstable or erroneous results. However, the simulations with the symmetry-preserving discretization are stable without artificial dissipation, and even though the results obtained on the coarser grids C and D are not perfectly accurate, they are definitely acceptable. The error in the friction Reynolds number is below 3% , and the slope of the mean velocity profile computed on the grids C and D in the log layer deviates only slightly from the results of the DNS: in particular, the log layer is somewhat shorter.
For channel flow, no-model simulations with the symmetry-preserving discretization can give more accurate results than simulations with an LES model on coarse grids [49]. In the current simulations, the relatively accurate predictions of mean flow velocities on a coarse grid D seem to go hand in hand with an over-prediction of the stream-wise velocity fluctuations near the wall of the channel. A similar trade-off of accuracy in mean flow and overshoot of the near-wall flow fluctuations has also been observed in under-resolved simulations with the incompressible symmetry-preserving discretization [100].

Decaying Grid Turbulence
To assess the applicability of the symmetry-preserving simulation method to under-resolved flow, simulations of the decaying grid turbulence experiment by Comte-Bellot and Corrsin [14] have been performed. They generate turbulence by placing a grid with a mesh spacing of M = 5.08 cm in a flow of mean velocity U 0 = 1000 cm s −1 [14]. As the grid turbulence is convected with the mean flow, its intensity gradually decreases. Energy spectra are measured at three stations 42M, 98M, and 171M down-stream of the grid.

Set Up of Simulations
The simulations can be simplified by considering the flow inside a small box that moves along with the mean flow. This box of turbulence is assumed to pass the grid at t = 0 s . Thus, the turbulence in the box is expected to match the measured energy spectra at t = 42M∕U 0 , Table 2 The computed friction Reynolds numbers in simulations with the symmetry-preserving discretization on a grid ranging from a fine grid (A) to a relatively coarse grid (D). Also the relative difference with the friction Reynolds number from a DNS [57]    The initial condition of the simulations is set equal to the flow at the first measurement station. The initial condition is generated by fitting a real-valued and divergence-free velocity field with randomized phases to the energy spectrum measured at the first station [40]. The random phases are adjusted by performing a preliminary simulation from t = 0 to t = 42M∕U 0 with the generated initial condition. Then the amplitudes of the Fourier modes of the solution at t = 42M∕U 0 are rescaled to the desired spectrum as in [33]. The resulting rescaled field is used as the initial condition at time t = 42M∕U 0 .

Simulation Results
The simulations have been performed with the fourth-order symmetry-preserving discretization on a uniform 64 3 computational grid. This grid is too coarse to capture the full energy cascade. Nonetheless, the simulations are stable without artificial dissipation. Figure 3 shows the results as plots of the energy decay and the energy spectra at the times corresponding to the experimental measurements. The reference energy levels in the energy decay plot are obtained by fitting the unfiltered experimental spectra to the present computational grid and computing the resolved energy.
In contrast with the accurate results obtained in underresolved simulations of channel flow, the under-resolved simulation of decaying grid turbulence with the symmetrypreserving discretization disagrees with the experimental measurements. The initial energy decay is considerably smaller than in the experiment, which leads to over-prediction of the total kinetic energy at the second and third measurements station. The computed energy spectra show a considerable accumulation of kinetic energy near the grid cut-off.
The kinetic energy of the turbulent structures near the scale of a grid cell is on average transferred to smaller subgrid scales. However, the symmetry-preserving discretization conserves this energy, so that it is trapped in the simulation, and has to be distributed over the resolved scales. This gives a smaller energy decay rate compared to the experimental data, and the results are unsatisfactory due to pile-up of kinetic energy at the scale of a grid cell. In the sequel (Sect. 6), it is shown that the accuracy of a simulation of decaying grid turbulence on a coarse grid can be improved considerably by using an LES model.

DNS of Flow Over a Delta Wing
The above examples are often-studied academic model problems to validate turbulent flow simulation. In the sections to follow, we will extend the above discussion of low-dissipation discretization methods with low-dissipation turbulence models. The same test cases will there be studied in conjunction with various modern models (Sects. 6 and 7).
But first, to illustrate the current state-of-the-art of DNS, we demonstrate the applicability of the developed discretization methods to large-scale (direct numerical) simulation  The computed kinetic energy decay (a) and energy spectra at times 42M∕U 0 , 98M∕U 0 , and 171M∕U 0 (b) in a highly underresolved simulation of decaying grid turbulence with the symmetrypreserving discretization. The vertical dots correspond to the onedimensional point-to-point oscillation (the grid cut-off) 1 3 of practical turbulent flows. Therefore simulations of the subsonic transitional flow over a simple delta wing at Re = 150, 000 have been performed. In particular, the slender sharp-edged delta wing used in the experiments by Riley and Lowson [69] is studied. The simulations have a high computational complexity, and have been performed on the Dutch national supercomputer at SARA in Amsterdam.

The Delta Wing and Its Aerodynamics
The focus of this simulation is on the transition to turbulence of the flow above a delta wing. Therefore, the chord Reynolds number, Mach number, and angle of attack have been selected to exclude other aerodynamic phenomena.
The wing has a root chord length of c = 471 mm and a thickness of t = 11.5mm. The sweep angle is = 85 • , which makes the delta wing very slender. The bevel angle is 30 • . The simulations have been performed at a relatively small angle of attack = 12.5 • . The Reynolds number Re c = ∞ u ∞ c∕ ∞ based on the chord length c of the delta wing is set to 150,000; in [71] also simulations at Re c = 211, 200 are presented. The free-stream Mach number is M = 0.3 . At these parameter values, the vortical flow structures do not break down above the wing, as the angle of attack is relatively small [44,50]. Also, shock waves are absent [3].
The flow over the delta wing is dominated by a system of conical vortices above the delta wing [3,69]. This vortex system is formed by the shear layer that separates at the leading edge of the delta wing. Under the influence of pressure differences, this shear layer rolls up into two large conical counter-rotating flow structures (see Fig. 4), known as primary vortices. An important aerodynamic property of the primary vortex is that the flow velocity in the core of the primary vortex can be significantly higher than the free-stream velocity. This reduces the pressure in the primary vortex and along the upper surface of the wing, thereby increasing lift.
An interesting effect of the primary vortex is its break-up into discrete sub-vortices and transition to full turbulence. Figure 4 shows axial slices of the instantaneous vorticity magnitude obtained in simulations on a fine grid (see below). Figure 5 shows iso-surfaces of the Q-criterion (the second invariant of the velocity gradient). The figures depict the time-averaged vortical tubes as pairs of sub-structures: a thin sub-structure with a high vorticity and a sub-structure with low vorticity. These vortical tubes co-rotate with the flow, just as the steady sub-vortices observed in the LDV measurements by Riley and Lowson [69]. In the simulations by Visbal and Gordnier [104], time-averaged sub-vortices that co-rotate with the flow have also been observed.
The separated shear layer becomes unsteady at approximately x ≈ 0.6c (≈ 280mm) and has broken into discrete sub-vortices at x ≈ 0.8c (≈ 370mm). For comparison, in the experiments of Riley and Lowson [69], full turbulence is observed downstream of approximately x ≈ 0.86c (≈ 400 mm). The sub-vortices deform as they further travel along the shear layer, and at x = c the deformations have caused so many irregularities in the flow that identification of the discrete sub-vortices becomes challenging.

Numerical Method
A cubic computational domain with a length of 21 chord lengths has been used. At the boundaries of the computational domain, far-field boundary conditions based on Riemann invariants are applied. The initial condition of the simulations is obtained by performing a RaNS simulation with a k-model. No perturbations are added to the initial condition.
The origin of the coordinate system is at the apex of the delta wing: the x-axis is aligned with the chord line, the z-axis is normal to the upper surface, and the y-axis is The grids used are approximately isotropic throughout the primary vortex downstream from the leading edge region (after x ≈ 65mm), as this is the most interesting region of the flow. Also the boundary layers are captured with sufficient resolution. To give an impression, the dimensions of the grid cells of the fine grid in the primary vortex increase approximately linearly from x = 0.08mm, y = 0.05mm, and z = 0.10 mm at the end of the leading edge wedge block to x = 0.57mm, y = 0.40mm, and z = 0.44 mm at the trailing edge.
The simulations have been performed with the fourthorder accurate dispersion-relation-preserving finite-volume discretization from Sect. 4.1.3 [37]. As the grid is still a bit coarse for a genuine DNS, the use of some background artificial or modeled dissipation is appropriate [38], although from the stability point-of-view this is not necessary. Sixthorder artificial dissipation with k 6 = 1∕8 preserves the fourth-order accuracy of the scheme, and localizes the dissipation at the small turbulent structures [38].
Explicit time stepping is used for the simulations. The time step size normalized by the free-stream flow velocity and the chord length of the wing is t u ∞ ∕c = 8.0 × 10 −6 on the coarse grid, t u ∞ ∕c = 4.5 × 10 −6 on the medium grid, and t u ∞ ∕c = 2.2 × 10 −6 on the fine grid. The simulations on the coarse and medium grid are performed from time tu ∞ ∕c = 0 to tu ∞ ∕c = 20 , which corresponds to 20 convective time units. The simulation on the fine grid is performed for 17 convective time units. After 2 convective time units the flow has transitioned to turbulence, and the collection of flow statistics starts. More details on these simulations can be found in Rozema's PhD thesis [71].

Grid Convergence
The quality of the numerical solution is assessed by studying the convergence upon grid refinement of velocity profiles and flow statistics. The results of the simulations are used to study the development of the separated shear layer. Comparison with the experimental measurements by Riley and Lowson [69] can be found in [71].
Velocity Profiles For an accurate simulation, grid convergence of the time-averaged velocity field is expected. Figure 6 shows the time-averaged axial velocity on a vertical line through the core of the primary vortex for the three grids. Although the agreement is not perfect, overall the time-averaged velocity obtained on the coarse grid agrees with the time average obtained on the medium and fine grids.
Turbulent Flow Statistics A more challenging test is to monitor the behavior of the turbulent flow statistics upon grid refinement. They have been recorded for 18 convective time units on the coarse and medium grid, and for 15 convective time units on the fine grid. These time intervals were found long enough for the statistics to settle [71]. Figure 7 shows the turbulent kinetic energy on lines through the separated shear layer, just above the leading edge of the wing at z = 7.6 mm , obtained on different grids. Although the turbulent kinetic energy is more sensitive to the grid resolution than time averages of the velocity and pressure, the turbulent kinetic energy converges as the grid is refined: the turbulent kinetic energy obtained on the medium and fine grid agree accurately. The shear layer becomes unsteady at approximately the same axial location for the medium and the fine grid, around x ≈ 250mm, but the transition to unsteady flow seems to be delayed on the coarse grid. This indicates that, not surprisingly, the coarse

Eddy-Viscosity Models
When a direct numerical simulation (DNS) is not feasible, a (much) cheaper alternative is a large-eddy simulation (LES), e.g. [23,47]. In LES, the computational grid resolves only the large eddies in a flow, and the effect of the smaller scales is modeled. In the common mathematical explanation of LES, this is formalized by application of a spatial filter to a solution of the incompressible Navier-Stokes equations. In practice the filter is often related to the computational grid by setting the filter width equal to the local mesh spacing, so that the filtered solution ̄ can be captured on the computational grid with sufficient accuracy. The residual � = −̄ represents sub-grid scales which cannot be accurately captured on the computational grid.
The evolution equation for the filtered velocity field can be derived by filtering of the Navier-Stokes equations. For the incompressible Navier-Stokes equations and a filter that commutes with spatial derivation, this gives (after normalizing the density) the following equation that is used as the basis for LES: where the molecular diffusion is expressed in terms of the resolved rate-of-strain tensor S = (∇ + (∇ ) T )∕2 . To emphasize that the LES solution is different from the filtered solution of the Navier-Stokes equations, it is denoted by ≈̄ .
The sub-grid tensor ( ) approximates ⊗ −̄ ⊗̄ and is present because the non-linearity in the convective term does not commute with the spatial filter. The challenge of LES is to find a suitable model for the sub-grid tensor in terms of the resolved LES solution .
Currently, there is no consensus on a best LES model, or even on what a proper model should do [66]. Moreover, the results of a practical LES are not completely determined by the sub-grid model, but also by for example the used numerical method, the implementation of the sub-grid model, the computational grid, and the applied boundary conditions. This makes a meaningful comparison of LES models based on results from the literature difficult. Nonetheless, in this section results from some sub-grid models are given, with emphasis on models with low eddy dissipation.
The textbook explanation of the turbulent energy cascade suggests that the larger scales in turbulence (on average) transfer kinetic energy to the sub-grid scales, and therefore the effect of the sub-grid scales on the large eddies is essentially dissipative. Eddy-viscosity models mimic the dissipative nature of the sub-grid scales by adding an eddy viscosity e to the molecular viscosity . This is equivalent to selecting an LES model of the form ( ) = −2 e S , where S is the resolved rate-of-strain tensor.
The Smagorinsky Model The classical example of an eddy-viscosity model is the Smagorinsky model [80], which sets the eddy viscosity equal to where is the filter length and C s is the Smagorinsky constant. The Smagorinsky constant is approximately scale invariant in the inertial range of decaying homogeneous isotropic turbulence [46]. In the early days of LES, the Smagorinsky model was used successfully with a constant C s = 0.20-0.22 to perform simulations of decaying homogeneous isotropic turbulence [42]. However, this Smagorinsky constant is too large to obtain decent results in other flows, such as a channel flow [17,52]. Also, for the temporal mixing layer the Smagorinsky constant C s = 0.20 is too high, The turbulent kinetic energy along a line above the leading edge at z = 7.6 mm (a) and along a line through the separated shear layer above the vortex core (b) for various grid resolutions and delays the transition to turbulence [105]. Therefore, the constant was later lowered to C s ≈ 0.17 [43]. These experiences motivate the developments of LES models that give eddy dissipation more appropriately. Because such models generally aim at delivering smaller eddy dissipation, we call these models low-dissipation models. We will shortly describe several of them in the following subsection. Thereafter, we will introduce a new model that minimizes the amount of eddy viscosity on anisotropic grids.

Low-Dissipation Turbulence Models
Below, we will review a number of low-dissipation LES models, i.e. models that are economical in adding eddy dissipation. We start with some early models: the dynamic Smagorinsky model [22,43], Nicoud's WALE [58] and -models [59] and the Vreman model [106]. Thereafter, regularization models are presented [20,24,90], which do not add any eddy-viscosity at all, and as a consequence are found to not always produce sufficient dissipation. In Sect. 6.3 some results with these models are presented, thereby motivating the quest for models that try to restrict eddy viscosity to a bare minimum: QR [96] and AMD [75].
The Dynamic Smagorinsky Model The issues with the Smagorinsky model have been addressed in different ways. An important improvement is the dynamic procedure [22,43]. It assumes scale invariance of the Smagorinsky constant C s in the inertial range and introduces a coarser filter level by coarse-graining the LES solution with a test filter. By scale invariance the Smagorinsky constant should be equal at the two filter levels, and this can be used to estimate the value of C s (the resulting expression is rather complicated and will be omitted here).
There are some perceived disadvantages of the dynamic Smagorinsky model. The Smagorinsky model requires explicit application of the test filter and evaluation of the model on the filtered field, which increases its computational complexity compared to the static Smagorinsky model. Also, often spatial averaging and clipping is required for stability [45]. Nonetheless, the dynamic Smagorinsky model can give accurate results in simulations of homogeneous isotropic turbulence [53], turbulent channel flow [22], and a turbulent mixing layer [105].
The WALE Model An alternative way to improve upon the static Smagorinsky model is to replace the eddy viscosity in Eq. (26) with a function which appropriately adapts to the LES solution, for instance by switching off (i.e. giving no eddy dissipation) in transitional and laminar flow. An example of this is the wall-adapted WALE model, which by construction vanishes at a desired rate near solid walls, so that the use of an additional wall damping function is not required [58].

The Vreman Model
The idea of an eddy-viscosity model which switches off in laminar flow has been formalized mathematically by Vreman [106]. Vreman derives the flows for which the exact eddy dissipation vanishes, and constructs a model that switches off for the same flows.
The Vreman model gives positive eddy dissipation in flows where the exact eddy dissipation is negative (back scatter) and in solid body rotation [59]. This results in a model which is competitive with the dynamic Smagorinsky model in simulation of the turbulent mixing layer and turbulent channel flow, as we will see below.
Singular Value Models Another development is to base the eddy viscosity on the singular values i of the velocity gradient tensor, resulting in Nicoud's -model [59]. Here Another model based on singular values, but now on those of the rate-of-strain tensor S, i.e. the symmetric part of the velocity gradient tensor, is Verstappen's QR-model [96]. It will be discussed in Sect. 7, including its generalization to anisotropic grids AMD [75].
Regularization Models Regularization does not dissipate turbulent kinetic energy, but restricts the creation of turbulence by explicitly filtering of the convective terms. Several options exist: Historically [27], the first regularization (27a) was proposed by Leray in 1934 [41] as a mathematical tool to analytically study the regularity of a Navier-Stokes solution. Leray regularization has been used as an LES model for simulations of a temporal mixing layer [24]. where it is found to compare favorably with the dynamic Smagorinsky model. Another example of a regularization model used for practical LES is the Lagrange-averaged Navier-Stokes-model (27b) [20,25]. It preserves the Kelvin circulation theorem in the regularized equations.

Example: LES of Turbulent Channel Flow
Some results generated with the methods described above have been collected in this section, to give a rough indication their performance. More comparisons can be found in monographs like [66] and the vast literature on the subject. To allow comparison with the DNS methods discussed above, the earlier test case on turbulent channel flow from Sect. 5.1 is being discussed again.
Simulations of two weakly compressible channel flows at M b = 0.2 have been performed using various of the above turbulence models. They correspond to the DNSs at friction Reynolds numbers Re ≈ 180 and Re ≈ 395 by Moser et al. [57]. The set up is similar to the one in Sect. 5.1. The simulations have been performed without model, with the symmetry-preserving regularizations, with the -model [59] and with Vreman's [106] eddy-viscosity model. For simplicity, the filter width of the regularizations is chosen equal to the mesh spacing.
The bulk Reynolds number based on the half-height of the channel H is fixed at either Re b = 2800 or Re b = 6875 by prescribing a uniform body force. The simulations at the friction Reynolds number Re ≈ 180 are performed on the coarse grid D from Sect. 5.1, and the simulations at Re ≈ 395 are performed on the 64 3 grid E which stretches towards the wall. Characteristics of these grids are listed in Table 1.
The basis for all simulations is the fourth-order accurate dispersion-relation-preserving discretization of the compressible Navier-Stokes equations proposed in Sect. 4.1.3; no artificial dissipation is applied. The turbulence models ( and Vreman) are added as extra molecular viscosity; for the regularization models the appropriate filtering has been added.
For each of the simulations, time averages are recorded from 800 to 1600 time units. The computed friction Reynolds numbers are listed in Table 3. The averaged mean velocity and turbulent fluctuations normalized by the computed friction velocity are shown in Fig. 8 for Re ≈ 180 and Fig. 9 for Re ≈ 395 . The results obtained without model predict a friction Reynolds number which is slightly higher than the actual value. The normalized (using the friction velocity) mean velocity profiles obtained without model differ from results for the DNS by Moser et al. [57], especially for the simulation at Re ≈ 395 . Note that for the latter case the flow details are smaller, hence the grid is relatively coarser.
The simulations with Leray regularization lower the wall friction compared to the c 2 regularization and the simulation without model. Simulations with eddy-viscosity models are in general more dissipative and lower the friction Reynolds number compared to the regularization models. This effect seems stronger when the relative resolution of the grid is better (at lower bulk Reynolds number). At poorer resolution, the c 2 regularization and the simulation without model overpredict the wall friction. As a preliminary (and not really surprising) conclusion, it seems advantageous to reduce the influence of a turbulence model when a grid becomes finer   (i.e. better resolving). This observation is in line with the studies of the "error landscape" by Klein et al. [35]; see also [46,59].

Minimum Dissipation Eddy-Viscosity Models
Following the previous conclusions, which basically boil down to being careful with additional viscosity, we next consider eddy-viscosity models that give the minimum eddy dissipation required to prevent accumulation of kinetic energy in the LES solution. The first minimumdissipation model is the QR model [96]. This model depends on invariants of the resolved rate-of-strain tensor and switches off in laminar and transitional flow. It will be shown that it gives good results on isotropic grids, but insufficient dissipation on anisotropic grids. To address this flaw, a new minimum-dissipation model for anisotropic grids AMD [75] is proposed; the text below gives a summary of the modeling. Further considerations on minimizing eddy viscosity by means of invariants can be found in [97] and [92].

The QR Model
The QR model, introduced by Verstappen [96], is an eddyviscosity model which gives the minimum eddy dissipation required to remove sub-grid scales from the LES solution, and is based on the invariants Q and R of the rate-of-strain tensor. In the derivation of the QR model, the LES filter width is related to the computational grid from the outset. Also, the sub-grid scales are assumed to be periodic over a grid scale, and the eddy viscosity is assumed to be constant over a grid cell. The basic line of reasoning behind this model is a balance between the production of sub-grid scale energy by the convective term (related to R) and its dissipation through flow gradients (related to Q).

Derivation from the Minimum-Dissipation Condition
Minimum-dissipation models are derived from the assumption that the sub-grid scales of a LES solution should not become dynamically significant. This condition can be formalized by bounding the kinetic energy of the sub-grid scales ′ from above on each grid cell . Unfortunately, this condition cannot be used directly to derive an LES model in terms of the resolved scales, because the evolution equation for the kinetic energy of the sub-grid scales 1 2 | ′ | 2 depends on the sub-grid scales. Therefore we convert the above condition by bounding the energy of the sub-grid scales from above in terms of known quantities.
Assuming periodicity of the sub-grid scales, a Poincaré inequality bounds the energy of variations of a solution of the Laplace equation (i.e. deviations ′ from its mean) by the norm of its gradient where C is the Poincaré constant for the domain . Using (29) with = , the non-growth requirement for ′ in (28) can be transferred to a related non-growth requirement for the gradient ∇ :  where | ⋅ | 3×3 is the norm in ℝ 3×3 .
The evolution of the velocity gradient energy 1 2 |∇ | 2

3×3
can be expressed (using Einstein convention) by rewriting the Navier-Stokes equations as [96] The term − (∇ ) T ∇ ∶ S represents the creation of velocity gradient energy by the convective terms in the Navier-Stokes equations. For incompressible flow, it can be shown [10,96] that this term can be expressed in terms of the third invariant of the rate-of-strain tensor R( ) ≡ −S 2 ∶ S∕3 = − det (S) as Substitution of the eddy-viscosity assumption, and assuming that the molecular and eddy viscosity are constant over a grid cell, gives after integration Here |∇S| 2 3×3×3 ≡ k S ij k S ij is the norm in ℝ 3×3×3 , and is the outward-pointing unit normal vector along the boundary of the grid cell. In the derivation of a minimum-dissipation model, terms that can be rewritten to boundary integrals are ignored, because they model an exchange of velocity gradient energy with the environment of the grid cell, and not the creation or dissipation of energy inside the grid cell.
A minimum-dissipation model, obeying (30), delivers the minimum eddy dissipation required to cancel the generation of velocity gradient energy by the convective terms, i.e. it (just) satisfies "dissipation ≥ production". Thus, neglecting the molecular viscosity, the eddy viscosity e should satisfy The minimum amount of eddy viscosity that satisfies this condition is given by Finding a (sharp) lower bound for the denominator in (34) gives an amount of eddy viscosity that is larger than the strict minimum. Such an estimate can be obtained from (once more) the Poincaré inequality (29) by choosing = S . Introducing the second invariant of the velocity gradient Q( ) = 1 2 S ∶ S (for incompressible flow), such a Poincaré inequality reads The Poincaré constant C depends on the size of the grid cell. Assuming the flow field satisfies Laplace's equation, it is equal to the inverse of the smallest non-zero eigenvalue of the negative Laplace operator − ≡ −∇ ⋅ ∇ = ∇ T ∇ on the grid cell . Payne and Weinberger [62] have derived an analytic value of this Poincaré constant: C = 2 ∕ 2 for convex domains of diameter . We will see later (Sect. 7.1.2) that on elongated domains this constant can be decreased (i.e. sharpened). We will also see that it is relevant to choose this constant as small as possible, but not too small.
With the above estimate, which is used from right to left, from (34) an upper bound for the required minimum eddy viscosity is found: The approximate equality at the right-hand side is obtained by averaging the integrals over the grid cell and evaluating them using mid-point integration [2,96]. Finally, the QR model is obtained by clipping the numerator to positive values: The computational complexity of the QR model is low. Compared to the classical Smagorinsky model, the QR model only needs additional computation of the third invariant R( ) of the already-available tensor S. The theoretical properties of the QR model are promising. The secondinvariant Q( ) is positive, and therefore the model delivers eddy dissipation if the third invariant R( ) is positive. The above derivation demonstrates that R( ) is only positive in regions where sub-grid scales are produced.
It can be shown that R( ) vanishes in flows with zero exact eddy dissipation according to the analysis by Vreman for all possible LES filters [106]. In practice, this means that the QR model switches off for laminar and transitional flows. Also, the QR model switches off in two-dimensional flows such as solid body rotation. This is considered to be an appropriate property of an LES model [59].
Remark The QR model can be related to the Smagorinsky model by choosing the Smagorinsky constant C s in (26) according to QR e = C max {R( ), 0} Q( ) .

Choice of Filter Width on Anisotropic Grids
On anisotropic grids a choice has to be made for the filter width that appears in the Poincaré constant C . In [96] the filter width of the QR model is set according to This filter width is dominated by the mesh spacing in the fine grid direction. The results in Sect. 7.3.5 will demonstrate that this choice gives insufficient eddy dissipation to remove sub-grid scales in the coarse directions. A more conventional approximation for the filter width on an anisotropic grid is the geometric mean approximation The most robust choice (with the largest e ) is presumably to set the filter width equal to the mesh spacing in the coarse direction Below, in Figs. 17 and 18, we will present some results obtained with the proposed filter widths. This will demonstrate that it is desirable to develop a better model for computational grids that are not isotropic.

The Anisotropic Minimum-Dissipation (AMD) Model
The QR model has good theoretical properties on isotropic computational grids as we will show below. However, as mentioned above, the velocity gradient energy used in the derivation of the QR model is dominated by variations in the direction with the smallest mesh spacing, which leads to insufficient damping in the coarse direction. This flaw can be addressed by realizing that when there is sufficient resolution at a given (small) grid spacing, then in that direction there is no sub-grid scale activity to account for. This observation suggests to relate production and dissipation of sub-grid energy to the velocity gradient with respect to the grid, i.e. in computational space rather than in physical space. The anisotropic minimum-dissipation AMD model is based on this observation. Below we will present its derivation, which follows the derivation of the QR model discussed above. We will also demonstrate that this model is consistent with the gradient model on anisotropic grids.

The Derivation of the AMD Model for Rectangular Grids
In the above spirit, an alternative for the decay of sub-grid energy in (30) is obtained if the velocity gradient is scaled in each direction by the mesh spacing in that direction: where the scaled velocity gradient is denoted ∇ ij = ( x) i i v j (no summation). In this way we look at the variations of the solution in computational space rather than in physical space. Also this scaled velocity gradient energy is connected to the sub-grid scale energy (28) via a Poincaré inequality: where this time the Poincaré constant C equals the inverse of the smallest non-zero eigenvalue of the differential operator −∇ ⋅ ∇ = ∇ T ∇ on the grid cell . It is independent of the grid size; the latter has been mixed into the definition of the scaled velocity gradient energy. For a scaled Laplace operator it equals 1∕ 2 [62].
Similar to (31), the evolution equation for the scaled velocity gradient energy inside a grid cell can be expressed as where the skew-symmetry of the differential operator ( ⋅ ∇) for incompressible flow was used to rewrite the production of scaled velocity gradient energy by the convective term [55,100], just as in the derivation of the QR model. In contrast to Eq. (33), the above evolution equation expresses the production and dissipation of scaled velocity gradient energy in terms of the velocity gradients ∇ and ∇ instead of only the rate-of-strain tensor S.
Demanding that the generation of scaled velocity gradient energy inside a grid cell should (at least) be canceled by the eddy dissipation gives the condition As in the derivation of the QR model, using the scaled Poincaré inequality in (41), with replaced by ∇ , the integral on the left-hand side can be bounded from below. This then leads to an amount of eddy viscosity that is sufficient to satisfy the above criterion:

3
where the Poincaré constant is now independent of the grid size.
Introducing the symmetric tensor the eddy viscosity of the AMD model is defined by where the integrals over a grid cell are evaluated using the mid-point integration rule.
Below, in Sect. 7.2.2, it is shown that the tensor B is proportional to the leading-order term of the Taylor expansion of the exact sub-grid-stress tensor for a separable LES filter when the Poincaré constant is chosen as C = 1∕12 . With this choice, the contraction −B ∶ S is proportional to the eddy dissipation of the gradient model [12] on an anisotropic rectangular grid. Also, this value of C does not differ much from the Poincaré constant for the solution of a Laplace equation C = 1∕ 2 .

Remarks
• With some notational modifications, the above model can also be applied on curvilinear grids [71]. • The AMD model has been derived based on the equations in computational space, and not on the equations in physical space which is the usual way for turbulence modeling. This strategy falls in a wider philosophy of first discretizing the basic flow equations and only thereafter take the actions you want to take [94]. • There exists an interesting quantitative difference between the QR model and the AMD model for two-dimensional flow. The leading-order term of the exact sub-grid tensor gives no eddy dissipation for two-dimensional flow on isotropic grids, but may give eddy dissipation for two-dimensional flow on anisotropic grids. The derivation of the QR model implicitly assumes an isotropic grid, and switches off for all two-dimensional flows [96]. However, the AMD model follows the behavior of the exact sub-grid tensor, and does give eddy dissipation for certain two-dimensional flows on anisotropic grids.

Consistency with the Anisotropic Gradient Model
The dependence of the AMD model on the velocity gradient can be explained by deriving an eddy-viscosity model which gives the same eddy dissipation as the anisotropic gradient model [12]. For this model with a separable anisotropic , finite-volume filter, just averaging cell values, the Taylor expansion of the coarse-grained velocity field is Substitution hereof in the sub-grid stress tensor gives a Taylor expansion which can be expressed in terms of the scaled gradient ∇ as The leading-order term of this expansion is the gradient model [12] for a rectangular separable filter, given by where B is the symmetric tensor in Eq. (43). The eddy dissipation of this anisotropic gradient model is equal to To rigorously derive the AMD model from the gradient model, the sub-grid model should be proportional to the velocity gradient tensor instead of the rate-of-strain tensor with an eddy dissipation given by Equating this with (46) gives an eddy viscosity e equal to We recognize the eddy viscosity of the AMD model in Eq. (44) with C = 1∕12.
Thus, the AMD model gives eddy dissipation exactly if the gradient model gives eddy dissipation. This implies that the AMD model switches off for flows with zero exact eddy dissipation [106]. Hence, consistency with the leadingorder term of the exact sub-grid stress is a nice theoretical property. It is to be stressed that this consistency does not hold if −B ∶ S < 0 , i.e. the AMD model vanishes when the leading-order term of the gradient model would yield negative dissipation (and becomes unstable in practice) [105].

Discrete Corrections for the Model Constant
In most numerical methods, the derivative in the convective term as used in the numerators of the QR and AMD (45)  1 3 models is discretized using a central 2 x stencil. In contrast, the two first derivatives of the dissipative term, related to the denominator in (44) are effectively discretized on a x stencil to prevent decoupling. Another way of saying this, is that the Laplacian induced by the convective discretization has a twice larger stencil than the Laplacian induced by the diffusive discretization.
This discrete difference has to be accounted for in the above analytic reasoning. As a consequence, a correction has to be made to the length scale in the Poincaré inequality. As the length scale appears quadratically, this means that for a second-order central discretization of convection the model constant should be multiplied by a factor 2 2 = 4 . For a fourth-order discretization, a similar correction is required. It has been argued in [71] that this correction equals 2.832.
Thus, for second-and fourth-order discretization, the Poincaré constants compatible with the gradient model are adapted to If a comparison with box-filtered experimental results (instead of results filtered with the spectral filter) is to be made, then decent agreement is observed for decaying grid turbulence if the above unfiltered model constants are chosen √ 2 times as large [71]. I.e. for the second-order discretization one may choose for the box-filtered models

QR and AMD Results
To validate the proposed minimum-dissipation models, simulations of turbulent channel flow (see Sects. 5.1 and 6.3), decaying grid turbulence (see Sect. 5.2) and a temporal mixing layer are performed. The results obtained with the minimum-dissipation models are compared with results obtained with other methods described in this paper, like the Smagorinsky models [43,80] and the Vreman eddyviscosity model [106].

Turbulent Channel Flow
The studied channel flows are the classical simulations, also studied in Sects. 5.1 and 6.3, at friction Reynolds numbers of approximately 395 and 590 [57]. The bulk Reynolds number based on half the height of the channel Re b = u b H∕ is fixed at either 6875 or 10, 975 by imposing a constant bulk velocity u b through an isotropic pressure gradient. The computational grid has 64 cells in each direction, with an isotropic mesh spacing in the streamwise and span-wise directions. Table 4 lists the characteristics of the grid in wall units for both bulk Reynolds numbers.
Simulations are performed with the proposed minimumdissipation models. The results are compared with results of a no-model simulation, a simulation with the Vreman model [106], and a direct numerical simulation (DNS) on a much finer grid [57]. The minimum-dissipation models are used with the model constants proposed in Eq. (50). The grid is anisotropic, and therefore the QR model requires an approximation of the filter width. We use the filter width approximation proposed in the literature in Eq. (37) and the conventional geometric mean approximation in Eq. (38).
The friction Reynolds numbers computed in the channel flow simulations are listed in Table 5. The mean flow velocity and the turbulent fluctuations normalized by the computed friction velocity at both Reynolds numbers are shown in Figs. 10 and 11. For both channel flows, the simulations without an LES model predict a friction Reynolds number which considerably exceeds the actual Reynolds number. The Vreman model and the AMD model give good results for the studied channel flows. For these models, the error in the predicted friction Reynolds number is smaller than 3%, and the normalized mean flow velocity profiles agree accurately with the DNS.
The results obtained with the QR model are very sensitive to the approximation of the filter width. The QR model with the filter width approximation based on the geometric mean in Eq. (38), ≈ y 1∕3 , is too dissipative and gives a friction Reynolds number which is considerably smaller than the actual friction Reynolds number. The QR model with the smaller filter width approximation from the literature in Eq. (37), ≈ √ 3 y , accurately predicts the friction

3
Reynolds number and the normalized mean flow velocity profiles in channel flow simulations.

Decaying Grid Turbulence
To assess the applicability of the minimum-dissipation models to decaying turbulence, large-eddy simulations of an experiment by Comte-Bellot and Corrsin [14] are performed on a coarse 64 3 grid. The setup of the simulations is essentially the same as in Sect. 5 The energy decay and the energy spectra computed in simulations with the two model constants are shown in Fig. 12. The box-filtered reference data in the energy decay plot is obtained by fitting a velocity field on a 64 3 grid to the desired energy spectrum, and application of the boxfilter. The energy decay obtained with the AMD model with the constant for unfiltered results given in Eq. (50) agrees with the energy decay measured in the experiments. The computed energy spectra agree with the measured energy spectra at all wave numbers up to the one-dimensional pointto-point oscillation. These results indicate that the AMD model appropriately reflects the dissipative nature of subgrid scales in decaying grid turbulence.
In a similar way, the energy decay obtained with the AMD model with the model constant for box-filtered results in Eq. (51) closely agrees with the energy decay of the boxfiltered experimental measurements. The energy spectra   To check grid convergence of the AMD model, simulations of the experiment by Comte-Bellot and Corrsin [14] have also been performed on 96 3 and 128 3 grids. The model constant (50) for unfiltered results is used. The energy spectra obtained on both grids are shown in Fig. 13 and differ visibly only in the cut-off region.

Temporal Mixing Layer
To assess the minimum-dissipation models for transitional flow and for anisotropic grids, simulations of a weakly compressible temporal mixing layer at a high Reynolds number have been performed. A temporal mixing layer consists of two streams with opposite velocities. At the interface of the two streams, the velocity differences trigger a Kelvin-Helmholtz instability, which eventually causes transition to turbulence. The temporal mixing layer studied here is essentially the same as in [106]. All quantities are nondimensionalized by half the initial vorticity thickness of the mixing layer, the free stream velocity, the free stream pressure, and the free stream temperature. The non-dimensionalized coordinate x is aligned with the stream-wise direction, y with the direction normal to the mixing layer, and z with the span-wise direction. The initial non-dimensionalized velocity profile is given by a hyperbolic tangent To trigger transition to turbulence, random perturbations with an amplitude of 0.05 exp(−y 2 ∕4) are superimposed to the initial velocity field. The initial pressure is set equal to the free stream pressure p = 1 , and the initial non-dimensionalized temperature profile is set to where M denotes the free stream Mach number. The mixing layer is weakly compressible with a Mach number of M = 0.25 , and the Reynolds number based on half the initial vorticity thickness is 100,000.
The computational domain spans 90 times half of the initial vorticity thickness in each direction. The boundary conditions in the stream-wise and span-wise directions are periodic. At the boundaries in the direction normal to the u = tanh (y) , v = 0 , w = 0 . (b) Energy spectra Fig. 12 The computed kinetic energy decay and energy spectra at times 42M∕U 0 , 98M∕U 0 and 171M∕U 0 in simulations of decaying grid turbulence with the AMD model with the collocated secondorder accurate method. The model constant for unfiltered results is compared with the model constant for box-filtered results. Also, results obtained with the dynamic Smagorinsky model, the experimental measurements, and box-filtered experimental measurements are shown The simulations of the temporal mixing layer have been performed with the collocated fourth-order accurate method for compressible flow employing the proposed minimumdissipation models, with the model constants as proposed in (50). Results obtained with the minimum-dissipation models are compared with results obtained with the Vreman model [106], which is considered to be a proper model for the mixing layer. Also, the results are compared with results obtained with the classical Smagorinsky model [80]. For a comparison with the dynamic Smagorinsky model, see [75].

Mixing Layer on an Isotropic Grid
The first set of simulations is carried out on an isotropic grid with 90 cells in each direction. Figure 14 shows the evolution of the total kinetic energy in the computational domain and the momentum thickness of the mixing layer.
The results of the minimum-dissipation models closely agree with results of the Vreman model. This is desirable, because the Vreman model is known to give good results for this mixing layer [106]. The minimum-dissipation models and the Vreman model predict an approximately constant growth rate of the mixing layer after t = 60 , which suggests that the mixing layer is self-similar in the turbulent regime. The classical Smagorinsky model is overly dissipative in the transitional regime, and delays transition of the mixing layer. Also, the Smagorinsky model does not predict a linear growth of the mixing layer. These disadvantages of the Smagorinsky model in simulations of the mixing layer are sometimes resolved by using the Smagorinsky constant C s = 0.10 instead of C s = 0.17 [17]. However, this   Smagorinsky constant gives too little eddy dissipation in the turbulent regime of the mixing layer [105].
To assess the behavior of the minimum-dissipation models in the self-similar regime, the growth rate of the mixing layer and plots of the velocity fluctuations in the stream-wise direction are shown in Fig. 15. The minimum-dissipation models predict an approximately constant growth rate of the mixing layer from t = 80 to t = 160 . Plots of the stream-wise velocity fluctuations against the normalized normal coordinate computed with the AMD model indeed collapse. This demonstrates that the minimum-dissipation models appropriately capture the self-similar character of the temporal mixing layer.
The minimum-dissipation models deliver eddy dissipation for less flows (as characterized by the velocity gradient tensor) than the Vreman model. To check that the minimumdissipation models give sufficient eddy dissipation to prevent pile-up of kinetic energy, stream-and span-wise energy spectra have been recorded at the center plane y = 0 of the mixing layer. Figure 16 shows the kinetic energy spectra at t = 140 . The energy spectra closely resemble the E(k) ∼ k −5∕3 decay law: no considerable pile-up of kinetic energy is observed. In fact, the energy spectra obtained with the minimum-dissipation models closely agree with energy spectra obtained with the Vreman model.

Mixing Layer on an Anisotropic Grid
The above results demonstrate that the QR model and the AMD model give good results for the temporal mixing layer on an isotropic grid. To assess the minimum-dissipation models on anisotropic grids, the simulations from the previous section are repeated on a grid with dimensions 90 × 360 × 90 . This grid has the same mesh spacing in the stream-wise and span-wise directions, but a four-times smaller mesh spacing in the normal direction.
Simulations of the temporal mixing layer on the anisotropic grid are performed with the QR model, the AMD model, and the Vreman model. As in the previous section, the minimum-dissipation models are used with the model constants for unfiltered results in Eq. (50). On anisotropic grids, the model constant for the QR model requires an approximation of the filter width ; see Sect. 7.1.2. Here, the choice (37) gives a filter width = 1∕ √ 6 = 0.408 , which is considerably smaller than the mesh spacing in the coarser grid directions x = z = 1 . This raises the question whether the QR model gives sufficient eddy viscosity to dissipate sub-grid variations in the coarser directions for this filter width. The more conventional geometric mean approximation (38) gives a filter width of = 0.630 . The most robust choice by setting the filter width equal to the mesh spacing in the coarse direction (39) gives a filter width of = x = z = 1 for the present grid. In contrast, for the AMD model, the anisotropy of the LES filter is merged into the scaled velocity gradients that constitute the model: an ad-hoc approximation of the filter width is not required. Figure 17 shows the evolution of the total kinetic energy in the computational domain and the growth rate of momentum thickness of the temporal mixing layer. Just as on the isotropic grid, results obtained with the AMD model closely agree with results obtained with the Vreman model. On this anisotropic grid, transition of the temporal mixing layer occurs slightly earlier than on the coarser isotropic grid. This is not troublesome, because the transition to turbulence in simulations of a mixing layer is in general very sensitive to grid resolution, the accuracy of the used numerical method, and the perturbation of the initial condition [70].
The AMD model and the Vreman model predict an approximately constant growth rate of the mixing layer from (b) Span-wise spectra Fig. 16 The stream-wise and span-wise kinetic energy spectra at the center plane of the temporal mixing layer at t = 140 computed on an isotropic grid with the QR model, the AMD model, the Vreman model, and the Smagorinsky model 1 3 t = 60 until t = 160 . The results of this model are comparable to results obtained with the dynamic Smagorinsky model [75]. The results of the QR model on the anisotropic grid are very sensitive to the used approximation of the filter width. The filter width (37) does not delay the transition compared with the Vreman model and the AMD model, but gives a lower growth rate of the mixing layer in the turbulent regime. For the geometric mean filter width (38), the QR model gives results comparable to the Vreman model and the anisotropic dissipation model. The robust maximum filter width (39) delays the transition of the mixing layer compared to the Vreman model and the anisotropic dissipation model, and predicts a higher growth rate of the mixing layer in the turbulent regime.
The differences between the LES models and approximations of the filter width on anisotropic grids can further be studied by examination of the energy spectra in the turbulent regime. Figure 18 shows the stream-wise and span-wise energy spectra at the center plane y = 0 of the mixing layer at t = 140 . The Vreman model and the AMD model properly dissipate the energy of sub-grid scales on the anisotropic grid, and the obtained energy spectra closely resemble the E(k) ∼ k −5∕3 decay law.

Conclusions
In this paper, low-dissipation methods and models for the simulation of turbulent airflow have been studied. Both numerical discretization methods as well as turbulence models have been addressed. Their common theme is the  The evolution of the non-dimensionalized total kinetic energy and the growth rate of the momentum thickness in simulations of the temporal mixing layer on an anisotropic grid with the QR model with the filter width proposed in the literature in (37), the conventional geometric mean filter width in (38), and the robust maximum filter width in (39), the AMD model, and the Vreman model  The stream-wise and span-wise kinetic energy spectra at the center plane of the temporal mixing layer at t = 140 computed on an anisotropic grid with with the QR model with different filter width approximations, the AMD model, and the Vreman model 1 3 minimal use of additional diffusion, be it as part of the discretization or as part of the turbulence model.

Discretization Without Artificial Dissipation
The first part of the paper focuses on discretization methods without artificial dissipation. A new analysis of energy-conserving methods for compressible flow has been proposed using square-root variables. It concisely expresses the conservation properties of convective transport as a skew-symmetry of the convective terms. By transferring this skew-symmetry to the discrete simulation method, a class of energy-conserving symmetry-preserving simulation methods for collocated curvilinear grids can be derived. Main properties are: 1. The methods discretely conserve the primary variables mass, momentum and total energy, i.e. they are proper finite-volume methods for compressible flow. 2. The convective terms preserve the secondary variables kinetic and internal energy at the discrete level. This property ensures that artificial dissipation cannot overwhelm the eddy dissipation of an LES model. Moreover, it prevents spurious generation of kinetic energy by discrete convection and therefore improves the numerical stability of the method. 3. The energy-conserving methods can be applied for general (structured) curvilinear grids. Higher-order versions can be obtained through Richardson extrapolation. 4. The formulation in square-root variables allows for straightforward derivation of energy-conserving timeintegration methods. 5. The square-root variables also allow to define symmetrypreserving regularization turbulence models.
The developed simulation method has been assessed by performing simulations of channel flow and decaying grid turbulence. The method is stable without additional artificial dissipation for simulations of subsonic channel flow on under-resolved grids where it accurately captures the solution. Under-resolved simulations of decaying grid turbulence on coarse grids are also stable, but erroneously give pile-up of kinetic energy near the grid cut-off showing the need for a turbulence model. The symmetry-preserving method has been implemented in the simulation method Enflow of NLR [37], with which large-scale simulations of the transitional airflow over a delta wing have been performed. The challenge of these simulations is to accurately capture the development of the shear layer that separates at the leading edge of the wing. Comparison with experiments shows that the onset of unsteady flow and the onset of full turbulence are predicted accurately [71].

Low-Dissipation Large-Eddy Simulation Models
When the computational grid is not fine enough to accurately capture the smallest turbulent flow structures by means of DNS, coarse-grained turbulent flow models are required. Therefore, in the second part of the paper, an overview of low-dissipation LES models has been given. Regularization models [90,95], which generate no dissipation at all but which merely redistribute the energy over the scales, are useful on nearly-resolved grids (e.g. in channel flow simulations), but there is possibly pile-up of energy near the grid cut-off (e.g. in decaying turbulence). Some amount of energy dissipation appears to be necessary on under-resolved grids, but not too much.
The test cases show that the amount of eddy viscosity needed in an LES method is dependent on the resolution of the grid. This information forms the basis for a new class of minimum-dissipation models: the QR model [96] and the anisotropic minimum dissipation model AMD [1,75]. The former has been developed for isotropic grids, whereas the latter is a generalization to anisotropic grids. To deal with the anisotropy, the analytic derivation of the AMD method is carried out in transformed spatial coordinates similar to computational space (as opposed to physical space).
Both models are based on an estimate of the production of sub-grid energy, and give just enough dissipation to remove the sub-grid scales from the LES solution. They appropriately give eddy viscosity in regions of turbulent flow, but switch off in regions of laminar and transitional flow. The models are consistent with the exact sub-filter tensor: the QR method on isotropic grids, the AMD method also on anisotropic grids. This makes the AMD model a good starting point in the research of consistent and practical LES models.
The fundamental simplicity of minimum-dissipation models is elegant. However, their application could be limited to flows in which the energy dissipation of forward scatter down to unresolved scales is dominant. Combination with a (energyconserving) regularization model to model back scatter might then be useful. Also, the minimization ideas of AMD could be applied to other models, like the Bardina model, as done by Streher et al. [85]. A further extension could be to include non-dissipative terms in the AMD model, e.g. related to the skew-symmetric part of the velocity gradient, while staying consistent with the exact sub-filter tensor [78,79].