## Abstract

Since its introduction, in the early 1970s, large eddy simulations (LES) have advanced considerably, and their application is transitioning from the academic environment to industry. Several landmark developments can be identified over the past 40 years, such as the wall-resolved simulations of wall-bounded flows, the development of advanced models for the unresolved scales that adapt to the local flow conditions and the hybridization of LES with the solution of the Reynolds-averaged Navier–Stokes equations. Thanks to these advancements, LES is now in widespread use in the academic community and is an option available in most commercial flow-solvers. This paper will try to predict what algorithmic and modelling advancements are needed to make it even more robust and inexpensive, and which areas show the most promise.

## 1. Introduction

Turbulence occurs frequently in nature and technological applications and has been the subject of study for several centuries. In 1510, Leonardo da Vinci accompanied a drawing of the vortices shed behind an immersed, blunt obstacle with the following note:
Observe the motion of the water surface, which resembles that of hair, that has two motions: one due to the weight of the shaft, the other to the shape of the curls; thus, water has eddying motions, one part of which is due to the principal current, the other to the random and reverse motion (personal translation.)

Unfortunately, the mathematical tools were not available to enable him to follow up on this comment. Several centuries had to pass before Osborne Reynolds would formulate the decomposition foreshadowed by Leonardo's observations. Since then, however, much progress has been made in the study of turbulent flows. Statistical theories of turbulence have provided understanding of the scaling laws in various flow regimes, and experiments have given insights on the statistics and structure of turbulent flows. A significant contribution to this progress is the development of numerical solution techniques for turbulent flows, in which the dynamics of some or all of the turbulent eddies are captured by the simulation. We shall refer to them as ‘eddy-resolving’ methods; they include direct numerical simulations (DNS) and large eddy simulations (LES). While in DNS, all the turbulent motions are resolved accurately (i.e. the grid is finer than the smallest turbulent eddy), in LES only the contribution of the large, energy-carrying structures to momentum and energy transfer is computed exactly, and the effect of the small scales of turbulence is modelled. Since the small scales tend to be more homogeneous and universal, and less affected by the boundary conditions than the large ones, there is hope that their models can be simpler and require fewer adjustments when applied to different flows than similar models for the Reynolds-averaged Navier–Stokes (RANS) equations.

Eddy-resolving techniques have made available data that had never been measurable previously: multi-point, unobtrusive measurements of velocity, velocity gradients, pressure, passive scalars, etc. Databases obtained from simulations have provided new insights on the physics of turbulent flows. As prediction and control of turbulence becomes increasingly important, the need for accurate models has become more critical, and eddy-resolving simulations are used to make inroads into these problems. LES and derived methods, in particular, have been applied to a significantly wider range of applications than DNS: most commercial codes include an LES option, often together with a hybrid RANS/LES mode; increasingly, industrial users are identifying problems in which turbulence models are unable to predict the flow reliably. As computer power increases, and the algorithms available improve, LES will continue its transition from the academic environment towards the R&D laboratories. The transition, however, is not painless: apart from the computational resources needed, the advanced level of competence required to run LES is an obstacle to its widespread application. Grid generation, modelling, definition of the boundary conditions depend critically on the skills of the user: LES are not plug-and-play.

The purpose of this article is not to provide a complete overview of LES, but to highlight some of the issues that are likely to drive the development of this technique over the next 15 years. In §2, the fundamentals of LES will be outlined. Then, those challenges that are most pressing, in the writer's opinion, and some possible routes to meet them will be presented. Concluding remarks on the outlook of LES will end the paper.

## 2. Formulation

### (a) Governing equations

LES are based on the use of a filtering operation: a filtered (or resolved or large scale) variable, denoted by an overbar, is defined as [1]
2.1
where *D* is the entire domain and *G* is the ‘filter function’. The size of the smallest resolved eddies is related to the length scale of the smoothing operator, the ‘filter width’, which will be denoted here by Δ. The grid size, *h*, should be sufficiently fine to allow eddies of size Δ to be represented accurately. Thus, the filter width Δ plays a central role in the modelling of the unresolved stresses, both from the point of view of the method formulation and from that of its practical implementation.

If the filtering operation (2.1) is applied to the governing equations, one obtains the filtered equations of motion. For an incompressible flow of a Newtonian fluid, they take the form:
2.2
These equations govern the evolution of the large, energy-carrying scales of motion. The effect of the small scales appears through a residual stress term, , that must be modelled. *τ*_{ij} is often called the ‘subgrid-scale’ (SGS) stress, owing to the usual adoption of a filter width that is proportional to the grid size. A more appropriate definition is ‘subfilter-scale’ (SFS) stress [2], which will be used in this paper.

### (b) Subfilter-scale models

The similarity of the small scales, which only transfer energy to smaller eddies (energy cascade), and the fact that the dissipation is set by the large scales, are exploited by the models for the SFS stresses. Their main purpose is to reproduce the energy transfer accurately, at least in a statistical sense, mimicking the drain associated with the energy cascade. It may not be necessary for a model to represent the exact SFS stresses accurately at each point in space and time, but only to account for their global effect. Thus, most subfilter-scale models in use presently are eddy-viscosity models of the form
2.3
that relate the subfilter-scale stresses *τ*_{ij} to the large-scale strain-rate tensor . In most cases, the eddy viscosity *ν*_{T} is obtained algebraically to avoid solving additional equations and increase the cost of an already expensive calculation. Moreover, since the small scales tend to be more homogeneous and isotropic than the large ones, it is hoped that even simple models can describe their physics accurately. Finally, since the SFS stresses only account for a fraction of the total stresses, modelling errors should not affect the overall accuracy of the results as much as in the standard turbulence modelling approach.

The earliest attempts to parametrize the eddy viscosity [3,4] used the filter width (in practice, the grid size) as length scale, and determined the time scale by assuming that the SFS are in equilibrium, adjusting instantaneosuly to changes in the large-scale field. The model constant was found by integrating the turbulence spectrum, assuming inertial-range dynamics. The model thus derived failed in the presence of shear, or near a wall, and *ad hoc* adjustments (decreasing the model constant, adding damping or intermittency functions) had to be included to correct for non-equilibrium effects.

Newer models better exploit the fact that the LES contain more information than RANS turbulence model, and the spectral information contained in the large-scale field can be used to model the SFS eddies more accurately. Approaches of this type include non-eddy-viscosity models [5,6], or more advanced eddy-viscosity models. Examples of the latter are dynamic models (in which the coefficient is adjusted based on the spectral content of the smallest resolved eddies [7–9]), or the structure–function model [10,11], in which the eddy viscosity also strongly depends on the smallest resolved scale.

### (c) Resolution requirements

The first complete analysis of grid resolution requirements for LES of turbulent boundary layers was performed by Chapman [12]. He considered the inner and outer layers separately; in the outer layer, in which the integral scale is proportional to the boundary layer thickness *δ* (i.e. the important eddies scale like *δ*) he estimated that the number of points required to resolve a domain of dimensions *L*_{x}×*δ*×*L*_{z} is proportional to *Re*^{0.4}. The resolution of the inner layer is much more demanding: its dynamics are dominated by quasi-streamwise vortices [13] whose dimensions are constant in wall units (i.e. when normalized with the kinematic viscosity *ν* and the friction velocity *u*_{τ}=(*τ*_{w}/*ρ*)^{1/2}, where *τ*_{w} is the wall stress and *ρ* the fluid density). If the inner-layer eddies are resolved, the number of points required for the viscous sublayer is (*N*_{x}*N*_{y}*N*_{z})_{vs}∝*Re*^{1.8}. Choi & Moin [14] revise this estimate upwards, obtaining that (*N*_{x}*N*_{y}*N*_{z})_{vs}∝*Re*^{13/7} for wall-resolved LES, and ∝*Re*_{Lx} (where *Re*_{Lx} is the Reynolds number based on the length of the computational domain) for wall-modelled ones.

In LES of free shear layers the size of the large eddies, and hence the number of grid points required, is usually proportional to *Re* [2]. Simulations of massively separated flows may also be more approachable, since the large eddies have length scales that are dictated by the geometry and are not strongly dependent on *Re* (the flow up to separation, however, obeys the standard scalings).

## 3. Challenges

We have briefly summarized the main characteristics of LES. This section will focus on some of the challenges that need to be faced to advance this technique over the next 15 years. First, those associated with the computational cost of LES will be discussed. Then, some of the improvements that are most likely to result in more widespread and accurate application of LES will be presented. The focus is on three areas: first, the development of accurate algorithms for the simulation of flows in complex geometries. Second, the emergence of hybrid RANS/LES methods, their benefits and challenges and the modelling of the region near the solid boundary. Finally, we will discuss some modelling issues that have been, so far, neglected, and need to be addressed to make LES complete and less dependent on the user's competence. Space constraints limit us from including other important issues (numerics, boundary conditions, etc.). Some of them, however, are discussed elsewhere in this volume.

### (a) Computational cost

The first issue to be faced is that of the computational requirements of the technique. We are approaching the limits of Moore's Law, and there is no reason to think that the available computational power will continue to grow unboundedly, always exceeding the power required. It is important, thus, to identify the applications for which LES is best suited, and the areas in which computational savings can be achieved, or new algorithms can give improvements.

First, we need to reiterate the very significant cost of simulations in which the wall layer (the region closest to a solid boundary) is resolved. Chapman's [12] and Choi & Moin's [14] estimates reported above show that, even if Moore's Law were to hold for the next 15 years, wall-resolved LES remain a tool for laboratory-scale experiments. For *Re*>*o*(10^{5}), over 90% of the grid points are needed to resolve less than 10% of the computational domain. Spalart [15] estimates that LES of a full aircraft at flight Reynolds numbers will not be possible until 2070; Pope [16] estimates that the number of streaks on the wing of a Boeing 777 is of order 10^{8}: even their marginal resolution is not feasible. Even wall-modelled LES (WMLES), in which the streaks are not resolved and the entire wall layer is modelled, will not be feasible, for full aircraft configurations, within a 15-year time span [15].

These considerations point out that the most appropriate area of applications may be free shear layers at high Reynolds numbers [16]: in these problems, mass, momentum and energy transport are driven by the large scale of motion; as long as the filter width is well into the inertial subrange, accurate results can be obtained. There are, however, other areas in which LES is contributing to the understanding of the physics and development of improved models. As Moin & Mahesh [17] point out, achieving real-life Reynolds numbers is not always necessary for a simulation to yield useful information; while some phenomena depend critically on the Reynolds number (local isotropy, for instance), others can be profitably studied at lower Reynolds numbers. Massively separated flows, in which the size of the eddies depends on the geometry, and not on the Reynolds number, are also a profitable target area. Problems in which the Reynolds number is low or moderate (the flow in the low-pressure turbine of a turbofan engine) may be better suited for LES, but present their own set of challenges: although modern SFS models can take into account the changes in the spectrum during transition, the accuracy of LES decreases when intermittency becomes important [18]; even for DNS, moreover, the flow during the nonlinear breakdown depends strongly on grid resolution [19]. Problems in which the inflectional instability is the driving mechanism might be better targets for LES.

### (b) Complex geometries

As LES are applied to more complex geometries, curvilinear, adaptive and unstructured meshes become more widespread. Their use may introduce errors due to the non-commutativity of filtering and differentiation in non-uniform meshes [20,21], and to the difficulty in enforcing discrete conservation of energy. Unstructured methods have become increasingly popular in recent years, and some effort has been made to develop numerical schemes that conserve kinetic energy in the discrete sense [22,23]. Some evidence [24] indicates that the conservation properties of a scheme are more critical than its order of accuracy for LES.

Grids with local refinement (of the type illustrated in figure 1) may be needed in problems in which the physical length scales vary significantly (perhaps by orders of magnitude) over the domain. Consider for example a plane channel: the requirement to resolve the streaks results in meshes that, for *Re*_{τ}>500, are excessively fine in the channel centre, resulting in wasted resources. Embedded meshes were used by, for instance, by Kravchenko *et al.* [26], who were able to reach values of *Re*_{τ} higher by a factor of 25 than those achievable at that time using a structured grid. However, the sudden discontinuities of the grid spacing expected in locally refined meshes of the type illustrated in figure 1 have more significant consequences in LES, in which the grid is marginal and significant energy may be present at the cutoff wavenumber, than in DNS or in RANS solutions, in which the flow is smooth on the grid scale. In LES non-commutativity errors may be significant, and the SFS model may give unphysical jumps in the eddy viscosity if the filter width is proportional to the grid size. Recent studies of LES on discontinuous grids [27,28] indicate that the explicit filtering of the nonlinear term is beneficial when coupled with a filter width that does not depend on the grid size, to ensure that the size of the eddies crossing a grid discontinuity is large enough that they can be resolved on the coarse side, thus reducing numerical errors.

An alternative method to compute flows in complex geometries, and one that is particularly attractive in cases with moving or deforming boundaries, is the ‘immersed-boundary method’ (IBM), in which the boundaries of the body do not conform to the grid, and the governing equations are discretized on fixed meshes. The boundary conditions must then be applied on interfaces (possibly moving) internal to the grid, in such a manner that the accuracy and efficiency of the structured solvers is maintained. Reviews of these methods can be found in [29,30]. An example of such technique is shown in figure 2, a flow visualization in an accelerating boundary layer over sand-grain roughness [31], which is modelled by the IBM using the technique proposed by Scotti [32]. The eddies are visualized as isosurfaces of the second invariant of the velocity-gradient tensor
3.1
normalized by *u*_{τ} and *ν*. Significant differences in the turbulent structures can be seen as the roughness size is varied, which result in different flow statistics. Using the IBM, realistic roughness (such as that due to rust or erosion) can be simulated, and the three-dimensional information yielded by the LES can complement the experimental data (at higher Reynolds numbers), leading to improved models for industrial and geophysical applications over real-life geometries.

The IBM is rapidly enabling the simulation of geometries with a level of complexity that could not be achieved 10 years ago. In some cases (studies of animal motion, problems involving fluid/structure interactions, etc.), IBMs may be the only feasible way to treat the moving boundaries. The algorithms used by these methods are constantly being improved. Coupling these methods with local, adaptive mesh refinement may result in large savings in terms of computational cost and will certainly become more popular over the next 15 years.

### (c) Hybrid Reynolds-averaged Navier–Stokes large eddy simulation methods

One of the developments that have made LES more suitable to applications in industrial problems is its hybridization with methods for the solution of the RANS equations. The first technique of this type was the detached eddy simulation (DES), proposed by Spalart *et al.* [33], who pointed out that turbulence models for the RANS equations predict the mean velocity in thin, attached shear layers quite accurately. In regions of separated or recirculating flow, however, their predictions incur much larger errors, beyond what is acceptable in engineering practice; there, however, the eddies are much larger, and their size is determined by the geometry and is not strongly *Re*-dependent. They, therefore, suggested the use of a model that switches seamlessly from a RANS turbulence model, to one more suitable for LES (i.e. one in which the dissipation is much smaller, and the eddy viscosity depends on the grid size). In the original proposal, the attached boundary layer was modelled using the Spalart–Allmaras [34] model, whose length scale was decreased (and made proportional to the grid size) in the detached-flow region. Several variants of the original formulation are reviewed by Spalart [35].

Figure 3 illustrates the turbulent eddies predicted by this method, compared with standard RANS models. The flow is much more three-dimensional in the DES, compared to the three-dimensional unsteady RANS, whose spanwise scales do not appear to be related to the physics of the flow (and the solution is sensitive to the spanwise dimensions of the domain). In DES, however, we observe the characteristic braid vortices, as well as a range of small scales that become smaller as the grid is refined. Spalart [35] points out that a certain level of empiricism is present, especially in some of the modifications of the original formulation [33], and that it is difficult to determine a filtering operator implicitly implied by the model; furthermore, the method depends critically on the ability of the RANS approach to predict the separation, which may be inaccurate in flows with adverse pressure gradient, or unsteady separation (the flow over contoured bumps, for instance).

A number of other hybrid techniques have been applied to a variety of flows; space limitations prevent us from discussing them; a recent review can be found in [36]. Generally, they have been most successful in massively separated flows. Despite some shortcomings, hybrid methods are beginning to be applied in industrial R&D; their ability to predict accurately the transport due to the largest eddies results in reasonably accurate prediction of the aerodynamic noise, a subject of particular interest in the aerospace and automotive industry, and of the unsteady forces in massively separated flows; the cooling of electronic components (data centres, for instance) could be another important area of application for hybrid methods.

### (d) Wall-layer modelling

One key feature that determines the accuracy of LES is whether, and how well, the momentum- and energy-carrying eddies (of size of the order of the integral scale) are resolved, and the main limitation to this methodology comes from situations in which their accurate resolution becomes prohibitively expensive. Near solid walls the size of the eddies that govern momentum transfer depends strongly on the Reynolds number and resolving those eddies requires meshes approaching those of DNS. These scales of motion may be modelled, but a strong dependence on the model chosen is then introduced.

This limitation was recognized from the earliest applications of LES, which bypass the inner layer and model its effects in a global sense (wall-modelled LES, or WMLES). If the grid near the wall is too coarse to resolve the eddies that contribute to the momentum transport, and no-slip conditions are applied, an incorrect velocity profile will result, associated with the under-prediction of the wall stress. To correct it, one must model the transport of momentum by the inner-layer scales either implicitly, by relating the outer-layer velocity to the wall stress through an assumed velocity profile, or explicitly, by parametrizing their effect in the Reynolds-averaged sense. The critical assumption is that the near-wall grid is so coarse that it contains a very large () number of near-wall eddies, and that the time-step is much larger than their time scale, so that the filtering operation becomes equivalent to long-time averaging.

Early calculations [37,38] used approximate boundary conditions similar to the wall functions applied in RANS simulations with turbulence models. These boundary conditions assume the existence of an equilibrium layer in which the stress is constant, which results in a logarithmic velocity profile that can be used to relate the velocity in the outer layer to the wall stress. They are still in widespread use, especially in environmental and geophysical flows, in which the geometry is generally quite simple, and the Reynolds number extremely high; corrections can be made for buoyancy, roughness and other effects common in these applications. In the presence of strong favourable or adverse pressure gradients, in separated flows or in flows in which the mean velocity is highly three-dimensional, hybrid methods have found application to WMLES as well. Several strategies can be used to switch between RANS in the wall layer and LES in the outer part of the boundary layer, such as changing the length scale in the model from a RANS mixing length to one related to the grid size [39], assigning the RANS and LES zones *a priori* [40] or using a blending function to merge the SGS and RANS eddy viscosities [41,42]. A more complete review can be found in [36].

The numerical treatment of the near-wall layer introduces several sources of error in an LES. The three main ones are the inaccuracy of the wall-layer modelling assumptions, the inapplicability of the SFS models in the near-wall region, and the increased numerical errors that appear when the resolution becomes marginal, or the grid is suddenly refined near a solid boundary. These errors interact with each other, but a global treatment is difficult. Another important issue is the fact that, near the wall, the integral scale of the flow is proportional to the distance from the solid boundary, and Δ*y*=*y*. Thus, the basic assumption that the grid size is fine enough to resolve the integral scale of the flow, fails: near the wall the SFS model is required to resolve all the scales of motion, i.e. to behave like a turbulence model for the RANS equations. Hybrid methods do not suffer from this problem, since the RANS model is used to parametrize all the scales of motion in the near-wall region. A new source of error is introduced, however, at the interface between the RANS and LES regions, where the disparity between the long time- and length scales imposed by the inner-layer RANS and the smaller eddies required by the LES to support the momentum transport may result in an overestimation of the velocity gradient in the interface region, and an underestimation of the wall stress [43].

Given the complexities of the flow in the inner layer and its importance in establishing the production cycle [44], one cannot expect that modelling its effects on the outer flow may be easy or painless, or that equal accuracy can be achieved in wall-resolved or wall-modelled LES. It is, nonetheless, critical, over the next 15 years, to develop better wall-layer models and reduce the errors that can be expected from WMLES.

### (e) Subfilter-scale modelling

Modelling of the subfilter scales might not appear to be a pacing item. Several developments improve considerably on first-generation models [3–5,45]. Models that include more of the physical features of the smallest resolved scales have been shown to be accurate in a variety of flows, and a range of approaches is available that is likely to result in accurate prediction without increasing excessively the cost of the calculations [7–9,11,46]. However, several models have been proposed, in recent years, whose use significantly increases the computational cost of a calculation: they often require multiple filtering operations (which are difficult to perform in unstructured, multi-block or adaptive grids) while yielding only marginal improvements over existing models. Comparisons are often made on the same grid; a more consistent approach would be to compare for constant cost. Their benefits are unclear: in the author's experience, the first-order causes of discrepancies between LES and experiments are the grid resolution and the boundary conditions, as long as a reasonably accurate model (such as the ones cited above) is used.

Further development of SFS models appears to be more urgent in areas in which LES have not been as successful, those in which small-scale motions play a role. They include combustion and particle transport, for instance. Modelling the SFS effects on the flow, in these applications, will require different approaches than the ones used to model the momentum transport; for instance, PDF methods are more likely to be successful for the prediction of turbulent combustion [47]. This may be, in this author's opinion, one area in which modelling efforts should be focused in the near future.

As far as momentum transfer is concerned, however, there are two issues that need to be addressed: the completeness of the SFS model and its grid-dependence. The issue of completeness, discussed at length by Pope [16] stems from the fact that the model length scale in LES (the filter width Δ) is generally linked to the grid, rather than to turbulence quantities. The simulation then becomes user- and grid-dependent. Pope [16] proposes an ‘adaptive LES’ in which a resolution parameter (the percentage of unresolved energy, or dissipation, for instance) is the only specified quantity, and the grid (and filter width) are adapted to achieve this. Development of more physically consistent models, in which the length scale is related to the flow physics rather than to the grid, would be very beneficial in this sense, as well as allowing the achievement of grid-independent LES: presently, as the mesh is refined, statistics that depend on the small scales (even the Reynolds stresses) cannot be expected to reach grid convergence until the DNS limit is reached. If the filter width does not depend on the grid size, on the other hand, the limit of the LES if the grid is refined is a grid-converged LES, a conceptually more satisfying situation.

Figure 4 shows the dissipation length, (where is the turbulent kinetic energy and *ε* the dissipation rate of ) in two flows recently studied by the author's research group [48,49]. Even at the moderate Reynolds numbers of the simulations (both *o*(10^{4})), the integral scale changes by two orders of magnitude. If the filter width is taken as a fraction of the dissipation length scale, the percentage of energy carried by the resolved scales would be nearly constant. The grid could be then be coarsened in the regions where the integral scale is larger (as indicated in the figure) without compromising the accuracy of the simulation. In the wall jet case, for instance, the grid spacings in the far field can be 10–15 times larger (in all directions) than those near the jet. Current practice, when structured grids are used, only coarsens the mesh in one (at most two) directions.

## 4. Outlook

What will large eddy simulations look like in 2030? It is safe to assume that, between the increases in computational power that are still achievable, and the algorithmic developments, LES will be used in more complex configurations, and will yield reliable results for engineering analysis, if not design. The most significant advance might be achieved through the development of grid-optimized LES, which will be better able to tackle complex geometries. It will be characterized by a filter width that is some fraction of the integral scale, i.e. is defined by the physics of turbulence. A constant ratio between Δ and the integral scale will ensure that the LES resolve roughly the same percentage of the turbulent kinetic energy (and hence the momentum transport) everywhere. The ratio between filter width and grid size, Δ/*h* should be large to minimize numerical errors, and, ideally, remain nearly constant throughout the domain (although grid-generation constraints might make this difficult to achieve). The wall layer will, hopefully, be modelled accurately, so that the grid resolution is affordable, even at high Reynolds numbers.

Further model development is required to tackle problems in which small-scale motions are important (combustion, for instance). Judging from the past 15 years, progress in this area might still lag behind other areas of application. The computational cost of an SFS model must always be considered by the developer: some models proposed recently give marginal improvements in accuracy over calculations performed on the same grid with standard models, while increasing the cost substantially; they are also often tested only in very simple configurations, and their robustness is not established.

As Reynolds numbers closer to those in real-life applications can be reached, the information gleaned from LES can also be used to develop improved models for engineering applications. This has been always a stated objective of LES and DNS, but its achievement has been, so far, elusive. Most applications of LES to RANS models, for instance have compared model predictions (for each term of interest) to the results of the LES or DNS. This comparison is not always helpful, if the Reynolds numbers at which LES and DNS are performed are much lower than those ideally suited to turbulence models; furthermore, the solution of transport equations in the models includes dynamics that are absent from the kinematic comparisons most often seen in the literature. A better approach is that proposed by Parneix *et al.* [50], who suggest solving modified transport equations for the turbulence mode, in which one term is replaced by the ‘exact’ one obtained from the LES or DNS. Different paradigms may also be more successful (development of look-up tables for the model constants, for instance) but have not yet been explored.

In 2030, the full aircraft will still be beyond reach for pure LES; everything considered, that might not be the most useful application of the technique. Smaller problems, in which the turbulence physics are very complex (multiple scales, strong departures from equilibrium, massive separation, etc.) may remain the natural target for LES. The information generated by this method, with its ability to capture the most important turbulent eddies and the flow unsteadiness, can result in more complete understanding of the flow physics, more accurate predictions, and, eventually, the development of better engineering designs.

Niels Bohr used to say that predictions are very difficult, especially about the future. The only excuse for this author's attempt to make any is the fact that, by 2030, he will probably be retired, and the scorn of his colleagues for his lack of prophetic skills will not affect him any longer.

## Funding statements

This paper is financially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under the Discovery Grant program.

## Acknowledgements

The Canada Research Chair program is gratefully acknowledged.

## Footnotes

One contribution of 13 to a Theme Issue ‘Aerodynamics, computers and the environment’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.