## Abstract

Restrictions on computing power make direct numerical simulation too expensive for complex flows; thus, the development of accurate large eddy simulation (LES) methods, which are industrially applicable and efficient, is required. This paper reviews recent findings about the leading order dissipation rate associated with high-resolution methods and improvements to the standard schemes for use in highly turbulent flows. Results from implicit LES are presented for a broad range of flows and numerical schemes, ranging from the second-order monotone upstream-centered schemes for conservation laws to very high-order (up to ninth-order) weighted essentially non-oscillatory schemes.

## 1. Introduction

As current computational power does not allow direct numerical simulation (DNS) of complex flows, large eddy simulation (LES) has emerged as a viable alternative where the time-dependent behaviour of the flow must be resolved. Conventional LES, where an explicit subgrid model is added to the averaged Navier–Stokes equations, has been employed successfully in many prototype flows; however, it is known to provide excessive dissipation in flows where the growth of an initially small perturbation to fully turbulent flow must be resolved (Lesieur & Metais 1996; Pope 2000). It has been recognized that some numerical schemes gain good results in complex flows without the explicit addition of a subgrid model (Lesieur & Metais 1996). This occurs when the subgrid model is implicitly designed into the limiting method of the numerical scheme, based on the observation that an upwind numerical scheme can be rewritten as a central scheme plus a dissipative term (see Drikakis 2003; Drikakis & Rider 2004; Grinstein *et al.* 2007 and references therein). Such implicit subgrid models fall into the class of structural models, as there is no assumed form of the nature of the subgrid flow; thus, the subgrid model is entirely determined by the structure of the resolved flow (Sagaut 2001).

Efforts towards the development of a theoretical framework, which justifies the use of the truncation error of the numerical scheme as a subgrid model and guides structured development of numerical methods, were presented in Margolin & Rider (2005) and Margolin *et al.* (2006). The first important element of such a theoretical justification is to consider the requirements on the numerical scheme dependent on how well resolved the physical problem is. By considering numerical resolution, there are two clear categories of implicit LES (ILES).

The first category relates to the behaviour of the numerical solution when the simulation is well resolved. In this case, the grid resolution and numerical methods are combined in such a way that the large scales present in the problem are fully resolved by the numerical method with negligible influence of numerical viscosity. To ensure a reasonable result, the numerical scheme must also resolve accurately the vortices that interact the strongest with the large scales, i.e. the neighbouring modes. A final requirement exists in the case of transitional flows, where the numerical method must have adequate resolution to capture the initial perturbations that become unstable and trigger turbulence in the problem (a requirement for all LES). An excellent example of an ILES which falls in the category is that of Cohen *et al.* (2002), where the fundamental flow instability is simulated using very high grid resolutions, thus isolating the scales of interest from the dissipative scales.

If these three criteria are satisfied, then the level and form of the numerical dissipation should not significantly influence the problem. The basis of this assumption is an application of Kolmogorov's hypothesis, where it was shown that the appearance of a subinertial range relies only on adequate separation between the start of the inertial range and the dissipative scales, it does not rely on the form of the dissipation. In an analogous situation, as long as there is such a separation between the scales that influence the energetic scales and the scales dissipated by the numerical method, then an implicit LES should give a correct result. It should be stated that an initial assumption in the large eddy approximation is that the flow quantities of interest are controlled by the large scales—this is not always the case (e.g. reactive flows).

This first category of ILES has a relatively simple framework, so long as the user is aware of the resolution requirements of their particular numerical scheme, e.g. how many computational cells are required to capture a single mode/vortex. This approach suggests the use of very high-order schemes that focus the dissipation into a narrow band at high wavenumbers, thus allowing an independent regime to evolve. It would also be beneficial for improving the scheme's ability to be able to resolve initial perturbations in the flow field, which then either trigger or promote turbulent fluctuations.

The second category of ILES emerges when there is an insufficient separation of the large and dissipative scales. This occurs frequently in complex industrial flows where very LES is commonplace. In this case, it is desirable that the numerical methods provide a dissipative influence on the large scales, mimicking that produced by the scales that they would normally interact with if the full spectrum was present. This is the regime where both classical LES and ILES encounter difficulties, as the simulation is now sensitive to the form of the imposed subgrid model.

At this point, there are a plethora of classical LES models that can be applied to the problem, and an equally large number of ILES numerical methods. The required form of the dissipation is usually derived using a closure assumption, such as the eddy-viscosity hypothesis, and the coefficients calibrated to give the required results in classical test cases, such as homogeneous decaying turbulence. Interestingly, studies by Meyers *et al.* (2007) have shown that the optimum coefficients for LES are often functions of the choice of grid resolution and numerical filter size. Improvements in specification of the coefficients by using dynamic procedures often give improvements only in areas of the flow where the LES model is not strictly applicable, such as near walls, transitional regions and laminar flow (Pope 2004).

A clear uncertainty in this area is the universality of the approach when applied to many different flows. For example, it is not certain to what extent anisotropy persists into the beginnings of the inertial range, which is likely to have an influence on the form of the energy spectrum, the magnitude of second-order statistics and hence the dissipation rate (Casciola *et al.* 2007). In these situations, a model that is based on the assumption of homogeneous turbulence at high wavenumbers could not be expected to be accurate. This argument is equally applicable to ILES methods. Given the uncertainties involved, this second category of LES is certainly by far the most challenging. The vast majority of ILESs fall into this category, principally due to the massive computational requirement for highly resolved simulations. This includes open cavity flow (Larcheveque *et al.* 2003; Hahn & Drikakis 2005; Thornber & Drikakis 2008), geophysical flows (Smolarkiewicz & Margolin 1998; Margolin *et al.* 1999), delta (Gordnier & Visbal 2005) and highly swept wings (Hahn & Drikakis 2008) and low resolution (less than 128^{3}) decaying turbulence (Fureby *et al.* 1997; Porter *et al.* 1998; Fureby & Grinstein 2002; Margolin *et al.* 2002; Drikakis *et al.* 2006; Hickel *et al.* 2006; Thornber *et al.* 2007).

The theoretical analysis of incompressible flow simulations using modified equation analysis (MEA), particularly the influence of finite volume schemes, has been investigated by Drikakis & Rider (2004), Margolin & Rider (2005), Margolin *et al.* (2006) and Grinstein *et al.* (2007). These papers emphasized the importance of the finite volume approach in that it allows the natural evolution of volume averaged quantities, i.e. the filtered quantities necessary in LES. Numerical methods for incompressible flows can be fairly flexible in that they do not require strict monotonic behaviour to give numerical stability. Often, favourable numerical methods for turbulence do not conserve the positivity of passive scalars—something that requires additional monotonic constraints. However, the flexibility allows for a greater scope of optimization of the numerical stencil to provide favourable dissipative properties. An excellent example of this is the approximate deconvolution methods proposed by Hickel *et al.* (2006), where the weightings for the variable reconstruction are chosen to give behaviour that is as close as possible to a spectral eddy viscosity.

In compressible methods there is significantly less room for manoeuvre. The numerical method must be monotonic in the thermodynamic quantities (density and pressure) in order to remain stable, i.e. so that the simulation does not crash. This stability is usually added as explicit diffusion terms, such as in the Jameson scheme (Jameson *et al.* 1981), or through upwinding of the fluxes, as employed in the Godunov method (Godunov 1959). This is especially important in the simulation of multi-component flows where overshoots in the mass or volume fraction can lead to the prediction of negative densities or extremely unphysical equations of state (e.g. ratios of specific heats less than one for a perfect gas equation of state).

It has been highlighted that upwind methods are typically overly dissipative for simulations of homogeneous decaying turbulence when using second- and third-order methods (Garnier *et al.* 1999). Resorting to very high-order methods improves the situation; however, it does not cure the problem completely (Thornber *et al.* 2007). Higher order methods will give reasonable results in terms of growth of the length scales and kinetic energy dissipation rate; however, the flows are dissipative at high wavenumbers, leading to spectra that fall off sharply (Thornber *et al.* 2007).

Analysis of the source of the dissipation of turbulent kinetic energy (TKE) in upwind schemes reveals a remarkably simple relationship. The absolute dissipation of fluctuating kinetic energy is proportional to the temperature multiplied by the change of entropy (assuming an approximately isothermal flow; Thornber *et al.* 2008*a*). This neglects the ‘apparent’ dissipation of kinetic energy caused by isentropic transformation of kinetic energy to internal energy in the form of local compressions and expansions. Using MEA, the evolution of entropy can be derived for a given compressible numerical scheme. This method has been developed in Thornber *et al.* (2008*a*) and demonstrates that the overly dissipative behaviour observed in simulations of homogeneous decaying turbulence can be ascribed to numerical dissipation that is proportional to the speed of sound. This excess dissipation delays the appearance of a subinertial range until much higher grid resolutions.

Taking this into account, the reconstruction method can be modified to ensure that the dissipation rate becomes constant at low Mach (Thornber *et al.* 2008*b*). The key element to removing the adverse effects of low Mach number dissipation is to progressively central difference the velocity components as the Mach number tends towards zero, which removes the speed of sound dependence of numerical dissipation, as shown theoretically and numerically in Thornber *et al.* (2008*b*).

What remains is a dissipation rate proportional to the velocity increment cubed divided by the cell length in a similar functional form to the required dissipation rate for a turbulent subgrid model (∝*u*^{3}/*l*), which is constant in Mach. The modified numerical method is still monotone; hence, it maintains the positivity of scalar quantities and has no noticeable difference in stability. The following sections will examine the influence of the order of accuracy and implementation of the improved version of the numerical method for four test cases ranging from the canonical case of homogeneous decaying turbulence to computational aero-acoustics of cavities, flow over a swept wing and shock induced turbulent mixing. It will also highlight the benefits from moving to higher order of accuracy numerical methods as compared to the classical second-order upwind methods.

## 2. Homogeneous decaying turbulence

Turbulence remains the greatest challenge in fluid modelling. Many practical applications involve low Mach number turbulence, so it is important to test the ILES approach in this respect. Much work has been carried out both experimentally and numerically (using DNS) on simulating the decay of a homogeneous cube of turbulence Orzag & Patterson (1972), against the various theories that have been developed regarding turbulence—notably the well-established work of Kolmogorov (1941, 1962). Results have been mixed, with no real consensus on many key parameters. It is, therefore, instructive to examine where ILES falls within the range for a number of these measures of flow and how the improved numerical method affects these results.

The problem has been simulated in line with previous work (Thornber *et al.* 2007). In this section, the TKE spectra are presented, as this is usually a weakness in ILES using shock capturing methods. It should be noted, in addition to the results presented here, that analysis of the growth of the integral length scales and dissipation of TKE shows that they are within the bounds expected from experiment and theory.

Figure 1 shows the three-dimensional kinetic energy spectra for standard and modified Godunov methods as a function of wavenumber *k*. The abbreviations ‘MM’, ‘VL’ and ‘M5’ stand for second-order Minmod limiting van Leer limiting and fifth-order accurate monotone upstream-centered schemes for conservation laws (MUSCL), respectively. The addition of ‘+LM’ implies that the standard method has be modified using the low Mach correction. Here, we can now distinguish between the general dissipation related to the limiting method and the much more radical physical change in the structure of the flow gained by correcting for the low Mach number dissipation. Kolmogorov's theory suggests the inertial range of the cascade should follow *k*^{−5/3}, as represented in the results by the thick-dashed line. On this relatively coarse grid, the flow is conventionally under-resolved with the more dissipative methods; however, the improved method allows for small scales to be captured, giving a much better correlation with theory over the entire flow.

## 3. Deep, open cavity flow

Forestier *et al.* (2003) investigated the flow over an open cavity with a length to depth ratio of 1 : 2.4 at a Reynolds number based on a cavity length of 860 000 and Mach 0.8. They measured mean and fluctuating velocity components and the sound pressure level spectrum (SPL) at a single point located at *z*/*L*=−0.7. This case is ideal for the validation of compressible LES methods, as it includes a strong shear layer, low Mach number perturbations within the cavity, and strong acoustic waves. The fundamental acoustic frequencies and their harmonics can usually be described using an empirical formula proposed by Rossiter (1964). However, in this case, the fundamental frequencies are not described well by Rossiter's theory without modifying the coefficients away from the recommended values.

The results presented in this section were generated using a finite volume Godunov method employing the improved fifth-order in space reconstruction method (Thornber *et al.* 2008*b*) in conjunction with the Harten–Lax–van Leer-contact (HLLC) approximate Riemann solver (Toro 1997). Time stepping is achieved using a third-order accurate Runge-Kutta method (Spiteri & Ruuth 2002). Three grids were employed of 0.8, 1.4 and 3 million grid points, respectively, as used in Thornber & Drikakis (2008). An additional simulation was conducted at the coarse grid resolution using the standard fifth-order MUSCL scheme. It should be stated that, while the modified reconstruction method improves the resolution of fine scale features within the cavity and shear layer, the large coherent structures themselves are driven by a relatively high Mach number flow. Hence, the standard fifth-order method also gains good results in terms of mean and fluctuating velocity components, as it still captures the large eddies with reasonable accuracy on these grids. Figure 2 shows flow visualizations consisting of iso-surfaces of ‘*Q*’ criteria (Jeong & Hussain 1995),(3.1)where *u*_{i} are the Cartesian fluid velocities and *x*_{i} the Cartesian coordinates, which highlight the turbulent nature of the flow. Figure 2*b* shows the predicted and experimental sound pressure levels. Both sets of results demonstrate that the numerical approach employed here adequately represents the turbulent flow physics for both the mean flow velocities and Reynolds stresses.

Figure 3 shows typical results for the mean and fluctuating velocities over the shear layer for the three grid levels using the improved reconstruction method at *x*/*L*=0.6 (where *x*=0 is at the leading edge of the cavity). Comparisons with experiment at *x*/*L*=0.05, 0.2, 0.4, 0.8 and 0.95 give a similar level of agreement.

Indeed, the discrepancies present are due to the initial specification of the boundary layer at the inlet, which is slightly too large. This leads to a thicker boundary layer at the leading edge of the cavity, thus increasing the height of the centre of the shear layer, as can be seen from the points of maximum Reynolds stress in figure 3. Despite this, the results for appear to be improving at finer grids, implying that, to some extent, the Reynolds stresses depend on the accuracy of the representation of flow structures in the incoming boundary layer. The sound pressure levels are predicted to within 6 dB for all grid resolutions and within 2 per cent of the experimentally measured frequency.

## 4. Swept-wing flow

Swept and delta wings can be found in all modern aircraft travelling at transonic or supersonic speeds. At present, there are no theoretical models capable of predicting the characteristic leading edge separation with any degree of certainty, nor can the nature of the leading edge vortex breakdown observed in swept-wing flow be convincingly explained. The complexity of the flow field makes it also extremely challenging to predict with any numerical method.

The flow around a swept wing at an angle of attack of 9° has been simulated with the present ILES approach using high-resolution methods for a Reynolds number of 210 000 and a near incompressible Mach number of 0.3 (Hahn 2008). The computational grid employed comprises 12.7 million points, which is approximately half of that used in a conventional hybrid Reynolds-averaged Navier–Stokes (RANS)–LES of the identical case (Li & Leschziner 2007). Please note that the mesh employed in the ILES is nearly wall-resolved, *y*+ ranges from 1 to 5, whereas the grid for the hybrid RANS–LES yields values between 20 and 40. Furthermore, third-order accuracy in space and time integration have been obtained by a modified MUSCL reconstruction method (Zoltak & Drikakis 1998) and a third-order Runge–Kutta scheme with an extended stability region (Drikakis & Rider 2004), respectively.

The general physics of this problem has been captured in the ILES, as can be demonstrated by the instantaneous flow visualization shown in figure 4. The flow topology reveals separation along the leading edge of the swept wing, subsequent formation of the characteristic vortex roll-up and vortex breakdown near the wing tips. Moreover, the time-averaged velocity and TKE profiles in the fully separated, turbulent region near the wing tip (90% half-span and 50% local chord), figure 5, provide insight into the quality of the ILES. Here, data from the LES using the third-order high-resolution method (ILES; Hahn & Drikakis 2008) and from a hybrid RANS–LES (Li & Leschziner 2007) are compared against experiments (S. Zhang 2006, personal communication). Despite using more than 24 million grid points, the hybrid RANS–LES fails to predict the flow separation. Consequently, the results for the TKE disagree strongly with the experimental measurements. The velocity obtained by ILES, on the other hand, is nearly identical to the experiment. Although the magnitude of the TKE is slightly lower, the high-resolution method is able to predict the correct shape of the profile. Future work could be directed towards designing a hybrid RANS–ILES strategy (Tucker 2004) that enables relaxation of the stringent grid requirements in the near-wall region, while maintaining the accuracy of the simulation.

## 5. Turbulent mixing experiment: double bump

The ILES approach has also been applied to complex multi-component turbulent mixing problems. The test case presented in this section corresponds to the Atomic Weapons Establishment's shock-tube experiment featuring a Richtmyer–Meshkov instability arising from an initial large-scale two-dimensional double-bump perturbation at the interface between two gases (Holder *et al.* 2003).

Various reconstruction methods ranging from second to ninth order have been investigated in conjunction with the HLLC approximate Riemann solver (Toro 1997) and a third-order accurate total variation diminishing Runge–Kutta time-marching algorithm (Drikakis & Rider 2004). Spatial high resolution is achieved by using one of the following standard schemes: second-order MUSCL van Albada (VA; Toro 1997), third-order extended MUSCL van Albada (M3; Zoltak & Drikakis 1998), fifth-order MUSCL (M5; Kim & Kim 2005), ninth-order weighted essentially non-oscillatory (W9; Balsara & Shu 2000) or their improved counterparts VALM, M3LM, M5LM, W9LM, incorporating the low Mach modification (Thornber *et al.* 2008*b*). Furthermore, a five-equation model has been employed for representing the two different species (Allaire *et al.* 2002). The linear shock-tube problem is discretized on a Cartesian grid comprising 160×80×40 cells, which corresponds to a uniform grid spacing of *Δ*=1/4 cm. A successive grid refinement study up to *Δ*=1/16 cm has also been performed, but will not be discussed in detail.

The initial condition, shown as a two-dimensional slice in figure 6, is given by a SF_{6} region (dark) encased in air (light) with the large-scale perturbation on the right-hand interface. As the incident shock wave of Mach 1.26 passes from left to right, the membranes separating the two fluids rupture and the gases begin to mix. The shock wave reflected from the end wall decelerates the SF_{6} layer and a high degree of turbulent mixing can be observed; see the three-dimensional structures in figure 6 representing constant values of volume fraction.

Results from the present simulations can be compared qualitatively with the experiment (Holder *et al.* 2003), see figure 7. Here, the experimental data is represented by the black lines overlaid on top of the two-dimensional averaged volume fraction and SF_{6} density contours obtained with the improved fifth-order method (M5LM) at identical times. As shown in figure 7, the air/SF_{6} interface and the shock position are in very good agreement with the experiment. Plots from other reconstruction methods have been omitted because they all yield similar results, with minor differences regarding the location and the shape of the spike. However, the results obtained with different methods seem to converge as grid resolution improves.

The large-scale flow development in the simulations can also be validated quantitatively against experiment. For this purpose, the location of the bubble and the spike is presented in figure 8*a* as the distance between the end plate of the shock tube and the characteristic interface position. The interface position is given by the left edge and the right edge of the large bubbles and the tip of the spike, respectively. In accordance with the previous statement, all methods yield locations in good agreement with the large-scale features observed in experiment. Moreover, regardless of the reconstruction employed, the low Mach modification leads to improved results for the late-time spike position, i.e. *T*>2.5 ms. The grid refinement study revealed that the locations predicted by all methods seem to converge to the experimental values.

The development of integral TKE provided in figure 8*b* gives further evidence for the improvement at later times due to the low Mach modification. Its effect is reflected clearly by the increase in turbulent energy observed for all methods when compared against their standard counterparts. More specifically, the modified third-order method (M3LM) seems to result in a late-time TKE level close to the one produced by the standard fifth-order scheme (M5). A similar observation can be made for the modified van Albada method (VALM). With respect to grid refinement, TKE seems to converge to a level slightly above the W9LM data in figure 8*b* at late times.

## 6. Conclusions

This paper has discussed the theoretical framework justifying the use of the ILES method for complex turbulent flows. It has first highlighted the importance of choice of appropriate numerical methods that can capture shocks, track multiple materials and have minimal dissipation of high wavenumber features in turbulent flows. It has been demonstrated that the ILES methodology can capture canonical flows such as homogeneous decaying turbulence with the appropriate kinetic energy spectrum and kinetic energy decay rates. Simulations of more complex geometries show that the methods can be applied to aeroacoustics problems, swept-wing configurations and shock-induced turbulent mixing, providing excellent comparisons to experiment, especially considering high Reynolds numbers and complex flow physics.

## Acknowledgments

The authors would like to acknowledge the financial support from the EPSRC, Ministry of Defence and Atomic Weapons Establishment through the EPSRC (EP/C515153)-JGS(no. 971) project; EPSRC in relation to computational resources (project EP/D053994/1); and EPSRC, BAE SYSTEMS and Ministry of Defence through the Defence Aerospace Research Partnership consortium ‘Modelling and Simulation of Turbulence and Transition for Aerospace’.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Applied large eddy simulation’.

- © 2009 The Royal Society