Thomas B Gatski

Christopher L Rumsey

Rémi Manceau

Current Trends in Modeling Research For Turbulent Aerodynamic Flows

Keywords: turbulence, modeling, hybrid-methods, computational fluid dynamics

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 as well as 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, that is, separation control. This review gives some current priorities in traditional Reynolds-averaged modeling research as well as some methodologies being applied to a new class of turbulent flow control problems.

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 methodology for industrial prediction. Aerodynamic flows are characterized by a wide parameter range associated with, for example, Mach number and Reynolds number and are turbulent or are transitioning to or from turbulent. 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 optimized. It is clear that improving operating efficiencies are crucial in developing next generation aircraft and aerospace vehicles. For this reason, accurate and efficient prediction of turbulent flows has become an important topic.

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 for 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 Reynolds number (based on friction velocity) of O(10 3 ). Therefore, either all or part of those scales that cannot be directly computed need to modeled. The general task then is to develop 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 fluctuating part, an idea first proposed by Osborne Reynolds (1895) over a hundred years ago. This process, known as the Reynolds decomposition, assumes the instantaneous flow can be partitioned into a fluctuating part representing all the turbulent motion centered 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 completely resolve the turbulent motion, its drawback is that unknown single-point, high-order correlations appear in both the mean and 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 as well as 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 [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF], and extended a few years later by Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds number[END_REF]. 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 to it (e.g. [START_REF] Geurts | Modern Simulation Strategies for Turbulent Flow[END_REF] 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 modeled. For statistically unsteady flows, the need for alternatives to the RANS approach is more obvious since the foundation for RANS models is based on (spectral) equilibrium hypotheses that may no longer be valid (e.g., [START_REF] Carpy | Turbulence modelling of statistically periodic flows: synthetic jet into quiescent air[END_REF]. The need for such predictions has led to the development of hybrid methods capable of having a RANS-type behavior in the vicinity of a solid boundary and an LES-type behavior 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 behavior just described. The most often utilized approach to date is Article submitted to Royal Society the detached eddy simulation (DES) approach developed by Spalart and colleagues (Spalart et al.1997[START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF]. An insightful discussion of this topic, from an LES perspective, is given by [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. Hybrid methods can be essentially classified in 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, an small overlap region; or methods where the governing set of equations is smoothly transitionning from a RANS behavior to a LES behavior, 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 [START_REF] Bunge | Coupled RANS/LES methods[END_REF], although this terminology is ambiguous since both are using different models in different zones. Both kind 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 modeling remains a cornerstone: the hybrid models highly 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 is completely determined by the performance of the RANS model. Moreover, using a RANS-type model as subgrid-scale model can be favourable to account for complex phenomena such as heat transfer, effects of rotation or stratification, etc. In industrial applications, where a cost reduction is requested, the use of coarse meshes is not compatible with the hypotheses of a simple inertial behavior 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., [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF].

From a physical standpoint in RANS model development, the task is to characterize the turbulence. One obvious characterization is to correctly describe 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 amongst 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 [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] and the shear stress transport (SST) model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering aplications[END_REF], but in general purpose codes, different forms of high-Reynolds number Kε models, such as the standard Kε model [START_REF] Launder | The numerical computation of turbulent flows[END_REF], the RNG Kε model [START_REF] Yakhot | Renormalization group analysis of turbulence. 1. Basic theory[END_REF] and the Kε model of Shih et al. (1995a), often termed as "Realizable Kε" model, are still widely used. The main improvement of RNG and Shih et al. Kε model compared to the standard Kε model is the correction of the so-called stagnation point anomaly, through a modification of the ε equation (RNG model) or of 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 inexpensive cost. The (SST) model has gained increasing favor 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 free shear flows. An important feature of the SST model is the modification to the definition of the 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 and Rumsey (2001).

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

These increased demands on turbulence model performance have also highlighted inherent deficiencies in the models themselves. For example, the inability to consistently predict separation accurately or to properly account for pressure gradient effects have become apparent in trying to replicate multi-element airfoil flow fields at higher angles-of-attack [START_REF] Rumsey | Prediction of High-Lift Flows Using Turbulent Closure Models[END_REF]Rumsey & Gatski 2001). Other examples are the difficulties to correctly account for rotation effects [START_REF] Jakirlić | A new approach to modelling near-wall turbulence energy and stress dissipation[END_REF] and wall heat transfer [START_REF] Thielen | Predictions of flow and heat transfer in multiple impinging jets with an elliptic-blending second-moment closure[END_REF], 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.

The Filtering Process and Governing Equations

As detailed in the last section, 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 motions being described, and this is predicated on how the higher-order correlations are parameterized.

Article submitted to Royal Society

As is customary in all these formations, the flow variable f is decomposed into a filtered part, f , and a sub-filtered part, f ′ , given as

f = f + f ′ .
(2.1)

Generally, the filtering process can be defined as a subset of the general operation (e.g. [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF])

f (x, t) = G * f = G(x -x ′ , t -t ′ ) f (x ′ , t ′ ) dx ′ dt ′ . (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 Eq. (2.2) is given by

G(x, t) ≡ G T (x, t) = G(x) G T (t) = δ(x) 1 T H(T -t) H(t) , (2.3) such that f T (x, t) = G T * f = 1 T t t-T f (x, t ′ )dt ′ .
(2.4)

In such flows, the entire spectral range of scales is modeled so the sub-filtered part f ′ is a fluctuating quantity whose average is zero f ′ = 0, and the mean quantity f can be extracted from

E{f (x, t)} = lim T →∞ f T (x, T ) = lim T →∞ 1 T T 0 f (x, t) dt, (2.5)
Note that, since lim T →∞ G T (t) ≡ 0, 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 Eq. (2.5). As a practical matter, Reynolds averaging is evaluated as long-time averaging, where "long-time" is not infinity but rather a time sufficiently large compared to the turbulence time scale. Therefore, Reynolds averaging can be considered, within any predefined accuracy, as the convolution filter G T with a sufficiently large T . A useful property of this filtering is that the average of the product of two quantities is f g = f g + f ′ g ′ . 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 (cyclo-stationary), it is customary to use phase averaging, corresponding to a filter function given by

G T (x, t) = G(x) G T (t) (2.6) = δ(x) lim N →∞ 1 N + 1 N n=0 δ(t + nT ) , (2.7)
where T is the period of the cycle. In that case, the phase average of the sub-filtered part f ′ is zero, < f ′ >= f ′ = 0, and the filtered, or phase-averaged quantity < f > can be extracted from

< f (x, t) >= f (x, t) = G T * f = lim N →∞ 1 N + 1 N n=0 f (x, t + nT ),
(2.8)

Article submitted to Royal Society

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 cutoff 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 work [START_REF] Pruett | On Eulerian times-domain filtering for spatial large-eddy simulation[END_REF][START_REF] Pruett | The temporally filtered Navier-Stokes equations: Properties of the residual stress[END_REF] suggests it can be an effective alternative to the spatial filtering process. Causal time domain filters can be constructed that are analogues to the spatial filters. A simple example of causal filter is the top-hat filter given by Eq. (2.3) where T is the temporal filter width. Analogous to the spatial filtering approach, only a portion of the spectral range of scales is modeled with sub-filtered part f ′ . The filtered quantity f can then be extracted from Eq. (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 [START_REF] Pruett | The temporally filtered Navier-Stokes equations: Properties of the residual stress[END_REF], since the Reynolds-averaged function E{f } is the limit of the temporally-filtered function f T when the temporal filter width T 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 highspeed flow fields increases. It should be noted, however, that having a turbulent flow compressible does not mean 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 by using incompressible dynamics under adiabatic conditions and in the absence of shocks.

The filtering processes described above yields 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 modeling or parameterizations of the various terms requiring 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 [START_REF] Reynolds | On the dynamical theory of incompressible viscous fluids and the determination of the criterion[END_REF], or Favre variables [START_REF] Favre | Equations des gaz turbulents compressibles[END_REF]) is advantageous, since the equations take a more compact form and are structurally similar to their incompressible counterpart.

For a dependent variable f , the Favre mean is defined as

f = ρf ρ .
(2.9)

The instantaneous value f can then be decomposed into either the usual Reynolds

Article submitted to Royal Society averaged variables or the Favre-averaged variables

f = f + f ′ = f + f ′′ .
(2.10)

As might be expected, an extensive list of relations exist between the Reynoldsaveraged variables and the Favre variables [START_REF] Barre | Compressible, High Speed Flows[END_REF].

The equations for the mean density ρ and mean momentum ρũ i are given by

∂ρ ∂t + ∂ ∂x j (ρũ j ) = 0, (2.11) ρ Dũ i Dt = ∂(ρũ i ) ∂t + ∂ ∂x j (ũ j ρũ i ) = - ∂p ∂x i + ∂σ ij ∂x j - ∂(ρτ ij ) ∂x j , (2.12)
and the mean viscous stress tensor σ ij is

σ ij = 2µ S ij - 1 3 S kk δ ij ≃ 2µ Sij - 1 3 Skk δ ij , (2.13)
where µ is the mean molecular viscosity, and τ ij = u ′′ i u ′′ j is the Favre-averaged velocity correlation tensor. Equation (2.13) neglects contributions from µ ′ , and assumes that u i ≈ ũi . This assumed equality between the average velocities implies that the average fluctuating velocity u ′′ i is small since u iũi = u ′′ i . The equation for the mean total energy ρ Ẽ is

∂(ρ Ẽ) ∂t + ∂ ∂x j ũj ρ H = - ∂ ∂x j q j + ρE ′′ u ′′ j , (2.14)
where

Ẽ = c v T + ũi ũi 2 + u ′′ i u ′′ i 2 , (2.15)
H = Ẽ + p ρ , (2.16
)

q j = -k T ∂T ∂x j ≃ -k T ∂ T ∂x j , (2.17)
and

ρE ′′ u ′′ j = c p ρ u ′′ j T ′′ + ũi (ρτ ij -σ ij ) + ρu ′′ i u ′′ i u ′′ j 2 -σ ij u ′′ i -σ ′ ij u ′ i , (2.18)
In the approximation of (2.17), fluctuations in the thermal conductivity are neglected, and the Favre-averaged and Reynolds-averaged mean temperatures are taken as approximately equal. The equation of state in mean variables is

p = ρR T , (2.19)
or in terms of the mean total energy ρ Ẽ

p = (γ -1) ρ Ẽ - 1 2 ρ ũ2 + ṽ2 + w2 -ρK (2.20) Article submitted to Royal Society
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,

ρK = ρ τ ii 2 = ρ u ′′ i u ′′ i 2 , (2.21)
suggests a strong coupling between the mean equations and the turbulent transport equations.

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 Eq. (2.12) require closure through a specification of τ ij , and the mean energy equation Eq. (2.14) requires the closure of the turbulent heat flux ρc p u ′′ j T ′′ . At the second-moment level, closure is required for the velocity-pressure gradient correlation, tensor dissipation rate, the velocity-triple moments, and the turbulent mass flux u ′′ i = -ρ ′ u ′ /ρ. 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 [START_REF] Gatski | Constitutive Equations for Turbulent Flows[END_REF]. It is not possible in the space allotted here to discuss the modeling 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., [START_REF] Kreuzinger | Compressibility effects in the solenoidal dissipation rate equation: A priori assessment and modeling[END_REF]. It is safe to assume that the newer hybrid methodologies have not reached a level of maturity where detailed modeling 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 analyzing turbulent flows. For the Reynolds stress tensor, the anisotropy tensor b ij = (τ ij -2Kδ ij /3)/(2K), where K = τ ii /2, can be defined, and for the heat flux vector, for example, the scaled scalar flux vector u ′′ i T ′′ /(KT ′′2 ) 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.

Article submitted to Royal Society

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 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. 1995b;[START_REF] Craft | Development and application of a cubic eddy-viscosity model of turbulence[END_REF]. 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 modeling research for aerodynamic flows is on correlations associated with the turbulent velocity field in isothermal flows. This does not suggest that modeling issues associated with the scalar fluxes (e.g. heat and mass) are unimportant, rather it is a reflection on the limited scope of advanced modeling 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

u ′ i ∂p ′ ∂x j + u ′ j ∂p ′ ∂x i φij = -p ′ ∂u ′ i ∂x j + ∂u ′ j ∂x i -Πij + ∂ ∂x j (p ′ u ′ i ) + ∂ ∂x i p ′ u ′ j D p ij . (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 a rapid and slow part. The rapid part is modeled by assuming dynamic equilibrium conditions in both physical and spectral space, and the slow part is modeled 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 Eq. (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, increase 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 modeled away from solid-boundaries in the general form [START_REF] Chou | On velocity correlations and the solutions of the equations of turbulent fluctuation[END_REF][START_REF] Lumley | Computational modeling of turbulent flows[END_REF] as

Π ij = ǫA ij (b) + KM ijkl (b) ∂u k ∂x l , (3.2)
where A ij (b) and M ijkl (b) 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 ǫ is the isotropic dissipation rate. The isotropic dissipation rate is by the second law of thermodynamics always positive, and this Reynolds-averaged turbulent dissipation rate can be written as

ρǫ = 2µ [s ′ :s ′ ] - 1 3 {s ′ } 2 (3.3)
Article submitted to Royal Society where s ′ = s ′ ij = (u ′ i,j + u ′ j,i )/2, { } 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 [START_REF] Pantano | A study of compressibility effects in the high-speed turbulent shear layer using direct simulation[END_REF][START_REF] Pirozzoli | Direct numerical simulation and analysis of a spatially evolving supersonic boundary layer at M = 2.25[END_REF], and the first term can be related to the fluctuating enstrophy (e.g. Sarkar et al.1991).

In Eq. (3.2), the A ij (b) term is usually associated with the "slow" relaxation of the turbulence toward isotropy, and the M ijkl (b) 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 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

Π ij (r) = p ′ (x + r, t) ∂ ∂x j u ′′ i (x, t) + ∂ ∂x i u ′′ j (x, t) , (3.4a)
and its equivalent Fourier transform,

Π ij (r) = ι V d 3 x d 3 k ′ d 3 k ′′ e ιk ′ •(x+r) e ιk ′′ •x × p(k ′ ) k ′′ j ûi (k ′′ ) + k ′′ i ûj (k ′′ ) , (3.4b)
where r is the two-point separation distance (r → 0 for single point closures), û and p 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 pressurestrain rate correlation once an appropriate representation for the energy spectrum tensor is 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 M ijkl (b) given in Eq. (3.2) is related to the energy spectrum tensor through the simple relation

M ijpq (b) = k p k q k 2 E ij (b, k)d 3 k , (3.5)
where the energy spectrum tensor E ij (b, k) 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 four-term representation

E ij (k, b) = E(k) 4πk 2 δ ij - k i k j k 2 + E a (k) 8πk 2 k n k m k 2 b nm δ ij + k i k j k 2 + E a (k) 4πk 2 b ij -b in k n k j k 2 + k i k n k 2 b nj (3.6)
where E(k) is the isotropic energy spectral density and E a (k) is the anisotropic energy spectral density.

Article submitted to Royal Society

With the functional dependency outlined in Eq. (3.2), the proper representation for Π would be generated from the integrity basis given by the invariant combinations of b, S, and W (e.g. [START_REF] Gatski | Constitutive Equations for Turbulent Flows[END_REF]). Obviously, any representation composed of the full integrity basis would be unmanageably large. Even though the energy spectrum tensor given in Eq. (3.6) is only a linear function of the anisotropy tensor b, many representations for the rapid part of the pressure-strain rate correlation have used bases with quadratic and cubic terms in b (cf. [START_REF] Sjögren | Development and calibration of algebraic nonlinear models for terms in the Reynolds stress transport equations[END_REF] although they only retain 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 Eq. (3.2); however, 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.

(ii) Dissipation Rate Tensor

The other unknown correlation that appears in the differential Reynolds stress equations and that continues to be of modeling interest is the tensor dissipation rate ǫ ij . Even though it is possible to derive a transport equation for ρǫ ij from the fluctuating momentum equation by taking the moment (cf. Speziale 1991)

2µ ∂u ′ i ∂x k ∂ ∂x k N u ′ j + ∂u ′ j ∂x k ∂ ∂x k (N u ′ i ) = 0 , (3.7)
where N 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 modeling. 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 assumption may not be generally true [START_REF] Durbin | Local anisotropy in strained turbulence at high Reynolds numbers[END_REF]. Thus, early attempts at modeling the deviatoric part of the dissipation rate tensor focused on a coupling with the Reynolds stress anisotropy [START_REF] Hanjalić | Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence[END_REF][START_REF] Hallbäck | An algebraic model for nonisotropic turbulent dissipation rate in Reynolds stress closures[END_REF].

The quantity ρǫ represents more than a transformation of kinetic energy to internal energy since it can be further decomposed into

ρǫ = µ ∂u ′ i ∂x k ∂u ′ i ∂x k ρε + ∂ ∂x k µ ∂u ′ i u ′ k ∂x i ρǫ ′ . (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 [START_REF] Bradshaw | A note on turbulent energy dissipation in the viscous wall region[END_REF] 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,

ρT µ ≡ µ u ′ i ∂ 2 u ′ i ∂x j ∂x j , (3.9)
Article submitted to Royal Society which is not necessarily positive and requires modeling, that appears in the transport equation for the turbulent kinetic energy. It is customary to decompose this term into

ρT µ = µ ∂ 2 K ∂x j ∂x j ρD µ -ρε , (3.10)
where the first term, ρD µ is called turbulent diffusion. This decomposition has two very desirable properties: ρD µ does not require modeling, and ρε is always positive which facilitates its modeling. The variable ε is often called the homogeneous dissipation, since ρT µ = ρε = ρǫ 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 infinity 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 modeled transport equation [START_REF] Hanjalić | Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence[END_REF].

There have been recent attempts to account for dissipation rate anisotropy through a tensor dissipation rate equation (e.g., [START_REF] Oberlack | Non-isotropic dissipation in non-homogeneous turbulence[END_REF][START_REF] Speziale | Analysis and modelling of anisotropies in the dissipation rate of turbulence[END_REF], 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 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 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 at only a small increase in computational cost over the linear eddy viscosity models. The details of this representation procedure have been put forth previously (e.g. [START_REF] Gatski | Constitutive Equations for Turbulent Flows[END_REF][START_REF] Gatski | Extending the Weak-Equilibrium Condition for Algebraic Reynolds Stress Models to Rotating and Curved Flows[END_REF] 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 accounting for the effects of solid boundaries. This inherently brings into consideration the proper accounting of turbulence anisotropy into the modeling process. Although this effort has sometimes been called low Reynolds number modeling, it is more properly called near-wall modeling to distinguish it from the separate problem of transition prediction discussed in the next section.

Both the pressure-strain rate correlation and the anisotropic dissipation rate that appear in the τ ij transport equation require modification, and the pressurevelocity correlation cannot be neglected anymore. Along with these parameterizations, which focus on the stress transport equations, is the behavior of the dissipation rate equation itself in the vicinity of the wall. Overall, the issue of near-wall modeling is at least as complex as the issue of developing the 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 and co-workers 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 behavior at the wall (e.g. [START_REF] Speziale | Analysis and modelling of anisotropies in the dissipation rate of turbulence[END_REF], Jakirlić & Hanjalić 2002).

A major focus over the last decade has been on what is termed the "elliptic relaxation method" [START_REF] Durbin | Near-wall turbulence closure modeling without 'damping functions[END_REF][START_REF] Durbin | A Reynolds stress model for near-wall turbulence[END_REF][START_REF] Manceau | A new form of the elliptic relaxation equation to account for wall effects in RANS modeling[END_REF][START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier-Stokes modelling[END_REF]Manceau et al. 2002;[START_REF] Laurence | A robust formulation of the v 2 -f model[END_REF]. 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 worthwhile to briefly outline the formulation. The transport equation for τ ij can be written in the form

∂τ ij ∂t + u k ∂τ ij ∂x k = P ij + φ ij -ε ij + D t ij + ν ∂ 2 τ ij ∂x k ∂x k , (3.11)
where φ ij is not traceless (cf. Π ij ), since it includes pressure diffusion (see Eq. 3.1).

In the elliptic relaxation method, the variation of the dissipation rate anisotropy d ij (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 Eq. (3.11) can be written as

φ ij -ε ij = εKf ij - τ ij K ε , (3.12)
with the relaxation function f ij defined by

εKf ij = φ ij -2ε (d ij -b ij ) . (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 inherent in the original scaling. The system is closed through the solution of a relaxation equation for

f ij 1 -L 2 ∇ 2 f ij = 1 εK (Π ij + 2εb ij ) , (3.14)
Article submitted to Royal Society where Π ij can in principle be any of the high Reynolds forms found in the literature.

Although some authors have tested nonlinear models [START_REF] Wizman | Modeling near-wall effects in second-moment closures by elliptic relaxation[END_REF], terms linear in b ij are used in general. Moreover, despite the fact that pressure diffusion can readily be accounted for in Eq. (3.14), this term is usually considered 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εb (away from the wall the dissipation rate is assumed to be isotropic

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 proposed to handle near-wall effects was developed by [START_REF] Hanjalić | Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence[END_REF]. Figure 1 shows a comparison of results from the elliptic relaxation method (Manceau et al.2002) and the Reynolds stress model [START_REF] Hanjalić | Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence[END_REF] 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, is 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 behavior 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) 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, because of the huge additional numerical effort required compared to simpler closures (the v 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 does not only lie in the requirement of solving six additional differential equations Eq. (3.14), but also in the numerically problematic wall boundary conditions for the components f 22 , f 12 and f 23 :

f ij | w = -20 ν 2 ε 2 lim x2→0 τ ij x 4 2 , (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. [START_REF] Manceau | Elliptic blending model: A new near-wall Reynolds-stress turbulence closure[END_REF] have proposed the elliptic blending Reynoldsstress model (EB-RSM), which tries to preserve all the desirable properties of the elliptic relaxation model, in particular, those properties related 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 [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], φ ijε ij is modeled as

φ ij -ε ij = (1 -α 2 ) (φ w ij -ε w ij ) + α 2 (φ h ij -ε h ij ) , (3.16)
The purpose of the introduction of the near-wall form φ w ijε w ij is the reproduction of the asymptotic behavior of φ ijε ij , a role that is played by the boundary conditions Eq. (3.15) in the elliptic relaxation model. Here, the near-wall form of the dissipation tensor is chosen as ε w ij = τ ij ε/k, but any other choice is possible: the important point being that to satisfy the behavior of φ w ijε w ij , the model for φ w ij must be adapted to the model for ε w ij . In that case, φ ij must take the form

φ ij ∼ x2→0 φ w ij ∼ x2→0 Kε f ij | w = -5 ε K   -v ′ 2 /2 u ′ v ′ 0 u ′ v ′ v ′ 2 v ′ w ′ 0 v ′ w ′ -v ′ 2 /2   (3.17)
Contrary to boundary conditions that can be imposed in a local frame linked to the wall where they are applied, Eq. (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 the ones used in the elliptic relaxation model

α -L 2 ∇ 2 α = 1 , (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

α = 1 -exp - y L (3.19)
if L was taken as a constant. Therefore, α is a function of the ratio y/L=Wall Distance/Integral Length Scale that the model feels. Figure 2 shows isocontours of α 2 around the step corner in the case of a backstep flow at Re = 37500. 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.

Article submitted to Royal Society

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 φ w ij is as well. This requirement is automatically ensured if φ w ij is made a function of objective quantities. In order to correctly enforce the asymptotic behavior Eq. (3.17), 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

N = n ⊗ n = (n i n j ) (3.20)
is also objective, and a simple way to reproduce Eq. (3.17) is to use the model

φ w ij = -5 ε K τ N + Nτ - 1 2 [τ :N] (N + I) . (3.21)
This near-wall form is then transitioned to the high Reynolds number form

φ h ij -ε h ij = Π ij - 2 3 εδ ij . (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 behaviors is reduced from 6 to 1, as compared to 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 Fig. 3.

Another focus of near-wall research, that has and continues to be of interest, is the modeling 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 modeled as a whole, many authors have investigated separately the dissipation tensor. Provided that all the models are based on an isotropic dissipation far away from solid boundaries, the modeling challenge consists of providing a smooth transition towards an anisotropic near-wall formulation.

The decomposition Eq. (3.10) between viscous diffusion and homogeneous dissipation discussed previously is not without modeling 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 behavior of the dissipation tensor in the vicinity of the wall. In the standard decomposition, where the homogeneous dissipation tensor ε ij is defined by

ρT µ ij = µ ∂ 2 u ′ i u ′ j ∂x j ∂x j ρD µ ij -µ ∂u ′ i ∂x k ∂u ′ j ∂x k ρεij , (3.23)
the asymptotic behavior 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

d ij =   b 11 2 b 12 b 13 2 b 12 4 b 22 2 b 23 b 13 2 b 23 b 33   (3.24)
which shows that the two tensors are not proportional. Therefore, the standard way of modeling the dissipation anisotropy

d ij = f ε b ij , (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 the components 12, 23 and 22. The analysis of the transport equations of the two-point correlation tensor u ′ iA u ′ jB [START_REF] Jovanović | Statistical interpretation of the turbulent dissipation rate in wall-bounded flows[END_REF] suggests that the decomposition of viscous diffusion/homogeneous dissipation should be written as

ρT µ ij = 1 2 µ ∂ 2 u ′ i u ′ j ∂x j ∂x j D µ * ij -µ ∂u ′ i ∂x k ∂u ′ j ∂x k + 1 2 µ ∂ 2 u ′ i u ′ j ∂x j ∂x j ε * ij , (3.26)
Jakirlić & Hanjalić (2002) have pointed out that with this decomposition, it is much more justified to assume proportionality between the anisotropies

d * ij = f ε b ij , (3.
f ε = 1 - √ AE 2
, 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 better represents d * ij than d ij , although, quite surprisingly, the improvement is especially visible on the components 11 and 33.

These types of approaches, focusing either on terms involving the pressure fluctuations or on the tensor dissipation behavior, have yielded encouraging results in an effort to more accurately account for near-wall behavior, 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.

Some Current Challenges

The challenge 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 Re-laminarization 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 re-laminarizing 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 with some degree of accuracy, the flow in these transitioning regions.

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 [START_REF] Emmons | The laminar-turbulent transition in a boundary-layer -part 1[END_REF] (see also [START_REF] Dhawan | Some properties of the boundary layer during the transition from laminar to turbulent flow motion[END_REF]. The desire to account for transitional effects in RANS models has been given renewed emphasis and has been ongoing for some time [START_REF] Schmidt | Simulating boundary layer transition with low-Reynolds-number kε turbulence models: Part1 An evaluation of prediction characteristics[END_REF][START_REF] Stock | Feasibility study of e N transition prediction in Navier-Stokes methods for airfoils[END_REF][START_REF] Walters | A new model for boundary layer transition using a single-point RANS approach[END_REF]; see also Rumsey et al.2005 for a more complete list). These efforts have mainly focused on engineering approaches rather than on approaches more theoretically based. The latter approaches have received much less attention (e.g., Thacker et al.1999a, b;[START_REF] Jovanović | On the criterion for the determination transition onset and breakdown to turbulence in wall-bounded flows[END_REF]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 a priori input of transition or re-laminarization location, and properly characterizing the dynamics that such quantities should possess. For methods utilizing 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 bypass 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.1999a, 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 re-laminarizing 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 re-laminarize.

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 general predictive nature of such models. More sophisticated approaches have been employed, such as the correlation-based method of [START_REF] Langtry | Transition modeling for general CFD applications in aeronautics[END_REF], 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 investigated further. Criterion should be developed (e.g. [START_REF] Jovanović | On the criterion for the determination transition onset and breakdown to turbulence in wall-bounded flows[END_REF] 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 Article submitted to Royal Society describe the transitioning or re-laminarization processes. There are two dynamic characteristics that the transitional and re-laminarization fields should possess: the first is that the "disturbance" or eddy viscosity should be vanishing and second, that 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) and Petterson-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 utilize such formulations in predicting transitioning and re-laminarizing 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 behavior were termed "pseudo laminar." These results have a direct implication on the use of such models for the prediction of transitioning and re-laminarizing flows. Without proper calibration and evaluation, RANS type models can yield solutions qualitatively similar to transitioning flow behavior but with the incorrect dynamic characteristics.

As an example of this anomalous behavior, the flow over an RAE2822 airfoil at freestream 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 Fig. 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). 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 behavior 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 dependent on the character of the equations rather than on any physical calibration. Thus, the utilization of such equations in predicting flows containing transitioning and re-laminarizing 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 inherent in the form of the equations. Another topic of current interest is flow control, particularly for 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 been recent attempts to assess both the experimental measurement capabilities as well as the model prediction capabilities of representative benchmark flows (Rumsey et al.2006). In such flows, flow field characteristics not well understood or documented need to be investigated. These include the effects of fluid injection through synthetic jets into quiescent air and cross-flow, and the effects of steady and unsteady fluid injection and suction on separated flow regions. A "hump model" configuration [START_REF] Seifert | Active flow separation control on wall-mounted hump at high Reynolds numbers[END_REF] exemplifies the latter effect, and is also a realistic configuration for flow control on wings. Some recent results can be shown highlighting 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 that 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 control applied near the separation point, the separation extent can be either reduced or eliminated. Most RANS models appear to do a reasonable job 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 Fig. 8, using a Kω Explicit Albebraic 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 (freestream velocity was 34.6 m/s). 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 behavior is believed to be the fact that it dramatically underpredicts the turbulent shear stress in the separated region, as shown for typical results in Fig. 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 predicting the turbulence characteristics in the separated region. By resolving the larger eddies (rather than modeling 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 under-resolution in time or space can have a significant negative impact on results. Nonetheless, recent LES results (Sarić 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 and[START_REF] Shmilovich | Flow control for the systematic buildup of high lift systems[END_REF] and over circular trailing edges (e.g., [START_REF] Swanson | Numerical issues for circulation control calculations[END_REF], Slomski et al.2006). In both of 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 modeling. 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 Fig. 10. The experiment is due to Novak et al. (Novak et al.1987). Here, a steady jet issues out of a small slot located on the upper surface near x/c = 0.93. Because of 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 the jet to separate 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.

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, that 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 and still remains as yet an unsolved problem. With each level of success comes increased demands on accuracy and solution capabilities that at its limits remains 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.

Figure 1 .

 1 Figure 1. Comparison of wall proximity corrections to RANS closure for fully developed turbulent channel flow.

Figure 2 .

 2 Figure 2. Isocontours of α 2 in a backstep flow. The solid lines correspond to the levels 0.25, 0.5, 0.75 and 0.9.

Figure 3 .Figure 4 .

 34 Figure 3. Comparison of Reynolds stress anisotropies predicted by the Elliptic Relaxation Model (13 equations) and the Elliptic Blending Model (8 equation) Figure 4. Comparison of the model fεbij for the dissipation anisotropy with two possible definition dij and d * ij . A priori test in a channel flow at Reτ = 590 (Moser et al. 1999).

  component 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 either on the standard decomposition Eq. (3.23) or on the decomposition Eq. (3.26) proposed by Jakirlić & Hanjalić. The blending function is modeled as

FigureFigure 7 .

 7 Figure 5. RAE2822 airfoil solution showing static pressure contours Figure 6. Streamwise variation of skin-friction coefficient on RAE2822 airfoil lower surface for two different solution procedures

Figure 8 .

 8 Figure 8. Long-time-average streamlines for hump model with synthetic jet flow control.

Figure 9 .Figure 10 .

 910 Figure 9. Turbulent shear stress profiles at x/c = 0.8 at four times during the oscillation cycle, for hump model with synthetic jet flow control.

Article submitted to Royal Society

Craft, T. J. & Launder, B. E. 2001 Closure modelling near the two-component limit.In Closure Strategies for Modelling Turbulent and Transitional Flows (eds B. E.Article submitted to Royal Society