## Abstract

The engineering tools of choice for the computation of practical engineering flows have begun to migrate from those based on the traditional Reynolds-averaged Navier–Stokes approach to methodologies capable, in theory if not in practice, of accurately predicting some instantaneous scales of motion in the flow. The migration has largely been driven by both the success of Reynolds-averaged methods over a wide variety of flows and the inherent limitations of the method itself. Practitioners, emboldened by their ability to predict a wide variety of statistically steady equilibrium turbulent flows, have now turned their attention to flow control and non-equilibrium flows, i.e. separation control. This review gives some current priorities in traditional Reynolds-averaged modelling research as well as some methodologies being applied to a new class of turbulent flow control problem.

## 1. Introduction

In studying aerodynamic flows, new technologies in flow control, renewed interest in high-speed (sub-orbital) vehicles, drag and/or noise reduction and improved propulsion systems are now the driving research incentives. Although the usual Reynolds-averaged methodologies have well-known deficiencies, they do remain the front-line methodologies for industrial prediction. Aerodynamic flows are characterized by a wide parameter range associated with, for example, Mach numbers and Reynolds numbers are turbulent, or are transitioning to or from turbulence. Of course, the challenge from a practical engineering standpoint is to be able to predict these flows, so that design cycle times are reduced and designs are optimized. It is clear that improving operating efficiencies are crucial in the development of next generation aircraft and aerospace vehicles. For this reason, the accurate and efficient prediction of turbulent flows has become an important topic of research.

Even with the current level of computational power, the direct numerical simulation (DNS) of such complex turbulent flows is not feasible and will not be so in the foreseeable future. The problem is simply the inability to resolve all the component scales within the turbulent flow. Even with computational grids exceeding 10^{9} grid points, simple fully developed channel flow simulations are still limited to the Reynolds number (based on friction velocity) of *O*(10^{3}). Therefore, either all or some of those scales that cannot be directly computed need to be modelled. The general task then is to develop the means to partition the resolved and unresolved scales, and then to develop suitable models for them. Within the context of scale partitioning, two (related) approaches exist.

The first approach is offered by the partitioning of the flow field into a mean and a fluctuating part, an idea first proposed by Osborne Reynolds (1895) over 100 years ago. This process, known as the Reynolds decomposition, assumes that the instantaneous flow can be partitioned into a fluctuating part, representing all the turbulent motion centred about a statistical mean value. This partitioning and accompanying ensemble averaging leads to a set of Reynolds-averaged Navier–Stokes (RANS) equations. Although this process eliminates the need to resolve the turbulent motion completely, its drawback is that unknown single-point, high-order correlations appear in both the mean and the turbulent transport equations. The need to model these high-order correlations is the well-known closure problem. The RANS method is a robust, easy to use and cost-effective means of computing both the mean flow and the turbulent stresses (velocity second moments) and has been, overall, a useful flow prediction technology.

The second approach is a filtering approach and is most commonly known as the large eddy simulation (LES) method. Here, the flow field variables are partitioned into resolved and unresolved scales. The original ideas were proposed over 40 years ago by Smagorinsky (1963) and extended a few years later by Deardorff (1970). In these early works, the formal analogy between the unclosed subgrid-scale stresses in LES and the unclosed Reynolds stresses in the RANS approach was exploited. The LES method is becoming a popular method for flow field predictions due to the rapid increase in computational power. While not as computationally demanding as DNS, proper implementation of LES—especially in the vicinity of solid boundaries—does require extensive computational resources. In addition, subtle issues associated with filter size and subgrid-scale stress models are still being debated. What remains to be determined is whether LES can be used with sufficient confidence to provide the accuracy now being demanded of the more established RANS methods. It is clear, however, that as computational power increases, the LES approach as well as variants of it (e.g. Geurts 2001) will become more prevalent.

The increasing demand for more detailed information about the flow fields over complex aerodynamic configurations or in turbine engines, and the need to control such flows, has led to the need for methodologies capable of capturing some part of the instantaneous motion. These turbulent flows may still be statistically steady or stationary, but now only a portion of the turbulent motion needs to be modelled. For statistically unsteady flows, the need for alternatives to the RANS approach is more obvious since RANS models are based on (spectral) equilibrium hypotheses that may no longer be valid (e.g. Carpy & Manceau 2006). The need for such predictions has led to the development of hybrid methods capable of having a RANS-type behaviour in the vicinity of a solid boundary and an LES-type behaviour away from the boundary. Such methods can also be interpreted within the usual LES context but with specially designed subgrid-scale (SGS) models capable of the dual behaviour just described. The most often used approach to date is the detached eddy simulation (DES) approach developed by Spalart and colleagues (Spalart *et al*. 1997; Spalart 2000). An insightful discussion of this topic, from an LES perspective, is given by Piomelli & Balaras (2002). Hybrid methods can be essentially classified into two categories: methods where the domain is decomposed into zones where LES is performed and zones where RANS is performed, with a sharp interface or, possibly, a small overlap region; or methods where the governing set of equations is smoothly transitioning from a RANS behaviour to an LES behaviour, based on criteria updated during the computation. The models in the first category are often termed as *zonal* and the models in the second category as *non-zonal* (Bunge 2006), although this terminology is ambiguous since both use different models in different zones. Both kinds of models will be a topic of interest and active research in the future. Currently, while such methods are often employed, little research has been conducted into formulating a consistent mathematical theory where such blending of methodologies is employed. Although the presentation will be somewhat general at the outset to show how the various approaches, at least formally, have the same form-invariant equations, the distinction being in the closure of these equations, the focus of this review will be on the traditional RANS approaches. Indeed, even in the hybrid RANS/LES frame, RANS modelling remains a cornerstone: the hybrid models strongly rely on a RANS model in the near-wall region, such that the influence of pressure gradients on the boundary layer and the prediction of flow separation are completely determined by the performance of the RANS model. Moreover, using a RANS-type model as SGS model can be useful to account for complex phenomena such as heat transfer, effects of rotation or stratification, etc. In industrial applications, where a cost reduction is needed, the use of coarse meshes is incompatible with the hypotheses of simple inertial behaviour of the subgrid scales, and the use of RANS-like transport equations and/or complex constitutive relations is and will be an active research field (e.g. Chaouat & Schiestel 2005).

From a physical standpoint in RANS model development, the task is to characterize the turbulence. One obvious characterization is to describe correctly the evolution of representative turbulent velocity and length-scales, an idea that originated over 60 years ago (Kolmogorov 1942). This is the physical motivation behind the development of turbulent closure models. There is a hierarchy of turbulence models currently available for solving aerodynamic flow problems. These classes include the differential Reynolds stress models (DRSM), the algebraic Reynolds stress models (ARSM), the nonlinear eddy–viscosity models (NLEVM) and the linear eddy–viscosity models (LEVM). While other alternatives may exist, these form the basis for the majority of methodologies used, and of these, the linear eddy–viscosity models are the most common.

The most commonly used linear eddy–viscosity models in aeronautics are the Spalart and Allmaras (SA) model (Spalart & Allmaras 1994) and the shear stress transport (SST) model (Menter 1994), but in general purpose codes, different forms of high-Reynolds number *K*–*ϵ* models, such as the standard *K*–*ϵ* model (Launder & Spalding 1974), the RNG *K*–*ϵ* model (Yakhot & Orszag 1986) and the *K*–*ϵ* model of Shih *et al*. (1995*a*), often termed the ‘Realizable *K*–*ϵ*’ model, are still widely used. The main improvement of the RNG and Shih *et al*. *K*–*ϵ* model when compared with the standard *K*–*ϵ* model is the correction of the so-called *stagnation point anomaly*, through a modification of the *ϵ*-equation (RNG model) or the eddy viscosity (Shih *et al*.). The SA model is a one-equation model based primarily on empiricism and on dimensional analysis arguments. It is easily used with any type of grid: structured or unstructured, single block or multiple blocks. The SA model has become popular among industrial users due to its ease of implementation and relatively low cost. The SST model has gained increasing favour due primarily to its robust formulation and improved performance for separated flows. It is a blend of the original *K*–*ω* formulation near walls and a *K*–*ϵ* formulation in the outer region and in the free shear flows. An important feature of the SST model is the modification to the definition of eddy viscosity to account for the effect of the transport of the principal turbulent shear stress. More detail about the mathematical form of these models can be found in the review by Gatski & Rumsey (2002).

In aerodynamic flows where the turbulence is not confined to relatively narrow regions of the flow domain, such as on multi-element configurations at a high angle of attack, wing–body junctions or turbine blades where separation or near-separation conditions exists, the demands for both accurate computational modelling of the physical problem as well as the quality of the turbulence model increase significantly. The role of the CFD practitioner becomes proportionately more critical since poor modelling of the physical problem and/or poor choice of turbulence model can lead to correspondingly poor predictions.

These increased demands on the turbulence model performance have also highlighted inherent deficiencies in the models themselves. For example, their inability to consistently predict separation accurately or to account properly for pressure–gradient effects became apparent in trying to replicate multi-element aerofoil flow fields at higher angles of attack (Rumsey *et al*. 1998; Rumsey & Gatski 2001). Other examples are the difficulties in accounting correctly for rotational effects (Jakirlić *et al*. 2002) and wall heat transfer (Thielen *et al*. 2005), which are both relevant to turbomachinery. A new deficiency has also been identified that is not directly related to turbulence model prediction itself, but to the application of turbulence models in predicting aerodynamic flow fields. This deficiency is related to the turbulence model being properly sensitized to transition onset location. While turbulence models are calibrated for fully turbulent flow predictions, aerodynamic flows include regions where the flow may be laminar and transitioning to fully turbulent. In attempting to predict and replicate such flows, a single system of equations is used throughout the computational domain and these equations are expected to predict the entire flow field, which can include regions of laminar flow as well as fully turbulent flow. Some of these issues will be discussed further here.

## 2. The filtering process and governing equations

As detailed in §1, it is currently necessary for numerical calculation of practical engineering turbulent flow fields to solve a set of equations for flow variables that represent the motion of a limited spectral range of scales. This description holds true for LES, RANS formulations and any of the newly developed hybrid or composite methodologies currently being proposed. As such, the equations describing the filtered motions in any of these formulations are form invariant. They obviously differ with respect to the flow field motion being described, and this is predicated on how the higher-order correlations are parameterized.

As is customary in all these formations, the flow variable *f* is decomposed into a filtered part, , and a sub-filtered part, *f*′, given as(2.1)Generally, the filtering process can be defined as a subset of the general operation (e.g. Sagaut 2006)(2.2)Different forms for the convolution kernel can be associated with the various solution methodologies. For the Reynolds-averaged formulation, for example, stationarity is usually assumed and a long-time average then corresponds to an ensemble average. The filter function in equation (2.2) is given by(2.3)such that(2.4)In such flows, the entire spectral range of scales is modelled, so the sub-filtered part *f*′ is a fluctuating quantity whose average is zero, , and the mean quantity can be extracted from(2.5)Note that, since , Reynolds averaging cannot be expressed as a convolution filter, although the Reynolds-averaged function is the limit of a series of functions obtained by convolution, as shown by equation (2.5). As a practical matter, Reynolds averaging is evaluated as long-time averaging, where ‘long time’ is not infinity but rather a time that is sufficiently large when compared with the turbulence time scale. Therefore, Reynolds averaging can be considered, within any predefined accuracy, as the convolution filter with a sufficiently large . A useful property of this filtering is that the average of the product of two quantities is . This allows for the easy extraction of stationary correlation data from numerical simulations where the instantaneous values of *f* and *g* are computed.

For a flow statistically periodic in time (cyclostationary), it is customary to use phase averaging, corresponding to a filter function given by(2.6)(2.7)where *T* is the period of the cycle. In that case, the phase average of the sub-filtered part *f*′ is zero, , and the filtered, or phase-averaged quantity 〈*f*〉 can be extracted from(2.8)

For the most part, the large eddy simulation methodology has been based on spatial filtering. The convolution filters most commonly employed are the top-hat filter, the Gaussian filter and the spectral cut-off filter. In numerical algorithms requiring structured meshes applicable to multi-dimensional flows, such spatial filtering can become problematic since anisotropies in the filtering process need to be introduced. For unstructured grid schemes, the anisotropy issue is compounded. An alternative approach that helps to alleviate this problem is the Eulerian temporal filtering approach. (It is recognized that the problem cannot be completely solved since an inherent filtering occurs even on the spatial scales with a temporal filtering process.) This procedure has received little attention in the past, although some recent works (Pruett 2000; Pruett *et al*. 2003) suggest that it can be an effective alternative to the spatial filtering process. Causal time-domain filters can be constructed that are analogous to the spatial filters. A simple example of a causal filter is the top-hat filter given by equation (2.3), where is the temporal filter width. Analogous to the spatial filtering approach, only a portion of the spectral range of scales is modelled with the sub-filtered part *f*′. The filtered quantity can then be extracted from equation (2.4). It can thus be seen that, within the realm of temporal filtering, it is possible to develop a more rigorous linkage between the large eddy and Reynolds-averaged approaches (Pruett *et al*. 2003), since the Reynolds-averaged function *E*{*f*} is the limit of the temporally filtered function when the temporal filter width goes to infinity. In the frame of spatial filtering, such a linkage can only be established in homogeneous flows.

The aerodynamic flows of current relevance are compressible and the relevant governing equations, for whichever solution methodology is used, should be cast in the appropriate form. This is becoming even more important as the interest in high-speed flow fields increases. It should be noted, however, that having a compressible turbulent flow does not mean that the turbulence dynamics is necessarily compressible. For attached flows, mean or filtered flow fields over most of the supersonic regime (*M*<5) are described reasonably well using incompressible dynamics under adiabatic conditions and in the absence of shocks.

The filtering processes described above yield governing equations that are all formally equivalent. They differ, not in form, but in the models or parameterizations needed to close the equations. The formal similarity of these equations, irrespective of the filtering process invoked to derive them, is the foundation upon which hybrid or composite methods can be constructed. However, the modelling or parameterizations of the various terms that require closure is a difficult and challenging problem since the scales of motion being represented by each filtering procedure can be different. In addition, for compressible flows, it has been found that a rewriting of the equations using mass-weighted (Reynolds 1895) or Favre variables (Favre 1965) is advantageous, since the equations take a more compact form and are structurally similar to their incompressible counterparts.

For a dependent variable *f*, the Favre mean is defined as(2.9)The instantaneous value *f* can then be decomposed into either the usual Reynolds-averaged variables or the Favre-averaged variables(2.10)As might be expected, an extensive list of relations exists between the Reynolds-averaged variables and the Favre variables (Barre *et al*. 2001).

The equations for the mean density and mean momentum are given by(2.11)(2.12)and the mean viscous stress tensor is(2.13)where is the mean molecular viscosity and is the Favre-averaged velocity correlation tensor. Equation (2.13) neglects contributions from *μ*′ and assumes that . This assumed equality between the average velocities implies that the average fluctuating velocity is small since .

The equation for the mean total energy is(2.14)where(2.15)(2.16)(2.17)and(2.18)In the approximation of equation (2.17), fluctuations in the thermal conductivity are neglected, and the Favre- and Reynolds-averaged mean temperatures are taken as approximately equal. The equation of state in mean variables is(2.19)or in terms of the mean total energy (2.20)where *γ* is the ratio of specific heats (*c*_{p}/*c*_{v}) and *R* is the gas constant. The presence of the turbulent kinetic energy term *K*,(2.21)suggests a strong coupling between the mean equations and the turbulent transport equations.

## 3. Closure strategies for RANS

The closure problem arises at all levels of the statistical moment equations. At the mean flow level, the RANS momentum equations, equation (2.12), require closure through a specification of _{ij}, and the mean energy equation, equation (2.14), requires the closure of the turbulent heat flux . At the second-moment level, closure is required for the velocity–pressure gradient correlation, the tensor dissipation rate, the velocity-triple moments and the turbulent mass flux . Closure for these terms is achieved either through the solution of partial differential transport equations or by assuming a tensor polynomial expansion for the unknown correlation in terms of basis tensors formed from the independent tensors on which it depends (Gatski 2004). It is not possible in the space allotted here to discuss the modelling of all these unknown higher-order correlations. Only the correlations related to the velocity–pressure gradient correlation and the tensor dissipation rate will be discussed. These correlations are relevant to the solution of both incompressible and compressible flows and, in the case of the tensor dissipation rate, it has been found that the solenoidal (incompressible) dissipation rate need only be accounted for (e.g. Kreuzinger *et al*. in press). It is safe to assume that the newer hybrid methodologies have not reached a level of maturity where detailed modelling issues associated with many of these correlations have been addressed.

### (a) Levels of turbulence closure

It is often desirable to work with the turbulent anisotropy tensors or scaled scalar flux vectors in analysing turbulent flows. For the Reynolds stress tensor, the anisotropy tensor , where , can be defined, and for the heat flux vector, for example, the scaled scalar flux vector can be defined.

The differential Reynolds stress model requires the solution of six partial differential equations for the Reynolds stress components and a transport equation for a variable from which a length-scale can be obtained. For the most part, this scale-related equation is the turbulent kinetic energy dissipation rate equation. In addition, various scalar flux equations, each having three components, are also required. Associated in general with the scalar flux equation is the need for a scalar variance and corresponding scale equation for the scalar variance. In general, the full differential form at the second-moment level can be computationally challenging.

Subsets of this full differential set are often desired and can be obtained. Implicit algebraic Reynolds stress and scalar flux equations can be obtained directly from the differential forms by invoking what are termed weak equilibrium assumptions. Applicable models are developed through explicit polynomial expansions (representations) of the anisotropy tensor and scaled scalar flux vector. The proper choice of basis is dependent on the functional dependencies associated with the anisotropy tensor and scalar flux vector. Note that the equations for the Reynolds stress anisotropy *b*_{ij} are only frame-indifferent under the Galilean group of transformations. Under the more general Euclidean group, the equations exhibit a dependence on system rotations and translational accelerations. Transport equations are still required since the algebraic models developed simply yield an expression for the components of the Reynolds stress anisotropy tensor and the scalar (heat) flux vector. These transport equations are associated with the turbulent kinetic energy and the energy dissipation rate (or dissipation rate per kinetic energy), and scalar variance and scalar variance length-scale.

What are termed nonlinear eddy-viscosity models form a link between what one would call the higher-order models—the differential and algebraic Reynolds stress models and the lower-order linear eddy-viscosity models. The nonlinear eddy-viscosity models describe the turbulent Reynolds stress field by a polynomial expansion that is similar to the algebraic stress model; however, the expansion coefficients are determined from calibrations with experimental or numerical data, and on some physical consistency constraints (e.g. Shih *et al*. 1995*b*; Craft *et al*. 1996). This contrasts with the algebraic stress models where the expansion coefficients are extracted directly from the Reynolds stress transport equations. An analogous situation holds for the scaled scalar flux vector. Once again, transport equations are required as in the case of the algebraic anisotropy tensor and algebraic scalar flux vector representations.

At the linear eddy viscosity and diffusivity level of closure, the mean momentum and temperature equations are closed by using a Boussinesq-type approximation between the turbulent Reynolds stress and the mean strain rate tensors, and scaled scalar flux and gradient of the mean scalar property. Although computational resources have become less of a prohibitive factor in choosing turbulence models, many practitioners have continued to rely on the lower-order eddy-viscosity and diffusivity models. Much of this is due to robustness issues associated with the higher-order models as well as to lack of familiarity with the formulations.

### (b) Improved high-order correlation closures

The main focus of advanced RANS turbulence modelling research for aerodynamic flows is on correlations associated with the turbulent velocity field in isothermal flows. This does not suggest that modelling issues associated with the scalar fluxes (e.g. heat and mass) are unimportant, rather it is a reflection on the limited scope of advanced modelling activity currently being supported. For the turbulent velocity field, the differential Reynolds stress transport equations are the starting point for model development. The two unknown correlations that currently receive the most attention are the velocity–pressure gradient and tensor dissipation rate correlations. These terms are important because they represent the mechanisms responsible for the redistribution of normal and shear stress throughout the entire boundary layer flow including the near-wall region.

#### (i) Pressure–strain rate correlation

In the case of the velocity–pressure gradient correlation, it is customarily rewritten in terms of the pressure–strain rate correlation and the pressure–velocity correlation as(3.1)

Initially, this partitioning was intended to isolate the effects of the pressure transport terms, which predominantly contribute to the dynamics near solid boundaries. The pressure–strain rate correlation term *Π*_{ij} is then decomposed into rapid and slow parts. The rapid part is modelled by assuming dynamic equilibrium conditions in both physical and spectral space, and the slow part is modelled by assuming either a linear or nonlinear relationship with the Reynolds stress anisotropy tensor. The wall proximity corrections required will be discussed separately below. As might be expected, an extensive literature base deals with the full range of issues associated with developing specific models (Hanjalić & Jakirlić 2002). Almost all the research efforts have focused on the pressure–strain rate correlation since this term is of the same order as the production term and acts as a redistribution between the Reynolds stress components. In compressible flows without shocks and close to adiabatic conditions, the partitioning in equation (3.1) is applicable and for the most part the pressure–velocity term is neglected. In compressible flows with shocks, the pressure–velocity correlation has a non-negligible effect on the dynamics even away from solid boundaries and should be accounted for.

Away from solid boundaries, the pressure–strain rate correlation also acts to diminish the difference between the normal stress components. In the vicinity of the solid boundaries, its action, along with pressure–velocity correlation, increases the anisotropy of the stress field. This action, along with that of the tensor dissipation rate, is to enforce the two-component limit on the Reynolds stress tensor. In aerodynamic flows of interest, these inhomogeneous effects can be important, so it is clear why such terms have received so much attention.

In all of the commonly used second-moment closure models, the pressure–strain rate correlation *Π*_{ij} is modelled away from solid boundaries in the general form (Chou 1945; Lumley 1978) as(3.2)where _{ij}(** b**) and

_{ijkl}(

**) are tensor functions of the anisotropy tensor and are related to integrals over the flow volume derived from a pressure Poisson equation for incompressible flows or a convective pressure wave equation for compressible flows, and**

*b**ϵ*is the isotropic dissipation rate. By the second law of thermodynamics the isotropic dissipation rate is always positive, and this Reynolds-averaged turbulent dissipation rate can be written as(3.3)where , { } is the trace and [ : ] is the trace of the product. The second term on the right is associated with fluctuating dilatation and can be neglected (Pantano & Sarkar 2002; Pirozzoli

*et al*. 2004), and the first term can be related to the fluctuating enstrophy (e.g. Sarkar

*et al*. 1991).

In equation (3.2), the _{ij}(** b**) term is usually associated with the ‘slow’ relaxation of the turbulence towards isotropy and the

_{ijkl}(

**) term is usually associated with the ‘rapid’ response of the turbulence to imposed mean velocity gradients. This partitioning has its origins in incompressible flows where the turbulent pressure field**

*b**p*′ is split into slow

*p*′

^{(S)}and rapid

*p*′

^{(R)}parts. The

*p*′

^{(S)}part is the solution of a Poisson equation that only involves gradients of the turbulent velocity field; the

*p*′

^{(R)}part is the solution of a Poisson equation involving the mean velocity gradients. The terminology originated from the observation that the slow part will only adjust as the turbulence itself adjusts, but the rapid part will adjust instantly through the mean velocity gradient.

For statistically stationary turbulence, the starting point for the model development of the rapid part can be treated with the transform pair(3.4a)and its equivalent Fourier transform,(3.4b)where **r** is the two-point separation distance (** r**→0 for single-point closures), and are the transformed fluctuating velocity and pressure fields. From this equation, it is now possible to obtain an expression for the rapid part of the pressure–strain rate correlation once an appropriate representation for the energy spectrum tensor has been found. (The slow part of the pressure–strain rate correlation, representing turbulence–turbulence interactions, yields triple-velocity terms. Models for this slow contribution are derived in an alternative manner by simply relating it to the Reynolds stress anisotropy.)

For incompressible flows, the fourth-order tensor _{ijkl}(** b**) given in equation (3.2) is related to the energy spectrum tensor through the simple relation(3.5)where the energy spectrum tensor

*E*

_{ij}(

**,**

*b***) can be represented by a polynomial expansion involving both the Reynolds stress anisotropy and the wavenumber vector. The simplest representation that can be used is given by a four-term representation(3.6)where**

*k**E*(

*k*) is the isotropic energy spectral density and

*E*

_{a}(

*k*) is the anisotropic energy spectral density.

With the functional dependency outlined in equation (3.2), the proper representation for ** Π** would be generated from the integrity basis given by the invariant combinations of

**,**

*b***and**

*S***W**(e.g. Gatski 2004). Obviously, any representation composed of the full integrity basis would be unmanageably large. Even though the energy spectrum tensor given in equation (3.6) is only a linear function of the anisotropy tensor

**, many representations for the rapid part of the pressure–strain rate correlation have used bases with quadratic and cubic terms in**

*b***(cf. Sjögren & Johansson 2000) although they retain only a linear dependence on the mean velocity-gradient field. This linear dependence on the mean velocity gradient is consistent with the functional form for the pressure–strain rate correlation given in equation (3.2); however, the appearance of the higher-order powers of the anisotropy tensor is inconsistent with the linear dependence on the anisotropy tensor assumed for the energy spectrum tensor.**

*b*#### (ii) Dissipation rate tensor

The other unknown correlation that appears in the differential Reynolds stress equations and that continues to be of modelling interest is the tensor dissipation rate *ϵ*_{ij}. Even though it is possible to derive a transport equation for from the fluctuating momentum equation by taking the moment (cf. Speziale 1991)(3.7)where is the Navier–Stokes differential operator, the final form of this equation is quite complex and contains terms where very little, if any, information is available for closure modelling.

It has been assumed that, away from solid boundaries, the tensor dissipation rate becomes isotropic (*ϵ*_{ij}=2/3*ϵδ*_{ij}), although there have been suggestions that such an assumption may not generally be true (Durbin & Speziale 1991). Thus, early attempts at modelling the deviatoric part of the dissipation rate tensor focused on a coupling with the Reynolds stress anisotropy (Hanjalić & Launder 1976; Hallbäck *et al*. 1990).

The quantity represents more than a transformation of kinetic energy to internal energy since it can be further decomposed into(3.8)The second term, , is the divergence of a flux, and as such is purely a transport term since it redistributes energy from one region to the other, but does not contribute to the global time evolution of kinetic energy. This term has been evaluated by Bradshaw & Perot (1993) in a channel flow and represents less than 2% of the total.

The energy dissipation rate also plays a pivotal role in the dynamics of the turbulent kinetic energy. It can be extracted from the viscous term,(3.9)which is not necessarily positive and requires modelling, that appears in the transport equation for the turbulent kinetic energy. It is customary to decompose this term into(3.10)where the first term, , is called *turbulent diffusion*. This decomposition has two very desirable properties: does not require modelling and is always positive which facilitates its modelling.

The variable *ϵ* is often called the *homogeneous* dissipation, since in homogeneous flows, and *ϵ*′=*ϵ*−*ϵ* is sometimes called the *inhomogeneous* dissipation. It is worth emphasizing that the decomposition of homogeneous/inhomogeneous dissipation is completely artificial, since there is an infinite number of possible decompositions: for instance, any linear combination of the form *ϵ*+*βD*^{μ} can be called *homogeneous* dissipation, since it goes to *ϵ* in homogeneous flows. In the differential Reynolds stress formulation, where the isotropic assumption for the dissipation rate tensor is used, this isotropic dissipation rate is obtained from a modelled transport equation (Hanjalić & Launder 1976).

There have been recent attempts to account for dissipation rate anisotropy through a tensor dissipation rate equation (e.g. Oberlack 1997; Speziale & Gatski 1997; Jakirlić & Hanjalić 2002). While not explicitly focused on developing closures to account for wall proximity effects, these studies used tensor representations to extract either an explicit algebraic dissipation rate model or a differential model for the tensor dissipation. Unfortunately, extensive validation studies are difficult due to the lack of data available for a quantity such as *ϵ*_{ij}. Numerical simulations provide the best source of data for such quantities but are unfortunately limited to simpler flow geometries.

It is now recognized that the effects of anisotropies in the turbulence statistics should be accounted for in some way for many flows of practical interest. Obviously, the full differential models do this, but at a higher computational overhead cost than the linear eddy-viscosity models. One of the main reasons that algebraic vector and tensor polynomial representations have become popular is that such anisotropies can be accounted for with only a small increase in computational cost over the linear eddy-viscosity models. The details of this representation procedure have been presented previously (e.g. Gatski 2004; Gatski & Wallin 2004) and it is not necessary to repeat them here.

### (c) Modelling of near-wall turbulence

While the discussion up to this point has been focused on mathematical representations of the rapid pressure–strain rate correlation in the case of homogeneous flows, the practical application of such models to inhomogeneous flows requires the effects of solid boundaries to be accounted for. This inherently brings into consideration the proper accounting of turbulence anisotropy into the modelling process. Although this has sometimes been called low Reynolds number modelling, it is more properly called near-wall modelling to distinguish it from the separate problem of transition prediction to be discussed in §4.

Both the pressure–strain rate correlation and the anisotropic dissipation rate that appear in the *τ*_{ij} transport equation require modification, and the pressure–velocity correlation can no longer be neglected. Along with these parameterizations, which focus on the stress transport equations, is the behaviour of the dissipation rate equation itself in the vicinity of the wall. Overall, the issue of near-wall modelling is at least as complex as the issue of developing high-Reynolds number models, and unfortunately less precise.

Several attempts have been made to develop near-wall closure corrections for the Reynolds stress transport equations and the dissipation rate equation. For the turbulent stress transport equations, attempts have almost exclusively focused on extending the high-Reynolds number pressure–strain models. This approach has been taken by Launder *et al*. (1975) who have developed a methodology that enforces the two-component limit of the turbulent Reynolds-stress field as the solid boundary is approached (Craft & Launder 2001). For the dissipation rate equation, additional terms and modified coefficients have been proposed to account for the anisotropic near-wall effects and the correct limiting behaviour at the wall (e.g. Speziale & Gatski 1997; Jakirlić & Hanjalić 2002).

A major focus over the last decade has been on the ‘elliptic relaxation method’ (Durbin 1991, 1993; Manceau & Hanjalić 2000; Manceau *et al*. 2001, 2002; Laurence *et al*. 2005). Based on a theoretical analysis of the influence of the blocking of the wall, this approach introduces a tensor function representing the combined effects of a near-wall velocity–pressure gradient correlation and anisotropic dissipation rate. The tensor function asymptotes to a high-Reynolds number form away from solid boundaries through an elliptic relaxation equation.

It is worth briefly outlining the formulation. The transport equation for *τ*_{ij} can be written in the form(3.11)where *ϕ*_{ij} is not traceless (*cf*. *Π*_{ij}), since it includes pressure diffusion (see equation (3.1)). In the elliptic relaxation method, the variation of the dissipation rate anisotropy (*d*_{ij}=(*ϵ*_{ij}−2*ϵδ*_{ij}/3)/(2*ϵ*)) as the wall is approached is accomplished by a relaxation to its wall value, which is assumed to be equal to *b*_{ij}. With this assumption, *ϕ*_{ij}−*ϵ*_{ij} in equation (3.11) can be written as(3.12)with the relaxation function *f*_{ij} defined by(3.13)

The original scaling of the relaxation function *f*_{ij} was solely through the turbulent kinetic energy; however, Manceau *et al*. (2002) have shown that adding a dissipation rate factor *ϵ* to the scaling (*ϵKf*_{ij}) eliminates an unwanted amplification effect that is inherent in the original scaling. The system is closed through the solution of a relaxation equation for *f*_{ij}(3.14)where *Π*_{ij} can in principle be any of the high-Reynolds forms found in the literature. Although some authors have tested nonlinear models (Wizman *et al*. 1996), terms linear in *b*_{ij} are generally used. Moreover, despite the fact that pressure diffusion can readily be accounted for in equation (3.14), this term is usually considered to be negligible far from the wall, i.e. in the right-hand side. Thus, the elliptic relaxation equation is driven by the high-Reynolds number form of the pressure–strain rate correlation ** Π** and a contribution from the Reynolds stress anisotropy 2

*ϵ*

**(away from the wall, the dissipation rate is assumed to be isotropic**

*b**d*

_{ij}=0).

It is interesting to contrast the predictive capabilities of previous near-wall corrections with the performance of the more rigorous elliptic relaxation approach. One of the first Reynolds stress models that was proposed to handle near-wall effects was developed by Hanjalić & Launder (1976). Figure 1 shows a comparison of results from the elliptic relaxation method (Manceau *et al*. 2002) and the Reynolds stress model (Hanjalić & Launder 1976) with direct simulation data at *Re*_{τ}=180 for the Reynolds shear stress (in wall units). The results of the Speziale, Sarkar and Gatski (SSG) model (Speziale *et al*. 1991), without near-wall correction, are also included for comparison. At this low Reynolds number, the RSM approach based on a modification of the tensor dissipation rate and the introduction of an extra production term was not able to accurately predict the near-wall behaviour of the Reynolds shear stress. The elliptic relaxation method, on the other hand, did an excellent job at this Reynolds number, and as shown in Manceau *et al*. (2002), it also made excellent predictions over a much larger Reynolds number range.

However, Reynolds-stress models based on elliptic relaxation have not spread into industrial codes and applications, owing to the huge additional numerical effort required when compared with simpler closures (the *ν*^{2}–*f* model for instance, which includes elliptic relaxation in a linear eddy–viscosity framework, has become popular and is now available in several commercial codes). The reason lies not only in the requirement to solve six additional differential equations, equation (3.14), but also in the numerically problematic wall boundary conditions for the components *f*_{22}, *f*_{12} and *f*_{23}(3.15)where *x*_{2} is the wall-normal direction. (The boundary conditions for the other components are *f*_{11}=*f*_{33}=−(1/2)*f*_{22} and *f*_{13}=0.) One of the difficulties is that these boundary conditions must be applied in the local frame determined by the orientation of the wall, which couples the Reynolds stress equations via the boundary conditions when the wall is not aligned with the global reference frame.

Manceau & Hanjalić (2002) have proposed the elliptic blending Reynolds stress model (EB–RSM), which tries to preserve all the desirable properties of the elliptic relaxation model, in particular, those properties relating to the representation of the blocking effect of the wall, while decreasing the number of equations and the numerical stiffness of the boundary conditions. In the most recent version of the model (Manceau 2005), *ϕ*_{ij}−*ϵ*_{ij} is modelled as(3.16)

The purpose of the introduction of the near-wall form is the reproduction of the asymptotic behaviour of *ϕ*_{ij}−*ϵ*_{ij}, a role that is played by the boundary conditions equation (3.15) in the elliptic relaxation model. Here, the near-wall form of the dissipation tensor is chosen as , but any other choice is possible: the important point being that to satisfy the behaviour of , the model for must be adapted to the model for . In that case, *ϕ*_{ij} must take the form(3.17)

Contrary to boundary conditions that can be imposed in a local frame linked to the wall where they are applied, equation (3.16) is applied *inside* the flow domain, and requires information about the wall distance and the orientation of the wall. The wall distance is felt implicitly by the blending function *α*, which goes from 0 at the wall to 1 far from the wall, since it is obtained via an elliptic relaxation equation, very similar to those used in the elliptic relaxation model(3.18)with the boundary condition *α*=0 at the wall. In a semi-infinite domain bounded by a flat plate, the analytical solution of this equation would be(3.19)if *L* were taken as a constant. Therefore, *α* is a function of the ratio *y*/*L*=*wall distance*/*integral length-scale* of the model. Figure 2 shows isocontours of *α*^{2} around the step corner in the case of a backstep flow at *Re*=37 500. It can be seen that the width of the region over which the near-wall model is active is strongly dependent on the local flow conditions, via the length-scale *L*, which is much larger in the recirculation region than in the incoming boundary layer.

The second important issue is the orientation of the wall. It is crucial that the model is valid regardless of the relative orientations of the wall and the reference frame, i.e. the near-wall model must be objective. Indeed, since the *ϕ*_{ij} tensor is objective, it is crucial to ensure that is as well. This requirement is automatically ensured if is made a function of objective quantities. In order to enforce the asymptotic behaviour equation (3.17) correctly, it is seen that the model cannot be proportional to *τ*_{ij}, but must rather be able to distinguish between the components, depending on the orientation of the wall. Therefore, an *objective* tensor must be defined that provides information about the orientation of the wall. First, it can be noted that the vector ** n**=∇

*α*is normal to the wall in its vicinity, since the wall is the isocontour

*α*=0, and is objective, since it is the gradient of an objective scalar. Thus, the tensor(3.20)is also objective, and a simple way to reproduce equation (3.17) is to use the model(3.21)This near-wall form is then transitioned to the high-Reynolds number form(3.22)Here, the SSG model (Speziale

*et al*. 1991) is used for

**, although any high-Reynolds number model could be used.**

*Π*Thus, with such a blending approach, the distance to the wall and its orientation are felt by the model through *α* and the tensor N. Since no geometrical information is explicitly introduced in the equations, the model is easy to use in complex geometries. Moreover, the number of equations necessary to enforce correct near-wall behaviours is reduced from six to one, when compared with the elliptic relaxation model, and the wall boundary conditions are much simplified, which is crucial from a numerical stability point of view. It is worth pointing out that the model is able to reproduce the near-wall anisotropy of turbulence as well as or better than the full elliptic relaxation model, as shown in figure 3.

Another focus of near-wall research, that has been and continues to be of interest, is the modelling of the near-wall anisotropy *d*_{ij} of the dissipation tensor *ϵ*_{ij}. Contrary to the elliptic relaxation/blending strategies, discussed above, in which the difference *ϕ*_{ij}−*ϵ*_{ij} is modelled as a whole, many authors have investigated the dissipation tensor separately. Provided that all the models are based on an isotropic dissipation far away from the solid boundaries, the modelling challenge consists of providing a smooth transition towards an anisotropic near-wall formulation.

The decomposition equation (3.10) between viscous diffusion and homogeneous dissipation discussed previously is not without modelling implications for the dissipation tensor. Indeed, all the models are based on an algebraic formulation of the dissipation tensor, which ensures a smooth transition between a high-Reynolds form, applied in regions far away from solid boundaries, and a near-wall form, derived in order to reproduce the asymptotic behaviour of the dissipation tensor in the vicinity of a wall. In the standard decomposition, where the homogeneous dissipation tensor *ϵ*_{ij} is defined by(3.23)the asymptotic behaviour of the anisotropy of the dissipation *d*_{ij} is difficult to relate to the anisotropy of the Reynolds stress *b*_{ij}. Indeed, for a wall located in *x*_{2}=0, the components of *d*_{ij} and *b*_{ij} are asymptotically related by(3.24)which shows that the two tensors are not proportional. Therefore, the standard way of modelling the dissipation anisotropy(3.25)where *f*_{ϵ} is any function approaching unity at the wall and going to zero far away from the wall, does not provide the correct limiting anisotropy for components 12, 33 and 22 in figure 4.

The analysis of the transport equations of the two-point correlation tensor (Jovanović *et al*. 1995) suggests that the decomposition of viscous diffusion/homogeneous dissipation should be written as(3.26)Jakirlić & Hanjalić (2002) have pointed out that with this decomposition it is much more justified to assume proportionality between the anisotropies(3.27)since the asymptotic analysis leads to(3.28)although component 22 is still not exact. Figure 4 shows a comparison of the model *f*_{ϵ}*b*_{ij} with the dissipation anisotropy obtained from the DNS of Moser *et al*. (1999), based on either the standard decomposition equation (3.23) or the decomposition equation (3.26) proposed by Jakirlić & Hanjalić (2002). The blending function is modelled as , where *A* and *E* are the flatness parameters associated with the Reynolds stress and the dissipation, respectively (Jakirlić & Hanjalić 2002). It can be seen that the model represents better than *d*_{ij}, although, quite surprisingly, the improvement is especially visible on components 11 and 33.

These types of approaches, focusing on either the terms involving the pressure fluctuations or the tensor dissipation behaviour, have yielded encouraging results in an effort to account more accurately for near-wall behaviour and improve the robustness of such numerical predictions. As with the other closures discussed in this paper, only extensive validation tests will confirm their ability to handle a wide range of flow fields.

## 4. Some current challenges

The challenging problems to be briefly discussed here are by no means complete. Arguably, other application areas are of current high interest, including buffet onset and compressible mixing. However, the two rather broad areas discussed here have, in our view, been pacing items in the migration from traditional RANS methodologies to hybrid RANS–LES methods and ultimately to large-scale simulation methodologies.

### (a) Accounting for transition and relaminarization

Many aerodynamic flow fields of interest, where single-point closures based on the Reynolds-averaged Navier–Stokes methodology are applicable, describe regions where the flow is transitioning from laminar to turbulent or relaminarizing from the turbulent state. RANS models have routinely been applied to such regions and the accuracy of the flow predictions in these ‘non-turbulent’ regions have varied significantly depending on the models. Obviously, the RANS turbulence models were never calibrated to predict such transitioning regions; however, the current demands on RANS predictions of aerodynamic flows requires that such models predict the flow in these transitioning regions with some degree of accuracy.

At the outset, it is necessary to recognize that what is required of a RANS model is not the ability to predict the detailed transition dynamics *per se,* but rather that the RANS model be properly sensitized to predict the statistical characteristics of the flow in these transitioning regions. The basis for many of the models being used today originated with the ideas of Emmons (1951) (see also Dhawan & Narasimha 1958). The desire to account for transitional effects in RANS models has been given renewed emphasis and has been ongoing for some time (Schmidt & Patankar 1991; Stock & Haase 1999; Walters & Leylek 2004; see also Rumsey *et al*. (2005) for a more complete list). These efforts have mainly focused on engineering approaches rather than on approaches that are more theoretically based. The latter approaches have received much less attention (e.g. Thacker *et al*. 1999*a*,*b*; Jovanović & Pashtrapanska 2004; Rumsey *et al*. 2005) and as yet have not reached a point where extensive flow field predictions have been made. It is worth pointing out some characteristics that a predictive transition-sensitized turbulence model should possess. These include properly identifying the quantities actually being computed by such models, not requiring the *a priori* input of transition or relaminarization locations, and properly characterizing the dynamics that such quantities should possess. For methods using closed transport equations for dynamic variables, such as disturbance kinetic energy or disturbance velocity second moments (stresses), the intent is to compute the entire flow domain such as upstream, over and downstream of airfoils, wings, turbine blades or aerodynamic bodies. It is necessary to formulate equations capable of describing both viscous dominated and turbulent fluctuations within the same flow.

In general, the usual types of transition modes considered are natural and by-pass transition. In natural transition, the fluid responds to a small random disturbance, which spatially and temporally evolves into a complex nonlinear process leading to a fully turbulent flow. The statistical correlations associated with these disturbances in the laminar transition region can be described by transport equations. While stability theory focuses on the description of deterministic quantities, the turbulent prediction methods such as RANS use equations for statistical correlations. Thus, these transitioning disturbance fields are represented by correlations of fluctuating quantities with an associated probability density function (Thacker *et al*. 1999*a*,*b*) and are governed by the same type of correlation transport equations as the RANS stresses in the fully turbulent region. In by-pass transition or in a relaminarizing turbulent flow, the flow dynamics are governed by either turbulent fluctuations in the free-stream (by-pass transition) or an existing turbulent field being subjected to some flow variation (such as pressure gradients) causing the flow to relaminarize.

Although the characterizations of the disturbance field just described validate the use of transport equations that are similar in form to the RANS equations, the commonly used engineering prediction methods have directly applied modified RANS models to the prediction of such transitioning flows. The most common and simple approach to handling a transition region on an aerodynamic body is to control the energy production term by having it ‘switched off’ in the transition region. Unfortunately, this *a priori* knowledge limits the range of applicability and the general predictive nature of such models. More sophisticated approaches have been employed, such as the correlation-based method of Langtry & Menter (2005), which is built strictly on local variables and solves two transport equations, one for intermittency and one for a transition onset criterion. These currently accepted engineering approaches to delimiting the transitioning regions of the flow need to be further investigated. Criterion should be developed (e.g. Jovanović & Pashtrapanska 2004) that would identify such transitioning zones based on some dynamic measure rather than a correlated input invoked by the user.

Having now identified the type of disturbance correlation fields being computed as well as the threshold criterion commonly used delimiting these transitioning zones, the next (more subtle) point is whether such correlation fields correctly describe the transitioning or relaminarization processes. There are two dynamic characteristics that the transitional and relaminarization fields should possess: the first is that the ‘disturbance’ or eddy viscosity should be vanishing and second, the ratio of the mean flow time-scale (represented by the mean shear *S*) to the disturbance field time-scale (represented by the kinetic energy *K* and kinetic energy dissipation rate *ϵ*) vanish. The former dynamic characteristic is easily met and is usually achieved in existing models by ‘turning off’ the production term in either the kinetic energy or second-moment equations in the transitioning region.

Recent studies by Rumsey *et al*. (2006*b*) and Pettersson Reif *et al*. (2006) have identified some distinguishing characteristics of several two-equation linear eddy-viscosity and algebraic stress models that can have an impact on any attempts to use such formulations in predicting transitioning and relaminarizing flows. These studies have shown through a dynamical systems analysis using nullclines in *K* and *ϵ* phase space that solution regions where disturbance (eddy) viscosities given by *K*^{2}/*ϵ*, for example, often do not predict the correct laminar limit in terms of the mean flow to turbulent time-scale ratio, which should be *ϵ*/*SK*→0 in the laminar regime. Such regions showing incorrect behaviour were termed ‘pseudo laminar’. These results have a direct implication on the use of such models for the prediction of transitioning and relaminarizing flows. Without proper calibration and evaluation, RANS-type models can yield solutions qualitatively similar to transitioning flow behaviour but with the incorrect dynamic characteristics.

As an example of this anomalous behaviour, the flow over an RAE2822 airfoil at free-stream Mach number *M*=0.75, angle of attack=2.72 and *Re*=6.2×10^{6} was computed. A plot showing the airfoil shape and resulting pressure contours for these conditions is given in figure 5. There is a strong shock wave present on the airfoil upper surface near 65% chord, whereas the flow on the lower surface remains subsonic. In this example, the initial and boundary conditions were kept fixed but the numerical solution method was changed (for details see Rumsey *et al*. 2006*b*). Figure 6 shows the skin-friction distribution on the lower surface of the airfoil obtained using two different numerical solution strategies for obtaining converged solutions. The converged results were completely different, with each suggesting a transition location in a different place. This inconsistent behaviour was due to the fact that particular forms of the *K*–*ϵ* model can converge to either a ‘pseudo-laminar’ state or the intended turbulent state, depending on initial conditions or other numerical parameters. Figure 7 shows the computed results along with the theoretical trajectories at a point in the boundary layer on a flat plate that eventually converged to a ‘pseudo-laminar’ result. Here, diffusion effects were negligible, and there is excellent agreement between theory and computation. The *SK*−*ϵ* phase plot shows that the degenerate fixed point at *SK*=*ϵ*=0 is reached.

These results indicate that any attempt at using transport equations for statistical correlations such as kinetic energy and associated dissipation rate can yield solutions that depend on the character of the equations rather than on any physical calibration. Thus, the use of such equations in predicting flows containing transitioning and relaminarizing regions can be problematic. A careful assessment of the solutions is needed to determine whether they are replicating the calibration characteristics or are a result of phase plane trajectories that are inherent in the form of the equations.

### (b) Separation and control

Another topic of current interest is flow control, particularly for the reduction or elimination of separation. By using suction, small jets or actuators on wings, for example, it is possible to delay or eliminate the onset of separation. There have recently been attempts to assess both the experimental measurement capabilities as well as the model prediction capabilities of representative benchmark flows (Rumsey *et al*. 2006*a*). In such flows, flow field characteristics that are not well understood or documented need to be investigated. These include the effects of fluid injection through synthetic jets into quiescent air and crossflow and the effects of steady and unsteady fluid injection and suction on separated flow regions. A ‘hump-model’ configuration (Seifert & Pack 2002) exemplifies the latter effect and is also a realistic configuration for flow control on wings. Some recent results can be shown to highlight the predictive capabilities of current methods. This flow has been studied extensively in both steady and time-dependent flow control applications using a wide variety of methods, which included RANS, hybrid RANS/LES, LES and (underresolved) DNS. The cases studied have included no-control, steady suction and synthetic jet (unsteady blowing/suction) control conditions.

The subsonic flow boundary layer in this case separates at the back side of the hump, near 65% chord. With no flow control, the flow reattaches near *x*/*c*=1.11. With either steady suction or synthetic jet flow controls applied near the separation point, the separation extent can be either reduced or eliminated. Most RANS models appear to do a reasonable job in predicting the separation point in these cases, but all models ranging from one-equation eddy viscosity through full Reynolds stress invariably predict too large a separation extent, either in steady state or in the mean. A typical example is shown in figure 8, using a *K*–*ω* explicit algebraic stress model. In this figure, time-averaged mean results are shown for the synthetic jet case, with oscillation frequency 138.5 Hz and peak velocity out of the slot of 26–27 m s^{−1} (free-stream velocity was 34.6 m s^{−1}). In the experiment (Greenblatt *et al*. 2005), the mean flow reattaches near *x*/*c*=1.0; whereas the RANS method predicts later reattachment. The reason for the poor RANS behaviour is believed to be that it dramatically underpredicts the turbulent shear stress in the separated region, as shown for typical results in figure 9. This underprediction leads to too little mixing, and hence a tendency to remain separated too long.

On the other hand, blended RANS–LES, LES and DNS appear to do a better job in predicting the turbulence characteristics in the separated region. By resolving the larger eddies (rather than modelling them), much of the turbulence dynamics in the separated region can be properly accounted for, yielding earlier reattachment. However, these methods are generally very time consuming to compute, and underresolution in time or space can have a significant negative impact on results. Nonetheless, recent LES results (Saric *et al*. 2006) suggest that careful application of this method is capable of capturing important effects seen in the hump experiment.

Other recent studies have focused on circulation control applications using Coanda-type jets over flaps (e.g. Lee-Rausch *et al*. 2006; Shmilovich & Yadlin 2006) and over circular trailing edges (e.g. Slomski *et al*. 2006; Swanson & Rumsey 2006). In both these types of aerodynamic flows, flow control can greatly enhance the maximum achievable lift. In the case of computing the Coanda flow over circular trailing edges, the flowfield—particularly regarding where the jet separates—has been found to be very sensitive to both numerical parameters as well as to turbulence modelling. For example, proper sensitization of the turbulence model to streamline curvature effects appears to be very important. Typical examples showing streamlines around a Coanda surface at the back of an airfoil are shown in figure 10. The experiment is due to Novak *et al*. (1987). Here, a steady jet issues out of a small slot located on the upper surface near *x*/*c*=0.93. Owing to the Coanda effect, the jet ‘sticks’ to the airfoil surface and draws the flow around to just beyond the back of the circular trailing edge, enhancing the circulation around the body. In this case, a one-equation linear eddy-viscosity turbulence model predicts that the jet separates from the Coanda surface too late, whereas the same model with a curvature correction predicts results in good agreement with experiment. However, because it is also very difficult to perform high-quality validation experiments (that provide sufficient information to evaluate and improve turbulence models) on this type of configuration, separation flow control remains an area of very active research both experimentally and computationally.

## 5. Concluding remarks

Significant strides have been made in the prediction of turbulent aerodynamic flows in the last quarter century. The availability and applicability of commercially available numerical solvers, which routinely solve very complex aerodynamic flow fields, has led some to believe that the technical discipline has reached full maturity. Understanding and ultimately predicting complex turbulent flow fields has been and still remains an unsolved problem. With each level of success come increased demands on accuracy and solution capabilities, which at its limits remain beyond current capabilities.

The topics addressed in this overview are but a brief glimpse into topical areas of research that are currently being investigated and which will help extend our ability to both understand and predict even more complex aerodynamic flows.

## Acknowledgments

The authors would like to thank A. Fadai-Ghotbi of the Laboratoire d'Etudes Aerodynamiques for providing the results used in figure 2.

## Footnotes

One contribution of 9 to a Theme Issue ‘Computational fluid dynamics in aerospace engineering’.

- © 2007 The Royal Society